Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	phylmd/calcul_fluxs_mod.F90	Lines:	68	70	97.1 %
Date:	2023-06-30 12:51:15	Branches:	56	86	65.1 %


!
! $Id: calcul_fluxs_mod.F90 3815 2021-02-01 14:30:57Z lguez $
!
MODULE calcul_fluxs_mod

  IMPLICIT NONE

CONTAINS
  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
       tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
       precip_rain, precip_snow, snow, qsurf, &
       radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
       fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
       tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)


    USE dimphy, ONLY : klon
    USE indice_sol_mod
    use sens_heat_rain_m, only: sens_heat_rain

    INCLUDE "clesphys.h"

! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
! une temperature de surface (au cas ou ok_veget = false)
!
! L. Fairhead 4/2000
!
! input:
!   knon         nombre de points a traiter
!   nisurf       surface a traiter
!   tsurf        temperature de surface
!   p1lay        pression 1er niveau (milieu de couche)
!   cal          capacite calorifique du sol
!   beta         evap reelle
!   cdragh       coefficient d'echange temperature
!   cdragq       coefficient d'echange evaporation
!   ps           pression au sol
!   precip_rain  precipitations liquides
!   precip_snow  precipitations solides
!   snow         champs hauteur de neige
!   runoff       runoff en cas de trop plein
!   petAcoef     coeff. A de la resolution de la CL pour t
!   peqAcoef     coeff. A de la resolution de la CL pour q
!   petBcoef     coeff. B de la resolution de la CL pour t
!   peqBcoef     coeff. B de la resolution de la CL pour q
!   radsol       rayonnement net aus sol (LW + SW)
!   dif_grnd     coeff. diffusion vers le sol profond
!
! output:
!   tsurf_new    temperature au sol
!   qsurf        humidite de l'air au dessus du sol
!   fluxsens     flux de chaleur sensible
!   fluxlat      flux de chaleur latente
!   dflux_s      derivee du flux de chaleur sensible / Ts
!   dflux_l      derivee du flux de chaleur latente  / Ts
!   sens_prec_liq flux sensible li� aux echanges de precipitations liquides
!   sens_prec_sol                                    precipitations solides
!   lat_prec_liq  flux latent li� aux echanges de precipitations liquides
!   lat_prec_sol                                  precipitations solides

    INCLUDE "YOETHF.h"
    INCLUDE "FCTTRE.h"
    INCLUDE "YOMCST.h"

! Parametres d'entree
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: knon, nisurf
    REAL   , INTENT(IN)                  :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
    REAL, DIMENSION(klon), INTENT(IN)    :: ps, q1lay
    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf, p1lay, cal, beta, cdragh,cdragq
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)    :: radsol, dif_grnd
    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay, u1lay, v1lay,gustiness
    REAL,                  INTENT(IN)    :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else)

    real, intent(in), optional:: rhoa(:) ! (knon)
    ! density of moist air  (kg / m3)

! Parametres entree-sorties
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT) :: snow  ! snow pas utile

! Parametres sorties
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new, evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l
    REAL, intent(out), OPTIONAL:: sens_prec_liq(:), sens_prec_sol(:) ! (knon)
    REAL, DIMENSION(klon), OPTIONAL      :: lat_prec_liq, lat_prec_sol

! Variables locales
!****************************************************************************************
    INTEGER                              :: i
    REAL, DIMENSION(klon)                :: zx_mh, zx_nh, zx_oh
    REAL, DIMENSION(klon)                :: zx_mq, zx_nq, zx_oq
    REAL, DIMENSION(klon)                :: zx_pkh, zx_dq_s_dt, zx_qsat
    REAL, DIMENSION(klon)                :: zx_sl, zx_coefh, zx_coefq, zx_wind
    REAL, DIMENSION(klon)                :: d_ts
    REAL                                 :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
    REAL                                 :: qsat_new, q1_new
    REAL, PARAMETER                      :: t_grnd = 271.35, t_coup = 273.15
    REAL, PARAMETER                      :: max_eau_sol = 150.0
    CHARACTER (len = 20)                 :: modname = 'calcul_fluxs'
    LOGICAL                              :: fonte_neige
    LOGICAL, SAVE                        :: check = .FALSE.
    !$OMP THREADPRIVATE(check)

