coupling_flaken.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_flake_n (CHF, DGO, D, DC, DMF, F, DST, SLT, &
                             HPROGRAM, HCOUPLING, PTSTEP, KYEAR, KMONTH, KDAY, PTIME,  &
                             KI, KSV, KSW, PTSUN, PZENITH, PZENITH2, PAZIM, PZREF,     &
                             PUREF, 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, PTSRAD, 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_FLAKE_n * - Driver for FLAKE scheme for lakes
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      V. Masson   05/2009 Implicitation of momentum fluxes
!!      B. Decharme 01/2010 Add XTT in water_flux
!!      V. Masson   11/2011 Ch limited to 1.E-7 in all cases and Cd coming from
!!                          Flake_interface routine if computed by flake
!!      B. Decharme 09/2012 New wind implicitation
!!      P. Le Moigne 10/2012 ECUME option for FLake. Remove wind threshold
!!      P. Le Moigne 04/2013 Remove ECUME option for FLake
!!      P. Le Moigne 04/2013 Chemistry, UPDATE_RAD_FLAKE
!!      B. Decharme  04/2013 New diag, new coupling variables
!!      P. Le Moigne 10/2014 Threshold on Cd when fluxes computed by FLake
!!------------------------------------------------------------------------------
!
!
USE modd_ch_flake_n, ONLY : ch_flake_t
USE modd_diag_n, ONLY : diag_options_t, diag_t
USE modd_diag_misc_flake_n, ONLY : diag_misc_flake_t
USE modd_flake_n, ONLY : flake_t
!
USE modd_dst_n, ONLY : dst_t
USE modd_slt_n, ONLY : slt_t
!
USE modd_reprod_oper, ONLY : cimplicit_wind
!
USE modd_csts,     ONLY : xrd, xcpd, xp00, xlvtt, xlstt, xkarman, xtt
USE modd_surf_par, ONLY : xundef
!
!                          
!
USE modd_slt_surf
USE modd_dst_surf
!
USE mode_dslt_surf
USE mode_thermos
! 
USE modi_water_flux
USE modi_add_forecast_to_date_surf
USE modi_diag_inline_flake_n 
USE modi_diag_misc_flake_n
USE modi_ch_aer_dep
USE modi_ch_dep_water
USE modi_dslt_dep
USE modi_flake_albedo
USE modi_update_rad_flake
USE modi_abor1_sfx
USE modi_flake_interface
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(ch_flake_t), INTENT(INOUT) :: CHF
TYPE(diag_options_t), INTENT(INOUT) :: DGO
TYPE(diag_t), INTENT(INOUT) :: D
TYPE(diag_t), INTENT(INOUT) :: DC
TYPE(diag_misc_flake_t), INTENT(INOUT) :: DMF
TYPE(flake_t), INTENT(INOUT) :: F
!
TYPE(dst_t), INTENT(INOUT) :: DST
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)  :: PAZIM     ! azimuthal angle      (radian from North, clockwise)
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)  :: PCO2      ! CO2 concentration in the air          (kg/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                           (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) :: PTSRAD    ! 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   (m2s/kg)
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_B_COEF! needed if HCOUPLING='I' (m/s)
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
!
REAL, DIMENSION(KI,KSW) :: ZDIR_ALB   ! Direct albedo at time t , 
REAL, DIMENSION(KI,KSW) :: ZSCA_ALB   ! Diffuse albedo at time t
!
REAL, DIMENSION(KI)  :: ZALB     ! surface albedo
REAL, DIMENSION(KI)  :: ZSWE     ! snow water equivalent (kg.m-2)
!
REAL, DIMENSION(KI)  :: ZEXNA  ! Exner function at forcing level
REAL, DIMENSION(KI)  :: ZEXNS  ! Exner function at surface level
!
REAL, DIMENSION(KI)  :: ZWIND  ! Wind
REAL, DIMENSION(KI)  :: ZGLOBAL_SW    ! Solar radiation flux at the surface (W/m2) 
REAL, DIMENSION(KI)  :: ZQA    ! Air specific humidity (kg/kg)
!
REAL, DIMENSION(KI)  :: ZUSTAR ! friction velocity (m/s)
REAL, DIMENSION(KI)  :: ZUSTAR2! square of friction velocity (m2/s2)
REAL, DIMENSION(KI)  :: ZSFM   ! flux of momentum (Pa)
!
REAL, DIMENSION(KI)  :: ZRESA_WATER ! aerodynamical resistance
!
!salgado only for inline diagnostics - not used for the moment
!                                      flake don't have it
REAL, DIMENSION(KI)  :: ZCD    ! Drag coefficient
REAL, DIMENSION(KI)  :: ZCDN   ! Neutral Drag coefficient
REAL, DIMENSION(KI)  :: ZCH    ! Heat transfer coefficient
REAL, DIMENSION(KI)  :: ZCE    ! Heat transfer coefficient
REAL, DIMENSION(KI)  :: ZRI    ! Richardson number
REAL, DIMENSION(KI)  :: ZHU    ! Near surface relative humidity
REAL, DIMENSION(KI)  :: ZZ0H   ! heat roughness length
REAL, DIMENSION(KI)  :: ZQSAT  ! humidity at saturation
REAL, DIMENSION(KI)  :: ZTSTEP ! time-step
REAL, DIMENSION(KI)  :: ZLE    ! total latent heat flux (W/m2)
REAL, DIMENSION(KI)  :: ZLEI   ! sublimation heat flux (W/m2)
REAL, DIMENSION(KI)  :: ZSUBL  ! sublimation (kg/m2/s)
REAL, DIMENSION(KI)  :: ZLWUP  ! upward longwave flux at t
REAL, DIMENSION(KI)  :: ZTRAD  ! Radiative temperature at time t
REAL, DIMENSION(KI)  :: ZWORK  ! Work array
!
REAL                 :: ZCONVERTFACM0_SLT, ZCONVERTFACM0_DST
REAL                 :: ZCONVERTFACM3_SLT, ZCONVERTFACM3_DST
REAL                 :: ZCONVERTFACM6_SLT, ZCONVERTFACM6_DST
!
INTEGER              :: ISWB   ! number of shortwave spectral bands
INTEGER              :: JSWB   ! loop counter on shortwave spectral bands
!  
INTEGER              :: ILUOUT ! output logical unit
!
LOGICAL              :: GPWG        = .false.
LOGICAL              :: GHANDLE_SIC = .false. ! no sea-ice model
!
REAL                 :: ZEPS = 1.e-7
!
INTEGER :: IBEG, IEND
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------------
! Preliminaries:
!-------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('COUPLING_FLAKE_N',0,zhook_handle)
IF (htest/='OK') THEN
  CALL abor1_sfx('COUPLING_FLAKEN: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
!-------------------------------------------------------------------------------------
! Variables needed by flake:
!-------------------------------------------------------------------------------------
!
zcd(:) = xundef    
zcdn(:) = xundef
zch(:) = xundef
zri(:) = xundef
zresa_water(:) = xundef
zz0h(:) = xundef
zqsat(:) = xundef
zwork(:) = xundef
zalb(:) = xundef
zglobal_sw(:) = xundef
zswe(:) = xundef
!
zdir_alb(:,:) = xundef
zsca_alb(:,:) = xundef
zle(:)   = xundef
zlei(:)   = xundef
zsubl(:)   = xundef
zlwup(:)   = xundef
ztrad(:)   = xundef
!
ztstep(:) = ptstep
!
zexns(:)     = (pps(:)/xp00)**(xrd/xcpd)
zexna(:)     = (ppa(:)/xp00)**(xrd/xcpd)
!
!
zwind(:) = sqrt(pu(:)**2+pv(:)**2)
!
zqa(:) = pqa(:)/prhoa(:)
!
psfts(:,:) = 0.
!
zhu(:) = 1.
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Time evolution
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
f%TTIME%TIME = f%TTIME%TIME + ptstep
 CALL add_forecast_to_date_surf(f%TTIME%TDATE%YEAR,f%TTIME%TDATE%MONTH,f%TTIME%TDATE%DAY,f%TTIME%TIME)
!
!----------------------------------------
!
psfu(:) = 0.
psfv(:) = 0.
zsfm(:) = 0.
!
IF (f%CFLK_FLUX=='DEF  ') THEN
!
    CALL water_flux(f%XZ0, pta, zexna, prhoa, f%XTS, zexns, pqa, &
                    prain, psnow, xtt, zwind, pzref, puref,      &
                    pps, ghandle_sic, zqsat, psfth, psftq, zustar,&
                    zcd, zcdn, zch, zri, zresa_water, zz0h        )  
!                  
    WHERE (f%XTS(:)<xtt)
      zle(:) = psftq(:) * xlstt
      zlei(:) = psftq(:) * xlstt
      zsubl(:) = psftq(:)
    ELSEWHERE
      zle(:) = psftq(:) * xlvtt
      zlei(:) = 0.0
      zsubl(:) = 0.0
    END WHERE
!
    IF(cimplicit_wind=='OLD')THEN    
!     old implicitation (m2/s2)
      zustar2(:) = (zcd(:)*zwind(:)*ppew_b_coef(:))/            &
                   (1.0-prhoa(:)*zcd(:)*zwind(:)*ppew_a_coef(:)) 
    ELSE
!     new implicitation (m2/s2)            
      zustar2(:) = (zcd(:)*zwind(:)*(2.*ppew_b_coef(:)-zwind(:))) /&
                   (1.0-2.0*prhoa(:)*zcd(:)*zwind(:)*ppew_a_coef(:)) 
!                   
      zwork(:) = prhoa(:)*ppew_a_coef(:)*zustar2(:) + ppew_b_coef(:)
      zwork(:) = max(zwork(:),0.)
!
      WHERE(ppew_a_coef(:)/= 0.)
            zustar2(:) = max( ( zwork(:) - ppew_b_coef(:) ) / (prhoa(:)*ppew_a_coef(:)), 0.)
      ENDWHERE
!                   
    ENDIF
!    
    WHERE (zwind(:)>0.)            
      zsfm(:) = - prhoa(:) * zustar2(:)
      psfu(:) = zsfm(:) * pu(:) / zwind(:)
      psfv(:) = zsfm(:) * pv(:) / zwind(:)
    END WHERE
!
ELSE
   zustar(:) = f%XUSTAR(:)
   zz0h(:) = f%XZ0   (:)
ENDIF
!
!----------------------------------------
!radiative properties at t
!----------------------------------------
!
iswb = SIZE(psw_bands)
!
DO jswb=1,iswb 
  zdir_alb(:,jswb) = f%XDIR_ALB(:)
  zsca_alb(:,jswb) = f%XSCA_ALB(:)
END DO
!
ztrad = f%XTS
!
 CALL flake_albedo(pdir_sw,psca_sw,ksw,zdir_alb,zsca_alb,zglobal_sw,zalb)
!
 CALL flake_interface(f, ki, &
! Atmospheric forcing:
                      psnow, zglobal_sw, plw, puref, pzref, zwind, pta, zqa, pps, &
! Constant parameters
                      ztstep,      &
! surface albedo
                      zalb,        &
! Surface heat, momentum fluxes, and other diags
                      psfth, zle, zsfm, zz0h, zqsat, zri, zustar,     &
                      zcd, psftq, zlei, zsubl, zlwup, zswe,           &
! Flags              
                      ppew_a_coef, ppew_b_coef, prhoa, cimplicit_wind )
!
!-------------------------------------------------------------------------------------
!
! Momentum fluxes
!
IF (f%CFLK_FLUX=='FLAKE') THEN
   psfu = 0.
   psfv = 0.
  WHERE (zwind(:)>0.)
    psfu(:) = zsfm(:) * pu(:) / zwind(:)
    psfv(:) = zsfm(:) * pv(:) / zwind(:)
  END WHERE
  ! 
  ! ZUSTAR and ZRESA_WATER are not in Flake but are needed to the ch_* routines
  !
  zustar(:)       = sqrt(abs(zsfm(:))/prhoa(:))
  zexns(:)       = (pps(:)/xp00)**(xrd/xcpd)
  zexna(:)       = (ppa(:)/xp00)**(xrd/xcpd)
  zresa_water=2.e4
  WHERE (psfth/=0.)
  zresa_water(:) = xcpd * prhoa(:) * (f%XTS(:) - pta(:) * zexns(:)/zexna(:)) &
                     / (psfth(:) * zexns(:))  
  END WHERE
!
ENDIF
!                               
f%XUSTAR(:) = zustar(:)
!
! CO2 flux
!
psfco2(:)       =  0.0    ! Assumes no CO2 emission over water bodies
!
!-------------------------------------------------------------------------------------
! Scalar fluxes:
!-------------------------------------------------------------------------------------
!
!
!salgado The scalar fluxes are computed as in watflux
IF (chf%SVF%NBEQ>0) THEN
  !
  IF (chf%CCH_DRY_DEP == "WES89") THEN
    !
    ibeg = chf%SVF%NSV_CHSBEG
    iend = chf%SVF%NSV_CHSEND
    !
    CALL ch_dep_water  (zresa_water, zustar, pta, ztrad, psv(:,ibeg:iend), &
                        chf%SVF%CSV(ibeg:iend), chf%XDEP(:,1:chf%SVF%NBEQ) )  

    psfts(:,ibeg:iend) = - psv(:,ibeg:iend) * chf%XDEP(:,1:chf%SVF%NBEQ) 
    ! 
    IF (chf%SVF%NAEREQ > 0 ) THEN
      ibeg = chf%SVF%NSV_AERBEG
      iend = chf%SVF%NSV_AEREND
      !
      CALL ch_aer_dep(psv(:,ibeg:iend),psfts(:,ibeg:iend),zustar,zresa_water,pta,prhoa)     
      END IF

  ELSE
    ibeg = chf%SVF%NSV_AERBEG
    iend = chf%SVF%NSV_AEREND
    !
    psfts(:,ibeg:iend) =0.
    IF(ibeg.LT.iend) psfts(:,ibeg:iend) =0.
  ENDIF
  !
ENDIF

IF (chf%SVF%NDSTEQ>0) THEN
  !
  ibeg = chf%SVF%NSV_DSTBEG
  iend = chf%SVF%NSV_DSTEND
  !
  CALL dslt_dep(psv(:,ibeg:iend), psfts(:,ibeg:iend), zustar, zresa_water, pta, prhoa, &
                dst%XEMISSIG_DST, dst%XEMISRADIUS_DST, jpmode_dst, xdensity_dst, &
                xmolarweight_dst, zconvertfacm0_dst, zconvertfacm6_dst, zconvertfacm3_dst, &
                lvarsig_dst, lrgfix_dst, cvermod  )  

  CALL massflux2momentflux(         &
    psfts(:,ibeg:iend),             & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
    prhoa,                          & !I [kg/m3] air density
    dst%XEMISRADIUS_DST,                &!I [um] emitted radius for the modes (max 3)
    dst%XEMISSIG_DST,                   &!I [-] emitted sigma for the different modes (max 3)
    ndstmde,                        &
    zconvertfacm0_dst,              &
    zconvertfacm6_dst,              &
    zconvertfacm3_dst,              &
    lvarsig_dst, lrgfix_dst         )  
ENDIF


IF (chf%SVF%NSLTEQ>0) THEN
  !
  ibeg = chf%SVF%NSV_SLTBEG
  iend = chf%SVF%NSV_SLTEND
  !
  CALL dslt_dep(psv(:,ibeg:iend), psfts(:,ibeg:iend), zustar, zresa_water, pta, prhoa, &
                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(         &
    psfts(:,ibeg:iend),             & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
    prhoa,                          & !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

!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Inline diagnostics
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (f%CFLK_FLUX=='FLAKE') THEN  !compute some variables not present in FLake code
!
  zch = zeps
!
  WHERE (abs((f%XTS(:) - pta(:) * zexns(:)/zexna(:))) > 1.e-2 .AND. zwind(:)/=0.)
     zch = max(zeps,psfth / (xcpd * prhoa(:) * zwind(:) * (f%XTS(:) - pta(:) * zexns(:)/zexna(:))) * zexns(:))
  END WHERE
!
  zcdn(:) = (xkarman/log(puref(:)/f%XZ0(:)))**2
  zcd(:) = max(zeps,zcd(:))
!
ENDIF
!
 CALL diag_inline_flake_n(dgo, d, dc, f, &
                          ptstep, pta,  zqa, ppa, pps, prhoa, pu,    &
                          pv, pzref, puref, prain, psnow,            &
                          zcd, zcdn, zch, zri, zhu,                  &
                          zz0h, zqsat, psfth, psftq, psfu, psfv,     &
                          pdir_sw, psca_sw, plw, zdir_alb, zsca_alb, &
                          zle, zlei, zsubl, zlwup, zalb, zswe        )  
!
!-------------------------------------------------------------------------------------
!
 CALL diag_misc_flake_n(dmf, f)
!
!-------------------------------------------------------------------------------
!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
!-------------------------------------------------------------------------------
!
ptsurf(:) = f%XTS  (:)
pz0(:) = f%XZ0  (:)
pz0h(:) = zz0h(:)
pqsurf(:) = zqsat(:)
!
!-------------------------------------------------------------------------------------
!Radiative properties at time t+1 (see by the atmosphere) in order to close
!the energy budget between surfex and the atmosphere
!-------------------------------------------------------------------------------------
!
 CALL update_rad_flake(f,pzenith2,pdir_alb,psca_alb,pemis,ptsrad  )                       
!                         
IF (lhook) CALL dr_hook('COUPLING_FLAKE_N',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
END SUBROUTINE coupling_flake_n