coupling_watfluxn.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_watflux_n (WM, DST, SLT, &
                               hprogram, hcoupling, ptimec,                                   &
                 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,                                    &
                 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_WATFLUX_n * - Driver for WATER_FLUX scheme for lakes
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      B. Decharme 03/2009 TS_WATER could change during a run => ALB and EMIS
!!      V. Masson   05/2009 Implicitation of momentum fluxes
!!      B. Decharme 09/2009 Radiative properties at time t+1 in order to close
!                           the energy budget between surfex and the atmosphere 
!!      B. Decharme 01/2010 Add XTT
!!      B. Decharme 09/2012 New wind implicitation
!!      P. LeMoigne 04/2013 Wind implicitation displaced
!!      B. Decharme  04/2013 new coupling variables
!!----------------------------------------------------------------------------
!
!
!
USE modd_surfex_n, ONLY : watflux_model_t
USE modd_dst_n, ONLY : dst_t
USE modd_slt_n, ONLY : slt_t
USE modd_watflux_n, ONLY : watflux_t
!
USE modd_reprod_oper, ONLY : cimplicit_wind
!
USE modd_csts,       ONLY : xrd, xcpd, xp00, xtt, xday, xtts
USE modd_surf_par,   ONLY : xundef
USE modd_sfx_oasis,  ONLY : lcpl_seaice
USE modd_water_par
!
!
USE modi_water_flux
USE modi_add_forecast_to_date_surf
USE modi_diag_inline_watflux_n
USE modi_ch_aer_dep
USE modi_ch_dep_water
USE modi_dslt_dep
USE modi_update_rad_water
USE modi_interpol_ts_water_mth
!
USE mode_dslt_surf
USE modd_dst_surf
USE modd_slt_surf
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
USE modi_coupling_iceflux_n
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(watflux_model_t), INTENT(INOUT) :: wm
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
REAL,                INTENT(IN)  :: ptimec    ! cumulated time since beginning of simulation
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) :: 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   (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)  :: zexna      ! Exner function at forcing level
REAL, DIMENSION(KI)  :: zexns      ! Exner function at surface level
REAL, DIMENSION(KI)  :: zwind      ! Wind
REAL, DIMENSION(KI)  :: zcd        ! Drag coefficient
REAL, DIMENSION(KI)  :: zcdn       ! Neutral Drag coefficient
REAL, DIMENSION(KI)  :: zch        ! Heat transfer coefficient
REAL, DIMENSION(KI)  :: zri        ! Richardson number
REAL, DIMENSION(KI)  :: zhu        ! Near surface relative humidity
REAL, DIMENSION(KI)  :: zresa_water! aerodynamical resistance
REAL, DIMENSION(KI)  :: zustar     ! friction velocity (m/s)
REAL, DIMENSION(KI)  :: zustar2    ! square of friction velocity (m2/s2)
REAL, DIMENSION(KI)  :: zz0h       ! heat roughness length over water
REAL, DIMENSION(KI)  :: zqsat      ! humidity at saturation
REAL, DIMENSION(KI)  :: zqa        ! specific humidity (kg/kg)
REAL, DIMENSION(KI)  :: zemis      ! Emissivity at time t
REAL, DIMENSION(KI)  :: ztrad      ! Radiative temperature at time t
REAL, DIMENSION(KI)  :: zsfth_ice  ! Sea ice flux of heat
REAL, DIMENSION(KI)  :: zsftq_ice  ! Sea ice flux of ice sublimation
REAL, DIMENSION(KI)  :: zwork      ! Work array
!
REAL, DIMENSION(KI,KSW) :: zdir_alb   ! Direct albedo at time t
REAL, DIMENSION(KI,KSW) :: zsca_alb   ! Diffuse albedo at time t
!
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      
!
LOGICAL                    :: ghandle_sic = .false. ! no sea-ice model
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------------
! Preliminaries:
!-------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('COUPLING_WATFLUX_N',0,zhook_handle)
IF (htest/='OK') THEN
  CALL abor1_sfx('COUPLING_WATFLUXN: FATAL ERROR DURING ARGUMENT TRANSFER')        
