coupling_ideal_flux.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_ideal_flux (DGL, &
                                hprogram, hcoupling, ptimec,                                  &
                 ptstep, kyear, kmonth, kday, ptime, ki, ksv, ksw, ptsun, pzenith, pazim,    &
                 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_IDEAL_FLUX * - Computes the surface fluxes for the temperature, 
!!    vapor, horizontal components of the wind and the scalar variables.   
!!
!!    PURPOSE
!!    -------
!       Give prescribed values of the surface fluxes for the potential 
!     temperature, the vapor, the horizontal components of the wind and the 
!     scalar variables. These fluxes are unsteady when a diurnal cycle 
!     is taken into account.
!
!!**  METHOD
!!    ------
!!         A temporal interpolation is performed to recover the values of the 
!!    fluxes at every instant of the simulation. The different values of the
!!    prescribed fluxes are given at their declarations.
!!         For the wind, z0 can also be prescribed and the flux is determined 
!!    with a neutral drag coefficient.
!!
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!     (from J. Cuxart and J. Stein)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      Modified    09/2012  : J. Escobar , SIZE(PTA) not allowed without-interface , replace by KI
!!      B. Decharme  04/2013 new coupling variables
!!      P. Le Moigne 03/2015 add diagnostics IDEAL case
!!      R. ROEHRIG   06/2015 add PTSTEP in TEMP_FORC_DISTS
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
!
USE modd_diag_ideal_n, ONLY : diag_ideal_t
!
USE modd_csts,       ONLY : xrd, xcpd, xp00, xpi, xlvtt, xday, xkarman, xtt, &
                            xlstt, xstefan
USE modd_ideal_flux, ONLY : nforcf, nforct, xsfth, xsftq, xsfts, xsfco2, &
                            custartype, xustar, xz0, xalb, xemis, xtsrad, &
                            xtimef, xtimet 
USE modd_surf_par,   ONLY : xundef
!
USE mode_sbls
USE mode_thermos
!
USE modi_diag_inline_ideal_n
USE modi_surface_aero_cond
USE modi_surface_ri
USE modi_surface_cd
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*       0.1   declarations of arguments
! 
!
TYPE(diag_ideal_t), INTENT(INOUT) :: dgl
!
 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       (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)  :: pzs       ! atmospheric model orography           (m)
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
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
!
REAL, DIMENSION(KI)  :: zz0     ! roughness length
REAL, DIMENSION(KI)  :: zlmo    ! Monin-Obuhkov length
REAL, DIMENSION(KI)  :: ztha    ! air potential temperature
REAL, DIMENSION(KI)  :: zrva    ! water vapor mixing ratio
REAL, DIMENSION(KI)  :: zustar  ! friction velocity
REAL, DIMENSION(KI)  :: zwind   ! wind
REAL, DIMENSION(KI)  :: zq0     ! surface temperature flux (mK/s)
REAL, DIMENSION(KI)  :: ze0     ! surface vapor flux (mkg/kg/s)
REAL, DIMENSION(KI)  :: zqa     ! specific humidity (kg/kg)
REAL, DIMENSION(KI)  :: zdircoszw
REAL, DIMENSION(KI)  :: zexns, zexna
REAL, DIMENSION(KI)  :: zac, zra, zch, zcd, zcdn, zri
REAL, DIMENSION(KI)  :: zhu
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,KSW) :: zdir_alb   ! Direct albedo at time t , 
REAL, DIMENSION(KI,KSW) :: zsca_alb   ! Diffuse albedo at time t
!
INTEGER              :: iswb   ! number of shortwave spectral bands
INTEGER              :: jswb   ! loop counter on shortwave spectral bands
!
REAL                        :: zalpha  ! interpolation coefficient
INTEGER                     :: ihourf  ! number of hours since midnight
INTEGER                     :: ihourt
INTEGER                     :: jiter   ! convergence loop counter
INTEGER                     :: jsv     ! loop on scalar variables
!
LOGICAL                     :: gcall_lmo ! flag in non-neutral case
!
INTEGER                     :: iluout  ! output listing logical unit
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('COUPLING_IDEAL_FLUX',0,zhook_handle)
IF (htest/='OK') THEN
  CALL abor1_sfx('COUPLING_IDEAL_FLUX: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
!----------------------------------------------------------------------------------
zle(:)   = xundef
zlei(:)   = xundef
zsubl(:)   = xundef
zlwup(:)   = xundef
!
!*       2.    COMPUTE TIME
!              ------------
!
 CALL temp_forc_dists(ptimec,ptstep,nforcf,xtimef,ihourf,zalpha)
!
!----------------------------------------------------------------------------------
!
!*       3.    CONS. TEMPERATURE SURFACE FLUX
!              ------------------------------
!
psfth(:) = xsfth(ihourf) + ( xsfth(ihourf+1)-xsfth(ihourf) )*zalpha  
!
gcall_lmo = ( xsfth(ihourf) + ( xsfth(ihourf+1)-xsfth(ihourf) )*zalpha ) /=0.
!----------------------------------------------------------------------------------
!
!*       4.    CONS. MIXING RATIO SURFACE FLUX
!              -------------------------------
!
psftq(:) = xsftq(ihourf) + ( xsftq(ihourf+1)-xsftq(ihourf) )*zalpha
!
gcall_lmo = gcall_lmo .OR. ( xsftq(ihourf) + ( xsftq(ihourf+1)-xsftq(ihourf) )*zalpha ) /=0.
!----------------------------------------------------------------------------------
!
!*       5.    WIND SURFACE FLUX
!              -----------------
!
!*       5.1   wind
!
zwind(:) = sqrt(pu**2+pv**2)
!
!*       5.2   Compute the surface stresses
!
SELECT CASE (custartype)
!
!
  CASE ('USTAR')
    !  when u* is prescribed
    zustar(:) = xustar(ihourf) + ( xustar(ihourf+1)-xustar(ihourf) )*zalpha
    ! 
  CASE ('Z0   ')
    !
    !* spatialized roughness length
    !
    zz0(:) = xz0
    !
    !* water mixing ratio
    !
    zrva(:) = 0.
    zqa(:)  = pqa(:) / prhoa(:)
    !
    WHERE (zqa(:)/=0.) zrva(:) = 1./(1./zqa(:) - 1.)
    !
    !* air potential temperature
    ztha(:) = pta(:) * (xp00/ppa(:))**(xrd/xcpd)
    !
    !* cinematic surface fluxes
    zq0(:) = psfth(:) / xcpd / prhoa(:)
    ze0(:) = psftq(:)        / prhoa(:)
    !
    !
    !* neutral case, as guess
    zlmo(:) = xundef
    zustar(:) = ustar(zwind(:),pzref(:),zz0(:),zlmo(:))
    !
    !* iterations in non-neutral case
    IF (gcall_lmo) THEN
      zustar(:) = max( zustar(:), 0.01 )
      DO jiter=1,10
        zlmo(:) = lmo(zustar(:),ztha(:),zrva(:),zq0(:),ze0(:))
        zustar(:) = ustar(zwind(:),pzref(:),zz0(:),zlmo(:))
      END DO
    END IF
    !
    !
END SELECT
!
psfu = 0.
psfv = 0.
WHERE (zwind>0.)
  psfu = - prhoa * zustar**2 * pu / zwind
  psfv = - prhoa * zustar**2 * pv / zwind
END WHERE
!
!-------------------------------------------------------------------------------
!  
!*       6.    SCALAR VARIABLES FLUXES.
!              -----------------------
!
DO jsv=1,SIZE(psfts,2)
  psfts(:,jsv) = xsfts(ihourf,jsv) + ( xsfts(ihourf+1,jsv)-xsfts(ihourf,jsv) )*zalpha
END DO
!
!-------------------------------------------------------------------------------
!  
!*       7.    CO2 FLUXES.
!              ----------
!
psfco2(:) = xsfco2(ihourf) + ( xsfco2(ihourf+1)-xsfco2(ihourf) )*zalpha
!
!-------------------------------------------------------------------------------
!  
!*       8.    OTHER OUTPUTS (RADIATIVE QUANTITIES) SET TO A CONSTANT VALUE
!              ------------------------------------
!
 CALL temp_forc_dists(ptimec,ptstep,nforct,xtimet,ihourt,zalpha)
!
ptrad(:) = xtsrad(ihourt) + ( xtsrad(ihourt+1)-xtsrad(ihourt) )*zalpha
!
pdir_alb = xalb
psca_alb = xalb
pemis    = xemis
!  
ptsurf(:) = ptrad(:)
pz0(:) = xz0
pz0h(:) = xz0
pqsurf(:) = qsat(ptsurf(:),pps(:))
!
!-------------------------------------------------------------------------------
!
!*       9.    INLINE DIAGNOSTICS
!              ------------------
!
zdircoszw = 1.
!
zexns(:)     = (pps(:)/xp00)**(xrd/xcpd)
zexna(:)     = (ppa(:)/xp00)**(xrd/xcpd)
zqa(:)       = pqa(:) / prhoa(:)
!
zhu(:)=1.
!
WHERE (ptsurf(:)<xtt)
  zle(:) = psftq(:) * xlstt
  zlei(:) = psftq(:) * xlstt
  zsubl(:) = psftq(:)
ELSEWHERE
  zle(:) = psftq(:) * xlvtt
  zlei(:) = 0.0
  zsubl(:) = 0.0
END WHERE
!
zlwup(:)=(1.-pemis(:))*plw(:)+pemis(:)*xstefan*ptsurf(:)**4
!
 CALL surface_ri(ptsurf,pqsurf,zexns,zexna,pta,pqa,  &
                pzref, puref, zdircoszw,zwind,zri)

 CALL surface_aero_cond(zri, pzref, puref, zwind, pz0, pz0h, zac, zra, zch)

 CALL surface_cd(zri, pzref, puref, pz0, pz0h, zcd, zcdn)

 CALL diag_inline_ideal_n(dgl, ptstep, pta, ptsurf, zqa, ppa, pps, prhoa, pu,      &
                            pv, pzref, puref, prain, psnow,                  &
                            zcd, zcdn, zch, zri, zhu, pz0,                   &
                            pz0h, pqsurf, psfth, psftq, psfu, psfv,          &
                            pdir_sw, psca_sw, plw, pdir_alb, psca_alb,       &
                            zle, zlei, zsubl, zlwup                          )
!
IF (lhook) CALL dr_hook('COUPLING_IDEAL_FLUX',1,zhook_handle)
!
!-------------------------------------------------------------------------------
 CONTAINS
!
SUBROUTINE temp_forc_dists (PTIMEIN,PSTEP,KFORC,PTIMES,KHOUR,PALPHA)
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: ptimein
REAL, INTENT(IN) :: pstep
INTEGER, INTENT(IN) :: kforc
REAL, DIMENSION(:), INTENT(IN) :: ptimes
INTEGER, INTENT(OUT):: khour
REAL, INTENT(OUT):: palpha
!
INTEGER :: jt
REAL :: ztimein
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('COUPLING_IDEAL_FLUX:TEMP_FORC_DISTS',0,zhook_handle)
!
ztimein = ptimein
!
IF (ptimes(kforc)==xundef) THEN
  khour = 1
  palpha = 0.
ELSEIF (ztimein<ptimes(1).OR.ztimein>ptimes(kforc)) THEN
    WRITE(*,*) 'COUPLING_IDEAL_FLUX', ztimein, ptimes(1), ptimes(kforc)      
  CALL abor1_sfx("COUPLING_IDEAL_FLUX:TEMP_FORC_DISTS: PTIMEC OUT OF BOUNDS!!!")
ELSEIF (ztimein==ptimes(kforc)) THEN
  khour = kforc
  palpha = 0.  
ELSE
 DO jt = kforc,1,-1
   IF (ptimein.GE.ptimes(jt)) THEN
     khour = jt
     EXIT
   ENDIF
  ENDDO      
  palpha = (ptimein-ptimes(khour)) / (ptimes(khour+1)-ptimes(khour))
ENDIF
!
IF (lhook) CALL dr_hook('COUPLING_IDEAL_FLUX:TEMP_FORC_DISTS',1,zhook_handle)
!
END SUBROUTINE temp_forc_dists
!
END SUBROUTINE coupling_ideal_flux