! End definition
!****************************************************************************************

    IF (check) WRITE(*,*)'Entree ', modname,' surface = ',nisurf

    IF (check) THEN
       WRITE(*,*)' radsol (min, max)', &
            MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
    ENDIF

! Traitement neige et humidite du sol
!****************************************************************************************
!
!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
!!PB test
!!$    if (nisurf == is_oce) then
!!$      snow = 0.
!!$      qsol = max_eau_sol
!!$    else
!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
!!$    endif
!!$    IF (nisurf /= is_ter) qsol = max_eau_sol


!
! Initialisation
!****************************************************************************************
    evap = 0.
    fluxsens=0.
    fluxlat=0.
    dflux_s = 0.
    dflux_l = 0.
    if (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
    if (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
!
! zx_qs = qsat en kg/kg
!****************************************************************************************
    DO i = 1, knon
       zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
       IF (thermcep) THEN
          zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
          zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
          zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
          zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
          zx_qs=MIN(0.5,zx_qs)
          zcor=1./(1.-retv*zx_qs)
          zx_qs=zx_qs*zcor
          zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
               /RLVTT / zx_pkh(i)
       ELSE
          IF (tsurf(i).LT.t_coup) THEN
             zx_qs = qsats(tsurf(i)) / ps(i)
             zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
                  / zx_pkh(i)
          ELSE
             zx_qs = qsatl(tsurf(i)) / ps(i)
             zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
                  / zx_pkh(i)
          ENDIF
       ENDIF
       zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
       zx_qsat(i) = zx_qs
       zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)**2+v1lay(i)**2)
       zx_coefh(i) = cdragh(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
       zx_coefq(i) = cdragq(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
!      zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)**2+v1lay(i)**2) &
!                * p1lay(i)/(RD*t1lay(i))
!      zx_coefh(i) = cdragh(i) * zx_wind(i)
!      zx_coefq(i) = cdragq(i) * zx_wind(i)
    ENDDO


! === Calcul de la temperature de surface ===
! zx_sl = chaleur latente d'evaporation ou de sublimation
!****************************************************************************************

    DO i = 1, knon
       zx_sl(i) = RLVTT
       IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
    ENDDO


    DO i = 1, knon
! Q
       zx_oq(i) = 1. - (beta(i) * zx_coefq(i) * peqBcoef(i) * dtime)
       zx_mq(i) = beta(i) * zx_coefq(i) * &
            (peqAcoef(i) -             &
! conv num avec precedente version
            fqsat * zx_qsat(i) + fqsat * zx_dq_s_dt(i) * tsurf(i))  &
!           fqsat * ( zx_qsat(i) - zx_dq_s_dt(i) * tsurf(i)) ) &
            / zx_oq(i)
       zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
            / zx_oq(i)

! H
       zx_oh(i) = 1. - (zx_coefh(i) * petBcoef(i) * dtime)
       zx_mh(i) = zx_coefh(i) * petAcoef(i) / zx_oh(i)
       zx_nh(i) = - (zx_coefh(i) * RCPD * zx_pkh(i))/ zx_oh(i)

! Tsurface
       tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
            (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
            + dif_grnd(i) * t_grnd * dtime)/ &
            ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
            zx_nh(i) + zx_sl(i) * zx_nq(i)) &
            + dtime * dif_grnd(i))

!
! Y'a-t-il fonte de neige?
!
!    fonte_neige = (nisurf /= is_oce) .AND. &
!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
!     & .AND. (tsurf_new(i) >= RTT)
!    if (fonte_neige) tsurf_new(i) = RTT
       d_ts(i) = tsurf_new(i) - tsurf(i)
!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)

!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
       evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
       fluxlat(i) = - evap(i) * zx_sl(i)
       fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)