END IF
!-------------------------------------------------------------------------------------
!
zexna(:) = xundef
zexns(:) = xundef
zwind(:) = xundef
zcd(:) = xundef    
zcdn(:) = xundef
zch(:) = xundef
zri(:) = xundef
zhu(:) = xundef
zresa_water(:) = xundef
zustar(:) = xundef
zustar2(:) = xundef
zz0h(:) = xundef
zqsat(:) = xundef
zemis(:) = xundef
ztrad(:) = xundef
zwork(:) = xundef
zdir_alb(:,:) = xundef
zsca_alb(:,:) = xundef
!
IF(lcpl_seaice)THEN
  zsftq_ice(:) = xundef
  zsfth_ice(:) = xundef
ENDIF
!
!-------------------------------------------------------------------------------------
!
zexns(:)     = (pps(:)/xp00)**(xrd/xcpd)
zexna(:)     = (ppa(:)/xp00)**(xrd/xcpd)
!
zwind(:) = sqrt(pu(:)**2+pv(:)**2)
!
psfts(:,:) = 0.
!
zhu = 1.
!
zqa(:) = pqa(:) / prhoa(:)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Time evolution
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
wm%W%TTIME%TIME = wm%W%TTIME%TIME + ptstep
 CALL add_forecast_to_date_surf(wm%W%TTIME%TDATE%YEAR,wm%W%TTIME%TDATE%MONTH,wm%W%TTIME%TDATE%DAY,wm%W%TTIME%TIME)
!
!--------------------------------------------------------------------------------------
! Fluxes over water according to Charnock formulae
!--------------------------------------------------------------------------------------
!
!
 CALL water_flux(wm%W%XZ0, pta, zexna, prhoa, wm%W%XTS, zexns, zqa, prain,    &
                psnow, xtt, zwind, pzref, puref, pps, ghandle_sic, &
                zqsat, psfth, psftq, zustar, zcd, zcdn, zch, zri,  &
                zresa_water, zz0h                                  )  
!
!-------------------------------------------------------------------------------------
!radiative properties at t
!-------------------------------------------------------------------------------------
!
iswb = SIZE(psw_bands)
!
DO jswb=1,iswb
     zdir_alb(:,jswb) = wm%W%XDIR_ALB(:)
     zsca_alb(:,jswb) = wm%W%XSCA_ALB(:)
END DO
!
zemis  = wm%W%XEMIS
ztrad  = wm%W%XTS
!
!-------------------------------------------------------------------------------------
!Specific fields for GELATO when using earth system model 
!(intermediate step before explicit sea and ice fluxes comutation)
!-------------------------------------------------------------------------------------
!
IF(lcpl_seaice)THEN   
  CALL coupling_iceflux_n(ki, pta, zexna, prhoa, wm%W%XTICE, zexns,      &
                            zqa, prain, psnow, zwind, pzref, puref, &
                            pps, wm%W%XTS, xtt, zsfth_ice, zsftq_ice     )  
ENDIF
!
!-------------------------------------------------------------------------------------
! Scalar fluxes:
!-------------------------------------------------------------------------------------
!
!
IF (wm%CHW%SVW%NBEQ>0) THEN
  IF (wm%CHW%CCH_DRY_DEP == "WES89") THEN
    CALL ch_dep_water(zresa_water, zustar, pta, ztrad,      &
                          psv(:,wm%CHW%SVW%NSV_CHSBEG:wm%CHW%SVW%NSV_CHSEND),       &
                          wm%CHW%SVW%CSV(wm%CHW%SVW%NSV_CHSBEG:wm%CHW%SVW%NSV_CHSEND),         &
                          wm%CHW%XDEP(:,1:wm%CHW%SVW%NBEQ) )  

   psfts(:,wm%CHW%SVW%NSV_CHSBEG:wm%CHW%SVW%NSV_CHSEND) = - psv(:,wm%CHW%SVW%NSV_CHSBEG:wm%CHW%SVW%NSV_CHSEND)  &
                                               * wm%CHW%XDEP(:,1:wm%CHW%SVW%NBEQ)  
     IF (wm%CHW%SVW%NAEREQ > 0 ) THEN
        CALL ch_aer_dep(psv(:,wm%CHW%SVW%NSV_AERBEG:wm%CHW%SVW%NSV_AEREND),&
                          psfts(:,wm%CHW%SVW%NSV_AERBEG:wm%CHW%SVW%NSV_AEREND),&
                          zustar,zresa_water,pta,prhoa)     
      END IF

  ELSE
    psfts(:,wm%CHW%SVW%NSV_CHSBEG:wm%CHW%SVW%NSV_CHSEND) =0.
    IF(wm%CHW%SVW%NSV_AERBEG.LT.wm%CHW%SVW%NSV_AEREND) psfts(:,wm%CHW%SVW%NSV_AERBEG:wm%CHW%SVW%NSV_AEREND) =0.
  ENDIF
