coupling_isban.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ###############################################################################
SUBROUTINE coupling_isba_n (DTCO, UG, U, USS, NAG, CHI, NCHI, DTI, ID, NGB, GB,         &
                            ISS, NISS, IG, NIG, IO, S, K, NK, NP, NPE, NDST, SLT,       &
                            HPROGRAM, HCOUPLING, PTSTEP,  KYEAR, KMONTH, KDAY, PTIME,   &
                            KI, KSV, KSW, PTSUN, PZENITH, PZENITH2, PZREF, PUREF, PZS,  &
                            PU, PV, PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, &
                            PDIR_SW, PSCA_SW, PSW_BANDS, PPS, PPA, PSFTQ, PSFTH, PSFTS, &
                            PSFCO2, PSFU, PSFV, PTRAD, PDIR_ALB, PSCA_ALB, PEMIS,       &
                            PTSURF, PZ0, PZ0H, PQSURF, PPEW_A_COEF, PPEW_B_COEF,        &
                            PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, HTEST    )  
!     ###############################################################################
!
!!****  *COUPLING_ISBA_n * - Driver for ISBA time step   
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!! First, all actions dependant on each patch is donbe independantly
!!     (loop on patches)
!! Second, actions common to all patches (e.g. prescription of new vegetation)
!! Third, energy fluxes are averaged
!!
!! Nota that chemical fluxes are also treated.
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      P Le Moigne 11/2004 add new diagnostics for isba
!!      A.Bogatchev 09/2005 EBA snow option
!!      P Le Moigne 09/2005 AGS modifs of L. Jarlan
!!      P Le Moigne 02/2006 z0h with snow
!!      P.Le Moigne 06/2006 seeding and irrigation
!!      B. Decharme   2008  reset the subgrid topographic effect on the forcing
!!                          PSNV allways <= PSNG
!!                          News diag
!!                          Flooding scheme and allows TRIP variables coupling
!!      A.L. Gibelin 04/2009 : Add respiration diagnostics
!!      A.L. Gibelin 04/2009 : BIOMASS and RESP_BIOMASS arrays 
!!      A.L. Gibelin 04/2009 : TAU_WOOD for NCB option 
!!      A.L. Gibelin 05/2009 : Add carbon spinup
!!      A.L. Gibelin 06/2009 : Soil carbon variables for CNT option
!!      A.L. Gibelin 07/2009 : Suppress RDK and transform GPP as a diagnostic
!!      A.L. Gibelin 07/2009 : Suppress PPST and PPSTF as outputs
!!        S.Lafont   01/2011 : add PTSTEP as arg of diag_misc
!!       B.Decharme  09/2012 : Bug in hydro_glacier calculation with ES or Crocus
!!                             New wind implicitation
!!                             New soil carbon spinup and diag
!!                             Isba budget
!!      F. Bouttier  01/2013 : Apply random perturbations for ensembles
!!      B. Decharme  04/2013 new coupling variables
!!                           Subsurface runoff if SGH (DIF option only)
!!                   07/2013 Surface / Water table depth coupling
!!      P Samuelsson 10/2014 : MEB
!!      P. LeMoigne  12/2014 EBA scheme update
!!      R. Seferian  05/2015 : Add coupling fiels to vegetation_evol call
!!-------------------------------------------------------------------
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_agri_n, ONLY : agri_np_t
USE modd_ch_isba_n, ONLY : ch_isba_t, ch_isba_np_t
USE modd_data_isba_n, ONLY : data_isba_t
USE modd_surfex_n, ONLY : isba_diag_t
USE modd_gr_biog_n, ONLY : gr_biog_t, gr_biog_np_t
USE modd_sso_n, ONLY : sso_t, sso_np_t
USE modd_sfx_grid_n, ONLY : grid_t, grid_np_t
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_t, isba_p_t, isba_pe_t, isba_nk_t, isba_np_t, isba_npe_t
!
USE modd_dst_n, ONLY : dst_np_t
USE modd_slt_n, ONLY : slt_t
!
USE modd_reprod_oper, ONLY : cimplicit_wind
!
USE modd_csts,         ONLY : xrd, xrv, xp00, xcpd, xpi, xavogadro, xmd
USE modd_co2v_par,     ONLY : xmco2, xspin_co2
!
USE modd_surf_par,     ONLY : xundef
USE modd_snow_par,     ONLY : xz0sn
USE modd_type_date_surf
!
USE modd_surf_atm,    ONLY : lnosof
!
USE modd_dst_surf
USE modd_slt_surf
USE mode_dslt_surf
USE mode_meb
!
USE modd_agri,           ONLY : lagrip
USE modd_deepsoil,       ONLY : ldeepsoil
!
#ifdef TOPD
USE modd_coupling_topd,  ONLY : lcoupl_topd, nmaskt_patch
#endif
!
USE modi_irrigation_update
USE modi_add_forecast_to_date_surf
USE modi_z0eff
USE modi_isba
USE modi_average_flux
USE modi_average_phy
USE modi_average_rad
USE modi_average_diag_isba_n
USE modi_vegetation_evol
USE modi_vegetation_update
USE modi_carbon_evol
USE modi_subscale_z0eff
USE modi_soil_albedo
USE modi_albedo
USE modi_diag_inline_isba_n
USE modi_diag_evap_cumul_isba_n
USE modi_diag_misc_isba_n
USE modi_reproj_diag_isba_n
!
USE modi_update_rad_isba_n
USE modi_deepsoil_update
USE modi_isba_sgh_update
USE modi_isba_flood_properties
USE modi_diag_cpl_esm_isba
USE modi_hydro_glacier
USE modi_isba_albedo
USE modi_carbon_spinup  
USE modi_ch_aer_dep
USE modi_abor1_sfx
USE modi_average_diag_evap_isba_n
USE modi_average_diag_misc_isba_n
USE modi_ch_bvocem_n
USE modi_soilemisno_n
USE modi_ch_dep_isba
USE modi_dslt_dep
USE modi_coupling_dst_n
USE modi_coupling_surf_topd
USE modi_isba_budget_init
USE modi_isba_budget
USE modi_unpack_diag_patch_n
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
!
TYPE(agri_np_t), INTENT(INOUT) :: NAG
TYPE(ch_isba_t), INTENT(INOUT) :: CHI
TYPE(ch_isba_np_t), INTENT(INOUT) :: NCHI
TYPE(data_isba_t), INTENT(INOUT) :: DTI
TYPE(isba_diag_t), INTENT(INOUT) :: ID
TYPE(gr_biog_np_t), INTENT(INOUT) :: NGB
TYPE(gr_biog_t), INTENT(INOUT) :: GB
TYPE(sso_t), INTENT(INOUT) :: ISS 
TYPE(sso_np_t), INTENT(INOUT) :: NISS
TYPE(grid_t), INTENT(INOUT) :: IG
TYPE(grid_np_t), INTENT(INOUT) :: NIG
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_nk_t), INTENT(INOUT) :: NK
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) ::NPE
!
TYPE(dst_np_t), INTENT(INOUT) :: NDST
TYPE(slt_t), INTENT(INOUT) :: SLT
!
 CHARACTER(LEN=6),    INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=1),    INTENT(IN)  :: HCOUPLING ! type of coupling
                                              ! 'E' : explicit
                                              ! 'I' : implicit