! Derives des flux dF/dTs (W m-2 K-1):
       dflux_s(i) = zx_nh(i)
       dflux_l(i) = (zx_sl(i) * zx_nq(i))

! Nouvelle valeure de l'humidite au dessus du sol
       qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
       q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
       qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
!
! en cas de fonte de neige
!
!    if (fonte_neige) then
!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
!      bilan_f = max(0., bilan_f)
!      fq_fonte = bilan_f / zx_sl(i)
!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
!      qsol(i) = qsol(i) + (fq_fonte * dtime)
!    endif
!!$    if (nisurf == is_ter)  &
!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
!!$    qsol(i) = min(qsol(i), max_eau_sol)
!
! calcul de l'enthalpie des precipitations liquides et solides
       if (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
            = - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
            q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
       if (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0.
       ! On calcule par rapport a T=0
       !! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i)
       !! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i)

       if (PRESENT(lat_prec_liq))                    &
         lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT)
       if (PRESENT(lat_prec_sol))                    &
         lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
    ENDDO


!**************************************************************************
!
  END SUBROUTINE calcul_fluxs
!
!****************************************************************************************
!
  SUBROUTINE calcul_flux_wind(knon, dtime, &
       u0, v0, u1, v1, gustiness, cdrag_m, &
       AcoefU, AcoefV, BcoefU, BcoefV, &
       p1lay, t1lay, &
       flux_u1, flux_v1)

    USE dimphy
    INCLUDE "YOMCST.h"
    INCLUDE "clesphys.h"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: knon
    REAL, INTENT(IN)                     :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: u0, v0  ! u and v at niveau 0
    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness  ! u and v at niveau 1
    REAL, DIMENSION(klon), INTENT(IN)    :: cdrag_m ! cdrag pour momentum
    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay   ! pression 1er niveau (milieu de couche)
    REAL, DIMENSION(klon), INTENT(IN)    :: t1lay   ! temperature
! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1
    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v1

! Local variables
!****************************************************************************************
    INTEGER                              :: i
    REAL                                 :: mod_wind, buf

!****************************************************************************************
! Calculate the surface flux
!
!****************************************************************************************
    DO i=1,knon
       mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
       buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
       flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
       flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
    END DO

  END SUBROUTINE calcul_flux_wind
!
!****************************************************************************************
!
END MODULE calcul_fluxs_mod


Generated by: GCOVR (Version 4.2)