ENDIF
!
IF (wm%CHW%SVW%NDSTEQ>0) THEN
  CALL dslt_dep(psv(:,wm%CHW%SVW%NSV_DSTBEG:wm%CHW%SVW%NSV_DSTEND), psfts(:,wm%CHW%SVW%NSV_DSTBEG:wm%CHW%SVW%NSV_DSTEND),   &
                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(:,wm%CHW%SVW%NSV_DSTBEG:wm%CHW%SVW%NSV_DSTEND), & !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 (wm%CHW%SVW%NSLTEQ>0) THEN
  CALL dslt_dep(psv(:,wm%CHW%SVW%NSV_SLTBEG:wm%CHW%SVW%NSV_SLTEND), psfts(:,wm%CHW%SVW%NSV_SLTBEG:wm%CHW%SVW%NSV_SLTEND),   &
                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(:,wm%CHW%SVW%NSV_SLTBEG:wm%CHW%SVW%NSV_SLTEND), & !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
!
!-------------------------------------------------------------------------------------
! Outputs fluxes at time t+1
!-------------------------------------------------------------------------------------
!
! Momentum fluxes
!
IF(cimplicit_wind=='OLD')THEN
! old implicitation
  zustar2(:) = (zcd(:)*zwind(:)*ppew_b_coef(:)) /          &
               (1.0-prhoa(:)*zcd(:)*zwind(:)*ppew_a_coef(:))
ELSE
! new implicitation
  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
!
psfu = 0.
psfv = 0.
WHERE (zwind(:)>0.)
  psfu(:) = - prhoa(:) * zustar2(:) * pu(:) / zwind(:)
  psfv(:) = - prhoa(:) * zustar2(:) * pv(:) / zwind(:)
END WHERE
!
! CO2 flux
!
psfco2(:)       =  0.0    ! Assumes no CO2 emission over water bodies
!
!-------------------------------------------------------------------------------------
! Inline diagnostics at time t for TS and TRAD
!-------------------------------------------------------------------------------------
!
 CALL diag_inline_watflux_n(wm%DGW, wm%W, &
                            ptstep,pta, zqa, ppa, pps, prhoa, pu, pv, pzref,  &
                             puref, zcd, zcdn, zch, zri, zhu, zz0h, zqsat,     &
                             psfth, psftq, psfu, psfv, pdir_sw, psca_sw, plw,       &
                             zdir_alb, zsca_alb, zemis, ztrad, prain, psnow,        &
                             zsfth_ice, zsftq_ice                            )
!
!-------------------------------------------------------------------------------
! IMPOSED MONTHLY TS AT TIME t+1
!-------------------------------------------------------------------------------
!  
IF (wm%W%LINTERPOL_TS.AND.mod(wm%W%TTIME%TIME,xday) == 0.) THEN
   CALL interpol_ts_water_mth(wm%W, &
                              wm%W%TTIME%TDATE%YEAR,wm%W%TTIME%TDATE%MONTH,wm%W%TTIME%TDATE%DAY,wm%W%XTS)
ENDIF
!
!-------------------------------------------------------------------------------
!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(:) = wm%W%XTS  (:)
pz0(:) = wm%W%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_water(wm%W%CWAT_ALB,wm%W%XTS,pzenith2,xtt,wm%W%XEMIS,wm%W%XDIR_ALB, &
                       wm%W%XSCA_ALB,pdir_alb,psca_alb,pemis,ptrad    )  
!
!=======================================================================================
!
IF (lhook) CALL dr_hook('COUPLING_WATFLUX_N',1,zhook_handle)
!
END SUBROUTINE coupling_watflux_n