INTEGER,             INTENT(IN)  :: KYEAR     ! current year (UTC)
INTEGER,             INTENT(IN)  :: KMONTH    ! current month (UTC)
INTEGER,             INTENT(IN)  :: KDAY      ! current day (UTC)
REAL,                INTENT(IN)  :: PTIME     ! current time since midnight (UTC, s)
INTEGER,             INTENT(IN)  :: KI        ! number of points
INTEGER,             INTENT(IN)  :: KSV       ! number of scalars
INTEGER,             INTENT(IN)  :: KSW       ! number of short-wave spectral bands
REAL, DIMENSION(KI), INTENT(IN)  :: PTSUN     ! solar time                    (s from midnight)
REAL,                INTENT(IN)  :: PTSTEP    ! atmospheric time-step                 (s)
REAL, DIMENSION(KI), INTENT(IN)  :: PZREF     ! height of T,q forcing                 (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PUREF     ! height of wind forcing                (m)
!
REAL, DIMENSION(KI), INTENT(IN)  :: PTA       ! air temperature forcing               (K)
REAL, DIMENSION(KI), INTENT(IN)  :: PQA       ! air humidity forcing                  (kg/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PRHOA     ! air density                           (kg/m3)
REAL, DIMENSION(KI,KSV),INTENT(IN) :: PSV     ! scalar variables
!                                             ! chemistry:   first char. in HSV: '#'  (molecule/m3)
!   
 CHARACTER(LEN=6), DIMENSION(KSV),INTENT(IN):: HSV  ! name of all scalar variables!
REAL, DIMENSION(KI), INTENT(IN)  :: PU        ! zonal wind                            (m/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PV        ! meridian wind                         (m/s)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PDIR_SW ! direct  solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PSCA_SW ! diffuse solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KSW),INTENT(IN)  :: PSW_BANDS ! mean wavelength of each shortwave band (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH   ! zenithal angle at t  (radian from the vertical)
REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH2  ! zenithal angle at t+1(radian from the vertical)
REAL, DIMENSION(KI), INTENT(IN)  :: PLW       ! longwave radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI), INTENT(IN)  :: PPS       ! pressure at atmospheric model surface (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PPA       ! pressure at forcing level             (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PZS       ! atmospheric model orography           (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PCO2      ! CO2 concentration in the air          (kg_CO2/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PSNOW     ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PRAIN     ! liquid precipitation                  (kg/m2/s)
!
!
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTH     ! flux of heat                          (W/m2)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTQ     ! flux of water vapor                   (kg/m2/s)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFU      ! zonal momentum flux                   (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFV      ! meridian momentum flux                (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFCO2    ! flux of CO2 positive toward the atmosphere (m/s*kg_CO2/kg_air)
REAL, DIMENSION(KI,KSV),INTENT(OUT):: PSFTS   ! flux of scalar var.                   (kg/m2/s)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTRAD     ! radiative temperature                 (K)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PDIR_ALB! direct albedo for each spectral band  (-)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PSCA_ALB! diffuse albedo for each spectral band (-)
REAL, DIMENSION(KI), INTENT(OUT) :: PEMIS     ! emissivity                            (-)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTSURF    ! surface effective temperature         (K)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0       ! roughness length for momentum         (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0H      ! roughness length for heat             (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PQSURF    ! specific humidity at surface          (kg/kg)
!
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_A_COEF! implicit coefficients
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_B_COEF! needed if HCOUPLING='I'
REAL, DIMENSION(KI), INTENT(IN) :: PPET_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPET_B_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_B_COEF
 CHARACTER(LEN=2),    INTENT(IN) :: HTEST ! must be equal to 'OK'
!
!
!*      0.2    declarations of local variables
!
!* forcing variables
!
TYPE(isba_k_t), POINTER :: KK
TYPE(isba_p_t), POINTER :: PK
TYPE(isba_pe_t), POINTER :: PEK
TYPE(sso_t), POINTER :: ISSK
!
REAL, DIMENSION(KI)     :: ZWIND    ! lowest atmospheric level wind speed           (m/s)
REAL, DIMENSION(KI)     :: ZDIR     ! wind direction                        (rad from N clockwise)
REAL, DIMENSION(KI)     :: ZEXNA    ! Exner function at lowest atmospheric level    (-)
REAL, DIMENSION(KI)     :: ZEXNS    ! Exner function at surface                     (-)
REAL, DIMENSION(KI)     :: ZALFA    ! Wind direction                                (-)
REAL, DIMENSION(KI)     :: ZQA      ! specific humidity                             (kg/kg)
REAL, DIMENSION(KI)     :: ZCO2     ! CO2 concentration                             (kg/kg)
REAL                    :: ZSPINCO2 ! CO2 concentration                             (ppmv)
REAL, DIMENSION(KI)     :: ZPEQ_A_COEF ! specific humidity implicit
REAL, DIMENSION(KI)     :: ZPEQ_B_COEF ! coefficients (hum. in kg/kg)
!
INTEGER                 ::ISPINEND
!
! Patch outputs:
!
REAL, DIMENSION(KI,IO%NPATCH) :: ZSFTH_TILE     ! surface heat flux (W/m2)
REAL, DIMENSION(KI,IO%NPATCH) :: ZSFTQ_TILE     ! surface vapor flux (kg/m2/s)
REAL, DIMENSION(KI,IO%NPATCH) :: ZSFCO2_TILE    ! surface CO2 flux positive toward the atmosphere (m/s*kg_CO2/kg_air)
REAL, DIMENSION(KI,IO%NPATCH) :: ZSFU_TILE      ! zonal momentum flux
REAL, DIMENSION(KI,IO%NPATCH) :: ZSFV_TILE      ! meridian momentum flux
REAL, DIMENSION(KI,IO%NPATCH) :: ZTRAD_TILE     ! radiative surface temperature
REAL, DIMENSION(KI,IO%NPATCH) :: ZEMIS_TILE     ! emissivity
REAL, DIMENSION(KI,IO%NPATCH) :: ZTSURF_TILE    ! surface effective temperature
REAL, DIMENSION(KI,IO%NPATCH) :: ZZ0_TILE       ! roughness length for momentum
REAL, DIMENSION(KI,IO%NPATCH) :: ZZ0H_TILE      ! roughness length for heat
REAL, DIMENSION(KI,IO%NPATCH) :: ZQSURF_TILE    ! specific humidity at surface
REAL, DIMENSION(KI,KSW,IO%NPATCH) :: ZDIR_ALB_TILE  ! direct albedo
REAL, DIMENSION(KI,KSW,IO%NPATCH) :: ZSCA_ALB_TILE  ! diffuse albedo
REAL, DIMENSION(KI,KSV,IO%NPATCH) :: ZSFTS_TILE     ! scalar surface flux
!
REAL, DIMENSION(KI, IO%NPATCH) :: ZCPL_DRAIN     ! For the coupling with TRIP
REAL, DIMENSION(KI, IO%NPATCH) :: ZCPL_RUNOFF    ! For the coupling with TRIP
REAL, DIMENSION(KI, IO%NPATCH) :: ZCPL_EFLOOD    ! For the coupling with TRIP
REAL, DIMENSION(KI, IO%NPATCH) :: ZCPL_PFLOOD    ! For the coupling with TRIP
REAL, DIMENSION(KI, IO%NPATCH) :: ZCPL_IFLOOD    ! For the coupling with TRIP
REAL, DIMENSION(KI, IO%NPATCH) :: ZCPL_ICEFLUX
!
! for chemical computations
!
REAL, DIMENSION(KI, IO%NPATCH) :: ZSW_FORBIO
!
REAL                       :: ZCONVERTFACM0_SLT, ZCONVERTFACM0_DST
REAL                       :: ZCONVERTFACM3_SLT, ZCONVERTFACM3_DST
REAL                       :: ZCONVERTFACM6_SLT, ZCONVERTFACM6_DST
!
! dimensions and loop counters
!
INTEGER :: ISWB   ! number of spectral shortwave bands
INTEGER :: JSWB   ! loop on number of spectral shortwave bands
INTEGER :: JP ! loop on patches
INTEGER :: JSV, IDST, IMOMENT, II, IMASK, JI
INTEGER :: JLAYER, JMODE, JSV_IDX
!
! logical units
!
INTEGER :: JJ, IBEG, IEND, ISIZE
LOGICAL :: GUPDATED, GALB             ! T if VEGETATION_UPDATE has reset fields
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
! --------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('COUPLING_ISBA_N',0,zhook_handle)
IF (htest/='OK') THEN
  CALL abor1_sfx('COUPLING_ISBAN: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
! --------------------------------------------------------------------------------------
!
!*      1.     Initializations
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Allocations:
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
zsfth_tile(:,:)   = xundef
zsftq_tile(:,:)   = xundef
zsfco2_tile(:,:)   = xundef
zsfu_tile(:,:)   = xundef
zsfv_tile(:,:)   = xundef
ztrad_tile(:,:)   = xundef
zemis_tile(:,:)   = xundef
zdir_alb_tile(:,:,:) = xundef
zsca_alb_tile(:,:,:) = xundef
ztsurf_tile(:,:)   = xundef
zz0_tile(:,:)   = xundef
zz0h_tile(:,:)   = xundef
zqsurf_tile(:,:)   = xundef
!
zsfts_tile(:,:,:) = 0.
!
zcpl_drain(:,:)   = 0.0
zcpl_runoff(:,:)  = 0.0
zcpl_eflood(:,:)  = 0.0
zcpl_pflood(:,:)  = 0.0
zcpl_iflood(:,:)  = 0.0
zcpl_iceflux(:,:) = 0.0
!
zsw_forbio(:,:)   =  xundef
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Forcing Modifications:
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
zdir=0.
!
DO jj=1,SIZE(pqa) 
! specific humidity (conversion from kg/m3 to kg/kg)
!
  zqa(jj) = pqa(jj) / prhoa(jj)
  zpeq_a_coef(jj) = ppeq_a_coef(jj) / prhoa(jj)
  zpeq_b_coef(jj) = ppeq_b_coef(jj) / prhoa(jj)
!
  zco2(jj) = pco2(jj) / prhoa(jj)
!
! Other forcing variables depending on incoming forcing (argument list)JJ
!
  zexns(jj)   = (pps(jj)/xp00)**(xrd/xcpd)
  zexna(jj)   = (ppa(jj)/xp00)**(xrd/xcpd)
!
!* wind strength
!
  zwind(jj) = sqrt(pu(jj)**2+pv(jj)**2)
!
!* wind direction
!
  IF (zwind(jj)>0.)  zdir(jj)=atan2(pu(jj),pv(jj))
!
!* angle between z0eff J axis and wind direction (rad., clockwise)
!
  zalfa(jj) = zdir(jj) - iss%XZ0EFFJPDIR(jj) * xpi/180.

  IF (zalfa(jj)<-xpi) zalfa(jj) = zalfa(jj) + 2.*xpi
  IF (zalfa(jj)>=xpi) zalfa(jj) = zalfa(jj) - 2.*xpi
!
ENDDO
!
!* number of shortwave spectral bands
!
iswb = ksw
!
!* irrigation
!
IF (lagrip .AND. (io%CPHOTO=='NIT'.OR. io%CPHOTO=='NCB') ) THEN
   CALL irrigation_update(nag, npe, io%NPATCH, ptstep, kmonth, kday, ptime  )  
ENDIF
!
!* Actualization of the SGH variable (Fmu, Fsat)
!
 CALL isba_sgh_update(ig%XMESH_SIZE, io, s, k, nk, np, npe, prain )
!
!
!* Actualization of deep soil characteristics
!
IF (ldeepsoil) THEN
  DO jp = 1,io%NPATCH
    kk => nk%AL(jp)
    CALL deepsoil_update(kk%XTDEEP, kk%XGAMMAT, s%TTIME%TDATE%MONTH)
  ENDDO
ENDIF
!
!* Actualization of soil and wood carbon spinup
!
! During soil carbon spinup with ISBA-CC: 
!        (1) Atmospheric CO2 concentration fixed to Pre-industrial CO2 consentration XCO2_START
!        (2) Atmospheric CO2 concentration rampin up from XCO2_START to XCO2_END
!
IF(io%LSPINUPCARBS.OR.io%LSPINUPCARBW)THEN
!
  ispinend = io%NNBYEARSPINS-nint(io%NNBYEARSPINS*xspin_co2)
!  
  io%LAGRI_TO_GRASS = .false.
!
  IF ( io%LSPINUPCARBS .AND. (io%NNBYEARSOLD <= ispinend) ) THEN
!
   io%LAGRI_TO_GRASS = .true.
!
   zco2(:) = io%XCO2_START * 1.e-6 * xmco2 / xmd
!
  ELSEIF(io%LSPINUPCARBS .AND. (io%NNBYEARSOLD > ispinend) .AND. (io%NNBYEARSOLD <= io%NNBYEARSPINS) )THEN
!
   zspinco2 = io%XCO2_START + (io%XCO2_END-io%XCO2_START) * &
                REAL(IO%NNBYEARSOLD - ISPINEND) / REAL(io%nnbyearspins - ispinend)
!
   zco2(:) = zspinco2 * 1.e-6 * xmco2 / xmd
!
  ENDIF
!
  CALL carbon_spinup( s%TTIME, io )
!
ENDIF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Time evolution
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
s%TTIME%TIME = s%TTIME%TIME + ptstep
 CALL add_forecast_to_date_surf(s%TTIME%TDATE%YEAR,s%TTIME%TDATE%MONTH,s%TTIME%TDATE%DAY,s%TTIME%TIME)
!
! --------------------------------------------------------------------------------------
!
!*      2.     Physical evolution
!
! --------------------------------------------------------------------------------------
! Patch Dependent Calculations
! --------------------------------------------------------------------------------------
!
patch_loop: DO jp=1,io%NPATCH
  
  IF (np%AL(jp)%NSIZE_P == 0 ) cycle
!
! Pack dummy arguments for each patch:
!
#ifdef TOPD
  IF (lcoupl_topd) THEN
    nmaskt_patch(:) = 0
    nmaskt_patch(1:np%AL(jp)%NSIZE_P) = np%AL(jp)%NR_P(:)
  ENDIF
#endif
  CALL treat_patch(nk%AL(jp), np%AL(jp), npe%AL(jp), niss%AL(jp), nag%AL(jp), &
                   nig%AL(jp), nchi%AL(jp), ndst%AL(jp), id%ND%AL(jp), id%NDC%AL(jp), &
                   id%NDE%AL(jp), id%NDEC%AL(jp), id%NDM%AL(jp), ngb%AL(jp)  )
!
ENDDO patch_loop
!
! --------------------------------------------------------------------------------------
! SFX - RRM coupling update if used :
! --------------------------------------------------------------------------------------
!
IF(io%LCPL_RRM)THEN
  CALL diag_cpl_esm_isba(io, s, nk, np, ptstep, zcpl_drain, zcpl_runoff, &
                         zcpl_eflood, zcpl_pflood, zcpl_iflood, zcpl_iceflux  )  
ENDIF
!
! --------------------------------------------------------------------------------------
! Vegetation update (in case of non-interactive vegetation):
! Or
! Vegetation albedo only update (in case of interactive vegetation):
! --------------------------------------------------------------------------------------
!
gupdated=.false.
!
IF (io%LVEGUPD) THEN
  galb = .false. 
  IF (io%CPHOTO=='NIT'.OR.io%CPHOTO=='NCB') galb = .true.
  DO jp = 1,io%NPATCH
    CALL vegetation_update(dtco, dti, ig%NDIM, io, nk%AL(jp), np%AL(jp), npe%AL(jp), jp, &
                         ptstep, s%TTIME, s%XCOVER, s%LCOVER,  lagrip,                   &
                         'NAT', galb, niss%AL(jp), gupdated             )  
  ENDDO
!
ENDIF
!
IF(io%LPERTSURF.AND.gupdated) THEN
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    pek => npe%AL(jp)
    issk => niss%AL(jp)

    DO ji = 1,pk%NSIZE_P
      imask = pk%NR_P(ji)
      !
      ! random perturbation for ensembles:
      ! reset these fields to their original values, as in compute_isba_parameters
      pek%XVEG(ji) = s%XPERTVEG(imask)
      pek%XLAI(ji) = s%XPERTLAI(imask)
      pek%XCV (ji) = s%XPERTCV (imask)
      ! reapply original perturbation patterns
      IF(pek%XALBNIR(ji)/=xundef)   pek%XALBNIR(ji)   = pek%XALBNIR(ji) *( 1.+ s%XPERTALB(imask) )
      IF(pek%XALBVIS(ji)/=xundef)   pek%XALBVIS(ji)   = pek%XALBVIS(ji) *( 1.+ s%XPERTALB(imask) )
      IF(pek%XALBUV(ji)/=xundef)    pek%XALBUV (ji)   = pek%XALBUV (ji) *( 1.+ s%XPERTALB(imask) )
      IF(pek%XZ0(ji)/=xundef)       pek%XZ0(ji)       = pek%XZ0(ji)     *( 1.+ s%XPERTZ0(imask) )
      IF(issk%XZ0EFFIP(ji)/=xundef) issk%XZ0EFFIP(ji) = issk%XZ0EFFIP(ji)*( 1.+ s%XPERTZ0(imask) )
      IF(issk%XZ0EFFIM(ji)/=xundef) issk%XZ0EFFIM(ji) = issk%XZ0EFFIM(ji)*( 1.+ s%XPERTZ0(imask) )
      IF(issk%XZ0EFFJP(ji)/=xundef) issk%XZ0EFFJP(ji) = issk%XZ0EFFJP(ji)*( 1.+ s%XPERTZ0(imask) )
      IF(issk%XZ0EFFJM(ji)/=xundef) issk%XZ0EFFJM(ji) = issk%XZ0EFFJM(ji)*( 1.+ s%XPERTZ0(imask) )
    ENDDO
  ENDDO
ENDIF
!
! --------------------------------------------------------------------------------------
! Outputs for the atmospheric model or update the snow/flood fraction at time t+1
! --------------------------------------------------------------------------------------
! Grid box average fluxes/properties: Arguments and standard diagnostics at time t+1
!
 CALL average_flux(s%XPATCH, zsfth_tile, zsftq_tile, zsfts_tile, &
                   zsfco2_tile, zsfu_tile, zsfv_tile, psfth, psftq,&
                   psfts, psfco2, psfu, psfv    )  
!
!
!-------------------------------------------------------------------------------
!Physical properties see by the atmosphere in order to close the energy budget 
!between surfex and the atmosphere. All variables should be at t+1 but very 
!difficult to do. Maybe it will be done later. However, Ts is at time t+1
!-------------------------------------------------------------------------------
!   
 CALL average_phy(s%XPATCH, ztsurf_tile, zz0_tile, zz0h_tile, &
                  zqsurf_tile, puref, pzref, ptsurf, pz0, pz0h, pqsurf )
!
!-------------------------------------------------------------------------------------
!Radiative properties at time t+1 (see by the atmosphere) in order to close
!the energy budget between surfex and the atmosphere
!-------------------------------------------------------------------------------------
!
DO jp = 1,io%NPATCH
  CALL update_rad_isba_n(io, s, nk%AL(jp), np%AL(jp), npe%AL(jp), jp, pzenith2, psw_bands, &
                         zdir_alb_tile(:,:,jp), zsca_alb_tile(:,:,jp),                     &
                         zemis_tile(:,jp), pdir_sw, psca_sw  )
ENDDO
!
 CALL average_rad(s%XPATCH, zdir_alb_tile, zsca_alb_tile, zemis_tile, &
                  ztrad_tile, pdir_alb, psca_alb, s%XEMIS_NAT, s%XTSRAD_NAT  )  
!
pemis = s%XEMIS_NAT
ptrad = s%XTSRAD_NAT
!
!-------------------------------------------------------------------------------------
!
! Any additional diagnostics (stored in MODD_DIAG_ISBA_n)
!
 CALL average_diag_isba_n(id%O, id%D, id%DC, id%ND, id%NDC, np, io%NPATCH, &
                          id%O%LSURF_BUDGETC, io%LCANOPY, puref, pzref, psfco2, ptrad)
!
! Cumulated diagnostics (stored in MODD_DIAG_EVAP_ISBA_n)
!
 CALL average_diag_evap_isba_n(id%O%LSURF_BUDGETC, id%DE, id%DEC, id%NDE, id%NDEC, np,  &
                               io%NPATCH, io%LGLACIER, io%LMEB_PATCH, ptstep, prain, psnow)
!
! Miscellaneous diagnostics (stored in MODD_DIAG_MISC_ISBA_n)
!
 CALL average_diag_misc_isba_n(id%DM, id%NDM, io, np, npe)
!
!--------------------------------------------------------------------------------------
!
 CALL coupling_surf_topd(id%DE, id%DEC, id%DC, id%DM, ig, &
                         io, s, k, nk, np, npe, ug, u, hprogram, u%NDIM_FULL)
!
! --------------------------------------------------------------------------------------
! Snow/Flood fractions, albedo and emissivity update :
! --------------------------------------------------------------------------------------
!
! --------------------------------------------------------------------------------------
! Chemical fluxes :
! --------------------------------------------------------------------------------------
!
IF (chi%SVI%NBEQ>0 .AND. chi%LCH_BIO_FLUX) THEN
 CALL ch_bvocem_n(chi%SVI, ngb, gb, io, s, np, npe, zsw_forbio, prhoa, psfts)
ENDIF
!
!SOILNOX
IF (chi%LCH_NO_FLUX) THEN
  CALL soilemisno_n(gb, s, k, np, npe, pu, pv)
ENDIF
!
!==========================================================================================
!
IF (lhook) CALL dr_hook('COUPLING_ISBA_N',1,zhook_handle)
CONTAINS
!
!=======================================================================================
SUBROUTINE treat_patch(KK, PK, PEK, ISSK, AGK, GK, CHIK, DSTK, DK, DCK, DEK, DECK, DMK, GBK )
!
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_sso_n, ONLY : sso_t
USE modd_agri_n, ONLY : agri_t
USE modd_ch_isba_n, ONLY : ch_isba_t
USE modd_dst_n, ONLY : dst_t
USE modd_diag_n, ONLY : diag_t
USE modd_diag_evap_isba_n, ONLY : diag_evap_isba_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
USE modd_gr_biog_n, ONLY : gr_biog_t
!
IMPLICIT NONE
!
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(sso_t), INTENT(INOUT) :: ISSK
TYPE(agri_t), INTENT(INOUT) :: AGK
TYPE(grid_t), INTENT(INOUT) :: GK
TYPE(ch_isba_t), INTENT(INOUT) :: CHIK
TYPE(dst_t), INTENT(INOUT) :: DSTK
TYPE(diag_t), INTENT(INOUT) :: DK
TYPE(diag_t), INTENT(INOUT) :: DCK
TYPE(diag_evap_isba_t), INTENT(INOUT) :: DEK
TYPE(diag_evap_isba_t), INTENT(INOUT) :: DECK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
TYPE(gr_biog_t), INTENT(INOUT) :: GBK
!
REAL, DIMENSION(PK%NSIZE_P) :: ZP_ZREF    ! height of T,q forcing                 (m)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_UREF    ! height of wind forcing                (m)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_U       ! zonal wind                            (m/s)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_V       ! meridian wind                         (m/s)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_WIND    ! wind                                  (m/s)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_DIR     ! wind direction                        (rad from N clockwise)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_QA      ! air specific humidity forcing         (kg/kg)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_TA      ! air temperature forcing               (K)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_CO2     ! CO2 concentration in the air          (kg/kg)
REAL, DIMENSION(PK%NSIZE_P,KSV) :: ZP_SV      ! scalar concentration in the air       (kg/kg)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_ZENITH  ! zenithal angle        radian from the vertical)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PEW_A_COEF ! implicit coefficients
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PEW_B_COEF ! needed if HCOUPLING='I'
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PET_A_COEF
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PET_B_COEF
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PEQ_A_COEF
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PEQ_B_COEF
REAL, DIMENSION(PK%NSIZE_P) :: ZP_RAIN    ! liquid precipitation                  (kg/m2/s)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SNOW    ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_LW      ! longwave radiation (W/m2)
REAL, DIMENSION(PK%NSIZE_P,ISWB) :: ZP_DIR_SW  ! direct  solar radiation (W/m2)
REAL, DIMENSION(PK%NSIZE_P,ISWB) :: ZP_SCA_SW  ! diffuse solar radiation (W/m2)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PS      ! pressure at atmospheric model surface (Pa)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_PA      ! pressure at forcing level             (Pa)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_ZS      ! atmospheric model orography           (m)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFTQ    ! flux of water vapor <w'q'>            (kg.m-2.s-1)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFTH    ! flux of temperature <w'T'>            (W/m2)
REAL, DIMENSION(PK%NSIZE_P,KSV) :: ZP_SFTS    ! flux of scalar      <w'sv'>           (mkg/kg/s)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFCO2   ! flux of CO2 positive toward the atmosphere (m/s*kg_CO2/kg_air)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_USTAR   ! friction velocity                     (m/s)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFU     ! zonal momentum flux                   (pa)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFV     ! meridian momentum flux                (pa)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_TRAD    ! radiative temperature                 (K)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_TSURF   ! surface effective temperature (K)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_Z0      ! roughness length for momentum (m)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_Z0H     ! roughness length for heat     (m)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_QSURF   ! specific humidity at surface  (kg/kg)
!
!*  other forcing variables (packed for each patch)
!
REAL, DIMENSION(PK%NSIZE_P) :: ZP_RHOA    ! lowest atmospheric level air density          (kg/m3)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_EXNA    ! Exner function at lowest atmospheric level    (-)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_EXNS    ! Exner function at surface                     (-)
REAL, DIMENSION(PK%NSIZE_P) :: ZP_ALFA    ! Wind direction   (-)
!
!*  working variables (packed for each patch)
!
REAL, DIMENSION(PK%NSIZE_P)      :: ZP_ALBNIR_TVEG         ! total vegetation albedo in ir
REAL, DIMENSION(PK%NSIZE_P)      :: ZP_ALBNIR_TSOIL        ! total soil albedo in ir
REAL, DIMENSION(PK%NSIZE_P)      :: ZP_ALBVIS_TVEG         ! total vegetation albedo in vis
REAL, DIMENSION(PK%NSIZE_P)      :: ZP_ALBVIS_TSOIL        ! total soil albedo in vis
REAL, DIMENSION(PK%NSIZE_P) :: ZP_EMIS                      ! emissivity
REAL, DIMENSION(PK%NSIZE_P) :: ZP_GLOBAL_SW                 ! global incoming SW rad.
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SLOPE_COS                 ! typical slope in the grid cosine
!
REAL, DIMENSION(PK%NSIZE_P) :: ZP_Z0FLOOD  !Floodplain 
REAL, DIMENSION(PK%NSIZE_P) :: ZP_FFGNOS   !Floodplain fraction over the ground without snow
REAL, DIMENSION(PK%NSIZE_P) :: ZP_FFVNOS   !Floodplain fraction over vegetation without snow
!
REAL, DIMENSION(PK%NSIZE_P,IO%NNBIOMASS) :: ZP_RESP_BIOMASS_INST         ! instantaneous biomass respiration (kgCO2/kgair m/s)
!
!*  Aggregated coeffs for evaporative flux calculations
!
REAL, DIMENSION(PK%NSIZE_P) :: ZP_AC_AGG      ! aggregated aerodynamic resistance
REAL, DIMENSION(PK%NSIZE_P) :: ZP_HU_AGG      ! aggregated relative humidity
!
!*  For multi-energy balance
!
REAL, DIMENSION(PK%NSIZE_P) :: ZPALPHAN                     ! snow/canopy transition coefficient
REAL, DIMENSION(PK%NSIZE_P) :: ZSNOWDEPTH                   ! total snow depth
REAL, DIMENSION(PK%NSIZE_P) :: ZZ0G_WITHOUT_SNOW            ! roughness length for momentum at snow-free canopy floor
REAL, DIMENSION(PK%NSIZE_P) :: ZZ0_MEBV                     ! roughness length for momentum over MEB vegetation part of patch
REAL, DIMENSION(PK%NSIZE_P) :: ZZ0H_MEBV                    ! roughness length for heat over MEB vegetation part of path
REAL, DIMENSION(PK%NSIZE_P) :: ZZ0EFF_MEBV                  ! roughness length for momentum over MEB vegetation part of patch
REAL, DIMENSION(PK%NSIZE_P) :: ZZ0_MEBN                     ! roughness length for momentum over MEB snow part of patch
REAL, DIMENSION(PK%NSIZE_P) :: ZZ0H_MEBN                    ! roughness length for heat over MEB snow part of path
REAL, DIMENSION(PK%NSIZE_P) :: ZZ0EFF_MEBN                  ! roughness length for momentum over MEB snow part of patch
! Temporary
REAL, DIMENSION(PK%NSIZE_P) :: ZP_MEB_SCA_SW                ! diffuse incoming SW rad.
!
!*  ISBA water and energy budget
!
REAL, DIMENSION(PK%NSIZE_P) :: ZP_WG_INI
REAL, DIMENSION(PK%NSIZE_P) :: ZP_WGI_INI
REAL, DIMENSION(PK%NSIZE_P) :: ZP_WR_INI
REAL, DIMENSION(PK%NSIZE_P) :: ZP_SWE_INI
!
! miscellaneous
!
REAL, DIMENSION(PK%NSIZE_P)               :: ZP_DEEP_FLUX ! Flux at the bottom of the soil
REAL, DIMENSION(PK%NSIZE_P)               :: ZP_TDEEP_A   ! coefficient for implicitation of Tdeep
REAL, DIMENSION(PK%NSIZE_P)               :: ZIRRIG_GR    ! green roof ground irrigation rate 
!
! For multi-energy balance
LOGICAL :: GMEB  ! True if multi-energy balance should be used for the specific patch
!
INTEGER :: JJ, JI, JK
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('COUPLING_ISBA_n:TREAT_PATCH',0,zhook_handle)
!
!--------------------------------------------------------------------------------------
!
! Pack isba forcing outputs
!
IF (io%NPATCH==1) THEN
   zp_zenith(:)     = pzenith(:)
   zp_zref(:)       = pzref(:)
   zp_uref(:)       = puref(:)
   zp_wind(:)       = zwind(:)
   zp_u(:)          = pu(:)
   zp_v(:)          = pv(:)
   zp_dir(:)        = zdir(:)
   zp_qa(:)         = zqa(:)
   zp_ta(:)         = pta(:)
   zp_co2(:)        = zco2(:)
   zp_sv(:,:)       = psv(:,:)
   zp_pew_a_coef(:) = ppew_a_coef(:)
   zp_pew_b_coef(:) = ppew_b_coef(:)
   zp_pet_a_coef(:) = ppet_a_coef(:)
   zp_pet_b_coef(:) = ppet_b_coef(:)
   zp_peq_a_coef(:) = zpeq_a_coef(:)
   zp_peq_b_coef(:) = zpeq_b_coef(:)
   zp_rain(:)       = prain(:)
   zp_snow(:)       = psnow(:)
   zp_lw(:)         = plw(:)
   zp_dir_sw(:,:)   = pdir_sw(:,:)
   zp_sca_sw(:,:)   = psca_sw(:,:)
   zp_ps(:)         = pps(:)
   zp_pa(:)         = ppa(:)
   zp_zs(:)         = pzs(:)
!
   zp_rhoa(:)       = prhoa(:)
   zp_exna(:)       = zexna(:)
   zp_exns(:)       = zexns(:)
   zp_alfa(:)       = zalfa(:)
ELSE
!cdir nodep
!cdir unroll=8
  DO jj=1,pk%NSIZE_P
   ji = pk%NR_P(jj)
   zp_zenith(jj)     = pzenith(ji)
   zp_zref(jj)       = pzref(ji)
   zp_uref(jj)       = puref(ji)
   zp_wind(jj)       = zwind(ji)
   zp_u(jj)          = pu(ji)
   zp_v(jj)          = pv(ji)
   zp_dir(jj)        = zdir(ji)
   zp_qa(jj)         = zqa(ji)
   zp_ta(jj)         = pta(ji)
   zp_co2(jj)        = zco2(ji)
   zp_pew_a_coef(jj) = ppew_a_coef(ji)
   zp_pew_b_coef(jj) = ppew_b_coef(ji)
   zp_pet_a_coef(jj) = ppet_a_coef(ji)
   zp_pet_b_coef(jj) = ppet_b_coef(ji)
   zp_peq_a_coef(jj) = zpeq_a_coef(ji)
   zp_peq_b_coef(jj) = zpeq_b_coef(ji)
   zp_rain(jj)       = prain(ji)
   zp_snow(jj)       = psnow(ji)
   zp_lw(jj)         = plw(ji)
   zp_ps(jj)         = pps(ji)
   zp_pa(jj)         = ppa(ji)
   zp_zs(jj)         = pzs(ji)
!
   zp_rhoa(jj)       = prhoa(ji)
   zp_exna(jj)       = zexna(ji)
   zp_exns(jj)       = zexns(ji)
   zp_alfa(jj)       = zalfa(ji)
  ENDDO
!
  DO jk=1,ksv
!cdir nodep
!cdir unroll=8
    DO jj=1,pk%NSIZE_P
      ji=pk%NR_P(jj)
      zp_sv(jj,jk) = psv(ji,jk)
    ENDDO
  ENDDO
!
  DO jk=1,SIZE(pdir_sw,2)
!cdir nodep
!cdir unroll=8
    DO jj=1,pk%NSIZE_P
      ji=pk%NR_P(jj)
      zp_dir_sw(jj,jk) = pdir_sw(ji,jk)
      zp_sca_sw(jj,jk) = psca_sw(ji,jk)
    ENDDO
  ENDDO
!
ENDIF
!
!--------------------------------------------------------------------------------------
!
! For multi-energy balance
gmeb = io%LMEB_PATCH(jp)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Cosine of the slope typically encoutered in the grid mesh (including subgrid orography)
!  and orientation of this slope
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
zp_slope_cos(:) = 1./sqrt(1.+issk%XSSO_SLOPE(:)**2)
IF(lnosof) zp_slope_cos(:) = 1.0
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Snow fractions
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! now caculated at the initialization and at the end of the time step 
! (see update_frac_alb_emis_isban.f90) in order to close the energy budget
! between surfex and the atmosphere. This fact do not change the offline runs.
!
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! No implicitation of Tdeep
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
zp_tdeep_a = 0.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Flood properties 
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF(io%LFLOOD)THEN
  CALL isba_flood_properties(pek%XLAI, kk%XFFLOOD, kk%XFFROZEN, zp_z0flood, zp_ffgnos, zp_ffvnos)  
ELSE
  zp_z0flood = xundef
  zp_ffgnos  = 0.0
  zp_ffvnos  = 0.0
ENDIF
!
! For multi-energy balance
   IF(gmeb)THEN
     zsnowdepth(:) = sum(pek%TSNOW%WSNOW(:,:)/pek%TSNOW%RHO(:,:),2)
     zpalphan(:)  =mebpalphan(zsnowdepth,pek%XH_VEG(:))
   ENDIF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Surface Roughness lengths (m):
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!* effective roughness
!
 CALL z0eff(pek%TSNOW%SCHEME, gmeb, zp_alfa, zp_zref, zp_uref,             &
            pek%XZ0, issk%XZ0REL, pek%XPSN, zpalphan, pek%XZ0LITTER,       &
            pek%TSNOW%WSNOW(:,1), issk, kk%XFF, zp_z0flood, pk%XZ0_O_Z0H,  &
            dk%XZ0, dk%XZ0H, dk%XZ0EFF, zz0g_without_snow,                 &
            zz0_mebv, zz0h_mebv, zz0eff_mebv, zz0_mebn, zz0h_mebn, zz0eff_mebn )
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Shortwave computations for outputs (albedo for radiative scheme)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! now caculated at the initialization and at the end of the time step 
! (see update_frac_alb_emis_isban.f90) in order to close the energy budget
! between surfex and the atmosphere. This fact do not change the offline runs.
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Shortwave computations for ISBA inputs (global snow-free albedo)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! ISBA needs global incoming solar radiation: it currently does
! not distinguish between the scattered and direct components,
! or between different wavelengths.
!
!
!* Snow-free surface albedo for each wavelength
!
 CALL isba_albedo(pek, io%LTR_ML, gmeb, zp_dir_sw, zp_sca_sw,                &
                  psw_bands, iswb, kk%XALBF, kk%XFFV, kk%XFFG, zp_global_sw, &
                  zp_meb_sca_sw, zp_albnir_tveg, zp_albvis_tveg,             &
                  zp_albnir_tsoil, zp_albvis_tsoil                   )  
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Intialize computation of ISBA water and energy budget
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
 CALL isba_budget_init(id%DE%LWATER_BUDGET, io%CISBA, pek, pk%XDG, pk%XDZG, &
                      zp_wg_ini, zp_wgi_ini, zp_wr_ini, zp_swe_ini  )
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Over Natural Land Surfaces:
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
zirrig_gr(:)= 0.
!
 CALL isba(io, kk, pk, pek, gk, agk, dk, dek, dmk,                                            &
           s%TTIME, s%XPOI, s%XABC, gbk%XIACAN, gmeb, ptstep, cimplicit_wind,                 &
           zp_zref, zp_uref, zp_slope_cos, zp_ta, zp_qa, zp_exna, zp_rhoa,                    &
           zp_ps, zp_exns, zp_rain, zp_snow, zp_zenith, zp_meb_sca_sw, zp_global_sw, zp_lw,   &
           zp_wind, zp_pew_a_coef, zp_pew_b_coef, zp_pet_a_coef, zp_peq_a_coef,               &
           zp_pet_b_coef, zp_peq_b_coef, zp_albnir_tveg, zp_albvis_tveg, zp_albnir_tsoil,     &
           zp_albvis_tsoil, zpalphan, zz0g_without_snow, zz0_mebv, zz0h_mebv, zz0eff_mebv,    &
           zz0_mebn, zz0h_mebn, zz0eff_mebn, zp_tdeep_a, zp_co2, zp_ffgnos, zp_ffvnos,        &
           zp_emis, zp_ustar, zp_ac_agg, zp_hu_agg, zp_resp_biomass_inst, zp_deep_flux,       &
           zirrig_gr             )
!
zp_trad = dk%XTSRAD
dk%XLE  = pek%XLE
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Glacier : ice runoff flux (especally for Earth System Model)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF(io%LGLACIER) CALL hydro_glacier(ptstep, zp_snow, pek, dek%XICEFLUX)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Calculation of ISBA water and energy budget (and time tendencies of each reservoir)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
 CALL isba_budget(io, pk, pek, dek, id%DE%LWATER_BUDGET, ptstep, zp_wg_ini, zp_wgi_ini,  &
                 zp_wr_ini, zp_swe_ini, zp_rain, zp_snow, dk%XEVAP   )
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Evolution of soil albedo, when depending on surface soil wetness:
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (io%CALBEDO=='EVOL' .AND. io%LECOCLIMAP) THEN
  CALL soil_albedo(io%CALBEDO, kk%XWSAT(:,1),pek%XWG(:,1), kk, pek, "ALL")  
  !
  CALL albedo(io%CALBEDO, pek )  
END IF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Vegetation evolution for interactive LAI
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (io%CPHOTO=='NIT' .OR. io%CPHOTO=='NCB') THEN
  CALL vegetation_evol(io, dti, pk, pek, lagrip, ptstep, kmonth, kday, ptime, gk%XLAT, &
                       zp_rhoa, zp_co2, issk, zp_resp_biomass_inst,  &
                       ! add optional for accurate dependency to nitrogen
                       ! limitation
                        pswdir=zp_global_sw ) 
END IF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Diagnostic of respiration carbon fluxes and soil carbon evolution
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!ii
zp_sfco2(:) = 0.
dek%XRESP_ECO (:) = 0.
dek%XRESP_AUTO(:) = 0.
!
IF ( io%CPHOTO/='NON' .AND. io%CRESPSL/='NON' .AND. any(pek%XLAI(:)/=xundef) ) THEN
  CALL carbon_evol(io, kk, pk, pek, dek, ptstep, zp_rhoa, zp_resp_biomass_inst )  
  ! calculation of vegetation CO2 flux
  ! Positive toward the atmosphere
  zp_sfco2(:) = dek%XRESP_ECO(:) - dek%XGPP(:)  
END IF
!
IF ( io%CPHOTO/='NON') THEN
  dek%XGPP(:) = dek%XGPP(:) * zp_rhoa(:)
  dek%XRESP_ECO(:) = dek%XRESP_ECO(:) * zp_rhoa(:)
  dek%XRESP_AUTO(:) = dek%XRESP_AUTO(:) * zp_rhoa(:)
ENDIF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Reset effecitve roughness lentgh to its nominal value when snow has just disappeared
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
 CALL subscale_z0eff(issk,pek%XZ0(:),.false.,omask=(pek%TSNOW%WSNOW(:,1)==0. .AND. pek%XPSN(:)>0.)  )   
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Turbulent fluxes
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
zp_sfth(:) = dk%XH(:)
zp_sftq(:) = dk%XEVAP(:)

zp_sfu(:) = 0.
zp_sfv(:) = 0.
WHERE (zp_wind>0.)
  zp_sfu(:) = - zp_u(:)/zp_wind(:) * zp_ustar(:)**2 * zp_rhoa(:)
  zp_sfv(:) = - zp_v(:)/zp_wind(:) * zp_ustar(:)**2 * zp_rhoa(:)
END WHERE
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Scalar fluxes
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
zp_sfts(:,:) = 0.
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! --------------------------------------------------------------------------------------
! Chemical dry deposition :
! --------------------------------------------------------------------------------------
IF (chi%SVI%NBEQ>0) THEN
  IF( chi%CCH_DRY_DEP == "WES89") THEN

    ibeg = chi%SVI%NSV_CHSBEG
    iend = chi%SVI%NSV_CHSEND 
    isize = iend - ibeg + 1 

    CALL ch_dep_isba(kk, pk, pek, dk, dmk, chik, &
                     zp_ustar, zp_ta, zp_pa, zp_trad(:), isize )  
 
    zp_sfts(:,ibeg:iend) = - zp_sv(:,ibeg:iend) * chik%XDEP(:,1:chi%SVI%NBEQ)  

    IF (chi%SVI%NAEREQ > 0 ) THEN
       
      ibeg = chi%SVI%NSV_AERBEG
      iend = chi%SVI%NSV_AEREND
      CALL ch_aer_dep(zp_sv(:,ibeg:iend), zp_sfts(:,ibeg:iend), zp_ustar, pek%XRESA, zp_ta, zp_rhoa)     
    END IF
  ELSE

    ibeg = chi%SVI%NSV_AERBEG
    iend = chi%SVI%NSV_AEREND
    zp_sfts(:,ibeg:iend) = 0.
    zp_sfts(:,ibeg:iend) = 0.

  ENDIF
ENDIF
!
! --------------------------------------------------------------------------------------
! Dust deposition and emission:
! --------------------------------------------------------------------------------------
!
IF(chi%SVI%NDSTEQ>0)THEN

  ibeg = chi%SVI%NSV_DSTBEG
  iend = chi%SVI%NSV_DSTEND
  idst = iend - ibeg + 1

  CALL coupling_dst_n(dstk, kk, pk, pek, dk, &
            hprogram,                    &!I [char] Name of program
            pk%NSIZE_P,      &!I [nbr] number of points in patch
            idst,                        &!I [nbr] number of dust emissions variables
            zp_ps,                       &!I [Pa] surface pressure
            zp_qa,                       &!I [kg/kg] specific humidity
            zp_rhoa,                     &!I [kg/m3] atmospheric density
            zp_pa,                       &!I [K] Atmospheric pressure
            zp_ta,                       &!I [K] Atmospheric temperature
            zp_u,                        &!I [m/s] zonal wind at atmospheric height 
            zp_uref,                     &!I [m] reference height of wind
            zp_v,                        &!I [m/s] meridional wind at atmospheric height
            zp_zref,                     &!I [m] reference height of wind
            zp_sfts(:,ibeg:iend)  &!O [kg/m2/sec] flux of dust            
            )  
!
   IF (chi%SVI%NSV_AEREND > 0)  THEN ! case of dust/ anthropogenic aerosols coupling

     DO jmode=1,ndstmde
       !
       !Make index which is 0 for first mode, 3 for second, 6 for third etc
       IF (lvarsig_dst) THEN
         jsv_idx = (jmode-1)*3
       ELSE IF (lrgfix_dst) THEN
         jsv_idx = jmode-2
       ELSE
         jsv_idx = (jmode-1)*2
       END IF
       !
       DO jsv=1, size(hsv)
         IF ((trim(hsv(jsv)) == "@DSTI").AND.(jmode==3)) THEN 
           ! add dust flux and conversion kg/m2/s into molec.m2/s
           zp_sfts(:,jsv) = zp_sfts(:,jsv) + zp_sfts(:,ibeg-1+jsv_idx+2)*xavogadro/xmolarweight_dst
         END IF
         IF ( (trim(hsv(jsv)) == "@DSTJ").AND.(jmode==2)) THEN 
           ! add dust flux and conversion kg/m2/sec into molec.m2/s
           zp_sfts(:,jsv) = zp_sfts(:,jsv) + zp_sfts(:,ibeg-1+jsv_idx+2)*xavogadro/xmolarweight_dst
         END IF
       END DO
       !
     END DO
    END IF
!    
!Modify fluxes due to dry deposition, we introduce a negative flux where dust is lost
  CALL dslt_dep(zp_sv(:,ibeg:iend), zp_sfts(:,ibeg:iend), zp_ustar, pek%XRESA,        &
                zp_ta, zp_rhoa, dstk%XEMISSIG_DST, dstk%XEMISRADIUS_DST, jpmode_dst,  &
                xdensity_dst, xmolarweight_dst, zconvertfacm0_dst, zconvertfacm6_dst, &
                zconvertfacm3_dst, lvarsig_dst, lrgfix_dst, cvermod  )
!
!Transfer these fluxes to fluxes understandable by all moments
  CALL massflux2momentflux(           &
    zp_sfts(:,ibeg:iend),             & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
    zp_rhoa,                          & !I [kg/m3] air density
    dstk%XEMISRADIUS_DST,              & !I [um] emitted radius for the modes (max 3)
    dstk%XEMISSIG_DST,                 & !I [-] emitted sigma for the different modes (max 3)
    ndstmde,                          &
    zconvertfacm0_dst,                &
    zconvertfacm6_dst,                &
    zconvertfacm3_dst,                &
    lvarsig_dst, lrgfix_dst           )   
!
ENDIF !Check on CDSTYN
!
! --------------------------------------------------------------------------------------
! Sea Salt deposition
! --------------------------------------------------------------------------------------
!
IF (chi%SVI%NSLTEQ>0) THEN
  !
  ibeg = chi%SVI%NSV_SLTBEG
  iend = chi%SVI%NSV_SLTEND
  !
  CALL dslt_dep(zp_sv(:,ibeg:iend), zp_sfts(:,ibeg:iend), zp_ustar, pek%XRESA,         &
                zp_ta, zp_rhoa, slt%XEMISSIG_SLT, slt%XEMISRADIUS_SLT, jpmode_slt,     &
                xdensity_slt, xmolarweight_slt, zconvertfacm0_slt, zconvertfacm6_slt,  &
                zconvertfacm3_slt, lvarsig_slt, lrgfix_slt, cvermod  )  

  CALL massflux2momentflux(           &
    zp_sfts(:,ibeg:iend),             & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
    zp_rhoa,                          & !I [kg/m3] air density
    slt%XEMISRADIUS_SLT,              & !I [um] emitted radius for the modes (max 3)
    slt%XEMISSIG_SLT,                 & !I [-] emitted sigma for the different modes (max 3)
    nsltmde,                          &
    zconvertfacm0_slt,                &
    zconvertfacm6_slt,                &
    zconvertfacm3_slt,                &
    lvarsig_slt, lrgfix_slt         ) 
ENDIF !Check on CSLTYN
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Inline diagnostics
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
 CALL diag_inline_isba_n(id%O, kk, dk, io%LCANOPY, zp_ta, zp_qa, zp_pa,         &
                         zp_ps, zp_rhoa, zp_u, zp_v, zp_zref, zp_uref, zp_sfth, &
                         zp_sftq, zp_sfu, zp_sfv, zp_dir_sw, zp_sca_sw, zp_lw )  
!
!
!-------------------------------------------------------------------------------
!Physical properties see by the atmosphere in order to close the energy budget 
!between surfex and the atmosphere. All variables should be at t+1 but very 
!difficult to do. Maybe it will be done later. However, Ts can be at time t+1
!-------------------------------------------------------------------------------
!
zp_tsurf(:) = dk%XTS (:)
zp_z0(:) = dk%XZ0 (:)
zp_z0h(:) = dk%XZ0H(:)
zp_qsurf(:) = dk%XQS (:)
!
!-------------------------------------------------------------------------------
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Isba offline diagnostics for each patch
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
 CALL diag_evap_cumul_isba_n(id%O%LSURF_BUDGETC, id%DE, deck, dck, dek, dk, pek, &
                             io, ptstep, pk%NSIZE_P, jp, zp_rhoa)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Isba offline diagnostics for miscellaneous terms over each patch
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
 CALL diag_misc_isba_n(dmk, kk, pk, pek, agk, io, id%DM%LSURF_MISC_BUDGET, &
                       id%DM%LVOLUMETRIC_SNOWLIQ, ptstep, lagrip, ptime, pk%NSIZE_P )                  
!
 CALL reproj_diag_isba_n(dk, dek, dmk, pek, id%O%LSURF_BUDGET, id%DE%LSURF_EVAP_BUDGET, &
                         id%DE%LWATER_BUDGET, id%DM%LSURF_MISC_BUDGET, id%DM%LPROSNOW, &
                         io%LMEB_PATCH(jp), zp_slope_cos)
!
! Unpack ISBA diagnostics (modd_diag_isban) for each patch:ISIZE_MAX = MAXVAL(NSIZE_NATURE_P)

!  (MUST be done BEFORE UNPACK_ISBA_PATCH, because of XP_LE)
!
IF (pek%TSNOW%SCHEME=='3-L'.OR.pek%TSNOW%SCHEME=='CRO') THEN
  pek%TSNOW%TEMP(:,:) = dmk%XSNOWTEMP(:,:)
  pek%TSNOW%TS  (:)   = dmk%XSNOWTEMP(:,1)
ENDIF
!

 CALL unpack_diag_patch_n(io, dek, pk, pk%NR_P, pk%NSIZE_P, io%NPATCH, jp,    &
                          zcpl_drain, zcpl_runoff, zcpl_eflood, zcpl_pflood,  &
                          zcpl_iflood, zcpl_iceflux)  
!
!----------------------------------------------------------------------
!
! for further chemical biogenic emissions
!
IF (chi%SVI%NBEQ>0 .AND. chi%LCH_BIO_FLUX) THEN
  !
  DO jj=1,pk%NSIZE_P
    zsw_forbio(pk%NR_P(jj),jp) = 0.
  ENDDO
  !
  DO jswb=1,iswb
!cdir nodep
!cdir unroll=8
    DO jj=1,pk%NSIZE_P
      zsw_forbio(pk%NR_P(jj),jp) = zsw_forbio(pk%NR_P(jj),jp) + zp_dir_sw(jj,jswb) + zp_sca_sw(jj,jswb)  
    ENDDO
  ENDDO
  !
ENDIF
!----------------------------------------------------------------------
!
! Unpack output dummy arguments for each patch:
!
IF (io%NPATCH==1) THEN
   zsftq_tile(:,jp)  = zp_sftq(:)
   zsfth_tile(:,jp)  = zp_sfth(:)
   zsfts_tile(:,:,jp)= zp_sfts(:,:)
   zsfco2_tile(:,jp)  = zp_sfco2(:)
   zsfu_tile(:,jp)  = zp_sfu(:)
   zsfv_tile(:,jp)  = zp_sfv(:)
   ztrad_tile(:,jp)  = zp_trad(:)
   ztsurf_tile(:,jp)  = zp_tsurf(:)
   zz0_tile(:,jp)  = zp_z0(:)
   zz0h_tile(:,jp)  = zp_z0h(:)
   zqsurf_tile(:,jp)  = zp_qsurf(:)   
ELSE
!cdir nodep
!cdir unroll=8
 DO jj=1,pk%NSIZE_P
   ji = pk%NR_P(jj)
   zsftq_tile(ji,jp)  = zp_sftq(jj)
   zsfth_tile(ji,jp)  = zp_sfth(jj)
   zsfco2_tile(ji,jp)  = zp_sfco2(jj)
   zsfu_tile(ji,jp)  = zp_sfu(jj)
   zsfv_tile(ji,jp)  = zp_sfv(jj)
   ztrad_tile(ji,jp)  = zp_trad(jj)
   ztsurf_tile(ji,jp)  = zp_tsurf(jj)
   zz0_tile(ji,jp)  = zp_z0(jj)
   zz0h_tile(ji,jp)  = zp_z0h(jj)
   zqsurf_tile(ji,jp)  = zp_qsurf(jj)   
 ENDDO
!
!cdir nodep
!cdir unroll=8
  DO jk=1,SIZE(zp_sfts,2)
    DO jj=1,pk%NSIZE_P
      ji=pk%NR_P(jj)    
      zsfts_tile(ji,jk,jp)= zp_sfts(jj,jk)
    ENDDO
  ENDDO
ENDIF
!
!----------------------------------------------------------------------
!
! Get output dust flux if we are calculating dust
IF (ndstmde .GE. 1) imoment = int(idst / ndstmde)
IF (chi%SVI%NDSTEQ>0) THEN
  DO jsv = 1,ndstmde
    IF (imoment == 1) THEN
      dstk%XSFDST(:,jsv) = zsfts_tile(:,ndst_mdebeg+jsv-1,jp)
    ELSE
      dstk%XSFDST(:,jsv) = zsfts_tile(:,ndst_mdebeg+(jsv-1)*imoment+1,jp)
    END IF

    dstk%XSFDSTM(:,jsv) = dstk%XSFDSTM(:,jsv) + dstk%XSFDST(:,jsv) * ptstep
  ENDDO
ENDIF
!
IF (lhook) CALL dr_hook('COUPLING_ISBA_n:TREAT_PATCH',1,zhook_handle)
!
END SUBROUTINE treat_patch
!==========================================================================================
END SUBROUTINE coupling_isba_n