Line	Branch	Exec	Source
1			!
2			! $Id: calcul_fluxs_mod.F90 3815 2021-02-01 14:30:57Z lguez $
3			!
4			MODULE calcul_fluxs_mod
5
6			IMPLICIT NONE
7
8			CONTAINS
9		473954	SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
10			tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
11			precip_rain, precip_snow, snow, qsurf, &
12		1152	radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, &
13			fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
14			tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
15	✓✓✓✗ ✗✓✓✓ ✓✗✗✓ ✓✓✓✗ ✗✓	1152	sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
16
17
18			USE dimphy, ONLY : klon
19			USE indice_sol_mod
20			use sens_heat_rain_m, only: sens_heat_rain
21
22			INCLUDE "clesphys.h"
23
24			! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
25			! une temperature de surface (au cas ou ok_veget = false)
26			!
27			! L. Fairhead 4/2000
28			!
29			! input:
30			! knon nombre de points a traiter
31			! nisurf surface a traiter
32			! tsurf temperature de surface
33			! p1lay pression 1er niveau (milieu de couche)
34			! cal capacite calorifique du sol
35			! beta evap reelle
36			! cdragh coefficient d'echange temperature
37			! cdragq coefficient d'echange evaporation
38			! ps pression au sol
39			! precip_rain precipitations liquides
40			! precip_snow precipitations solides
41			! snow champs hauteur de neige
42			! runoff runoff en cas de trop plein
43			! petAcoef coeff. A de la resolution de la CL pour t
44			! peqAcoef coeff. A de la resolution de la CL pour q
45			! petBcoef coeff. B de la resolution de la CL pour t
46			! peqBcoef coeff. B de la resolution de la CL pour q
47			! radsol rayonnement net aus sol (LW + SW)
48			! dif_grnd coeff. diffusion vers le sol profond
49			!
50			! output:
51			! tsurf_new temperature au sol
52			! qsurf humidite de l'air au dessus du sol
53			! fluxsens flux de chaleur sensible
54			! fluxlat flux de chaleur latente
55			! dflux_s derivee du flux de chaleur sensible / Ts
56			! dflux_l derivee du flux de chaleur latente / Ts
57			! sens_prec_liq flux sensible li� aux echanges de precipitations liquides
58			! sens_prec_sol precipitations solides
59			! lat_prec_liq flux latent li� aux echanges de precipitations liquides
60			! lat_prec_sol precipitations solides
61
62			INCLUDE "YOETHF.h"
63			INCLUDE "FCTTRE.h"
64			INCLUDE "YOMCST.h"
65
66			! Parametres d'entree
67			!****************************************************************************************
68			INTEGER, INTENT(IN) :: knon, nisurf
69			REAL , INTENT(IN) :: dtime
70			REAL, DIMENSION(klon), INTENT(IN) :: petAcoef, peqAcoef
71			REAL, DIMENSION(klon), INTENT(IN) :: petBcoef, peqBcoef
72			REAL, DIMENSION(klon), INTENT(IN) :: ps, q1lay
73			REAL, DIMENSION(klon), INTENT(IN) :: tsurf, p1lay, cal, beta, cdragh,cdragq
74			REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
75			REAL, DIMENSION(klon), INTENT(IN) :: radsol, dif_grnd
76			REAL, DIMENSION(klon), INTENT(IN) :: t1lay, u1lay, v1lay,gustiness
77			REAL, INTENT(IN) :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else)
78
79			real, intent(in), optional:: rhoa(:) ! (knon)
80			! density of moist air (kg / m3)
81
82			! Parametres entree-sorties
83			!****************************************************************************************
84			REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! snow pas utile
85
86			! Parametres sorties
87			!****************************************************************************************
88			REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
89			REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new, evap, fluxsens, fluxlat
90			REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
91			REAL, intent(out), OPTIONAL:: sens_prec_liq(:), sens_prec_sol(:) ! (knon)
92			REAL, DIMENSION(klon), OPTIONAL :: lat_prec_liq, lat_prec_sol
93
94			! Variables locales
95			!****************************************************************************************
96			INTEGER :: i
97		2304	REAL, DIMENSION(klon) :: zx_mh, zx_nh, zx_oh
98		2304	REAL, DIMENSION(klon) :: zx_mq, zx_nq, zx_oq
99		2304	REAL, DIMENSION(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat
100		2304	REAL, DIMENSION(klon) :: zx_sl, zx_coefh, zx_coefq, zx_wind
101		2304	REAL, DIMENSION(klon) :: d_ts
102			REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
103			REAL :: qsat_new, q1_new
104			REAL, PARAMETER :: t_grnd = 271.35, t_coup = 273.15
105			REAL, PARAMETER :: max_eau_sol = 150.0
106			CHARACTER (len = 20) :: modname = 'calcul_fluxs'
107			LOGICAL :: fonte_neige
108			LOGICAL, SAVE :: check = .FALSE.
109			!$OMP THREADPRIVATE(check)
110
111			! End definition
112			!****************************************************************************************
113
114	✗✓	1152	IF (check) WRITE(,)'Entree ', modname,' surface = ',nisurf
115
116	✗✓	1152	IF (check) THEN
117			WRITE(,)' radsol (min, max)', &
118			MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
119			ENDIF
120
121			! Traitement neige et humidite du sol
122			!****************************************************************************************
123			!
124			!!$ WRITE(,)'test calcul_flux, surface ', nisurf
125			!!PB test
126			!!$ if (nisurf == is_oce) then
127			!!$ snow = 0.
128			!!$ qsol = max_eau_sol
129			!!$ else
130			!!$ where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
131			!!$ where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
132			!!$! snow = max(0.0, snow + (precip_snow - evap) * dtime)
133			!!$ where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
134			!!$ endif
135			!!$ IF (nisurf /= is_ter) qsol = max_eau_sol
136
137
138			!
139			! Initialisation
140			!****************************************************************************************
141	✓✓	1146240	evap = 0.
142	✓✓	1146240	fluxsens=0.
143	✓✓	1146240	fluxlat=0.
144	✓✓	1146240	dflux_s = 0.
145	✓✓	1146240	dflux_l = 0.
146	✓✓✓✓	573696	if (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
147	✓✓✓✓	573696	if (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
148			!
149			! zx_qs = qsat en kg/kg
150			!****************************************************************************************
151	✓✓	473954	DO i = 1, knon
152		472802	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
153			IF (thermcep) THEN
154		472802	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
155		472802	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
156		472802	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
157		472802	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
158		472802	zx_qs=MIN(0.5,zx_qs)
159		472802	zcor=1./(1.-retv*zx_qs)
160		472802	zx_qs=zx_qs*zcor
161			zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
162		472802	/RLVTT / zx_pkh(i)
163			ELSE
164			IF (tsurf(i).LT.t_coup) THEN
165			zx_qs = qsats(tsurf(i)) / ps(i)
166			zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
167			/ zx_pkh(i)
168			ELSE
169			zx_qs = qsatl(tsurf(i)) / ps(i)
170			zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
171			/ zx_pkh(i)
172			ENDIF
173			ENDIF
174		472802	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
175		472802	zx_qsat(i) = zx_qs
176		472802	zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)2+v1lay(i)2)
177		472802	zx_coefh(i) = cdragh(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
178		1152	zx_coefq(i) = cdragq(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
179			! zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)2+v1lay(i)2) &
180			! * p1lay(i)/(RD*t1lay(i))
181			! zx_coefh(i) = cdragh(i) * zx_wind(i)
182			! zx_coefq(i) = cdragq(i) * zx_wind(i)
183			ENDDO
184
185
186			! === Calcul de la temperature de surface ===
187			! zx_sl = chaleur latente d'evaporation ou de sublimation
188			!****************************************************************************************
189
190	✓✓	473954	DO i = 1, knon
191		472802	zx_sl(i) = RLVTT
192		1152	IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
193			ENDDO
194
195
196	✓✓	473954	DO i = 1, knon
197			! Q
198		472802	zx_oq(i) = 1. - (beta(i) * zx_coefq(i) * peqBcoef(i) * dtime)
199			zx_mq(i) = beta(i) * zx_coefq(i) * &
200			(peqAcoef(i) - &
201			! conv num avec precedente version
202			fqsat * zx_qsat(i) + fqsat * zx_dq_s_dt(i) * tsurf(i)) &
203			! fqsat * ( zx_qsat(i) - zx_dq_s_dt(i) * tsurf(i)) ) &
204		472802	/ zx_oq(i)
205			zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
206		472802	/ zx_oq(i)
207
208			! H
209		472802	zx_oh(i) = 1. - (zx_coefh(i) * petBcoef(i) * dtime)
210		472802	zx_mh(i) = zx_coefh(i) * petAcoef(i) / zx_oh(i)
211		472802	zx_nh(i) = - (zx_coefh(i) * RCPD * zx_pkh(i))/ zx_oh(i)
212
213			! Tsurface
214			tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
215			(radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
216			+ dif_grnd(i) * t_grnd * dtime)/ &
217			( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
218			zx_nh(i) + zx_sl(i) * zx_nq(i)) &
219		472802	+ dtime * dif_grnd(i))
220
221			!
222			! Y'a-t-il fonte de neige?
223			!
224			! fonte_neige = (nisurf /= is_oce) .AND. &
225			! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
226			! & .AND. (tsurf_new(i) >= RTT)
227			! if (fonte_neige) tsurf_new(i) = RTT
228		472802	d_ts(i) = tsurf_new(i) - tsurf(i)
229			! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
230			! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
231
232			!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
233			!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
234		472802	evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
235		472802	fluxlat(i) = - evap(i) * zx_sl(i)
236		472802	fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
237
238			! Derives des flux dF/dTs (W m-2 K-1):
239		472802	dflux_s(i) = zx_nh(i)
240		472802	dflux_l(i) = (zx_sl(i) * zx_nq(i))
241
242			! Nouvelle valeure de l'humidite au dessus du sol
243		472802	qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
244		472802	q1_new = peqAcoef(i) - peqBcoef(i)evap(i)dtime
245		472802	qsurf(i)=q1_new(1.-beta(i)) + beta(i)qsat_new
246			!
247			! en cas de fonte de neige
248			!
249			! if (fonte_neige) then
250			! bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
251			! & dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
252			! & RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
253			! bilan_f = max(0., bilan_f)
254			! fq_fonte = bilan_f / zx_sl(i)
255			! snow(i) = max(0., snow(i) - fq_fonte * dtime)
256			! qsol(i) = qsol(i) + (fq_fonte * dtime)
257			! endif
258			!!$ if (nisurf == is_ter) &
259			!!$ & run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
260			!!$ qsol(i) = min(qsol(i), max_eau_sol)
261			!
262			! calcul de l'enthalpie des precipitations liquides et solides
263	✓✗	280418	if (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
264			= - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
265		280418	q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
266	✓✓✓✗	472802	if (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0.
267			! On calcule par rapport a T=0
268			!! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i)
269			!! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i)
270
271	✓✓	472802	if (PRESENT(lat_prec_liq)) &
272		280418	lat_prec_liq(i) = precip_rain(i) * (RLVTT - RLVTT)
273	✓✓	472802	if (PRESENT(lat_prec_sol)) &
274		1152	lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
275			ENDDO
276
277
278			!**************************************************************************
279			!
280	✓✓✓✓	1728	END SUBROUTINE calcul_fluxs
281			!
282			!****************************************************************************************
283			!
284		1152	SUBROUTINE calcul_flux_wind(knon, dtime, &
285			u0, v0, u1, v1, gustiness, cdrag_m, &
286			AcoefU, AcoefV, BcoefU, BcoefV, &
287			p1lay, t1lay, &
288			flux_u1, flux_v1)
289
290	✓✓✓✓	1979242	USE dimphy
291			INCLUDE "YOMCST.h"
292			INCLUDE "clesphys.h"
293
294			! Input arguments
295			!****************************************************************************************
296			INTEGER, INTENT(IN) :: knon
297			REAL, INTENT(IN) :: dtime
298			REAL, DIMENSION(klon), INTENT(IN) :: u0, v0 ! u and v at niveau 0
299			REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! u and v at niveau 1
300			REAL, DIMENSION(klon), INTENT(IN) :: cdrag_m ! cdrag pour momentum
301			REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
302			REAL, DIMENSION(klon), INTENT(IN) :: p1lay ! pression 1er niveau (milieu de couche)
303			REAL, DIMENSION(klon), INTENT(IN) :: t1lay ! temperature
304			! Output arguments
305			!****************************************************************************************
306			REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1
307			REAL, DIMENSION(klon), INTENT(OUT) :: flux_v1
308
309			! Local variables
310			!****************************************************************************************
311			INTEGER :: i
312			REAL :: mod_wind, buf
313
314			!****************************************************************************************
315			! Calculate the surface flux
316			!
317			!****************************************************************************************
318	✓✓	473954	DO i=1,knon
319		472802	mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))2 + (v1(i)-v0(i))2)
320		472802	buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
321		472802	flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
322		473954	flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
323			END DO
324
325		1152	END SUBROUTINE calcul_flux_wind
326			!
327			!****************************************************************************************
328			!
329			END MODULE calcul_fluxs_mod