GCC Code Coverage Report


Directory: ./
File: phys/calcul_fluxs_mod.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 68 70 97.1%
Branches: 56 86 65.1%

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 790372 SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
10 tsurf, p1lay, cal, beta, cdragh, cdragq, ps, &
11 precip_rain, precip_snow, snow, qsurf, &
12 1920 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
12/18
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 960 times.
✓ Branch 2 taken 960 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 960 times.
✓ Branch 6 taken 960 times.
✓ Branch 7 taken 960 times.
✓ Branch 8 taken 960 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 960 times.
✓ Branch 12 taken 960 times.
✓ Branch 13 taken 960 times.
✓ Branch 14 taken 960 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 960 times.
1920 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 3840 REAL, DIMENSION(klon) :: zx_mh, zx_nh, zx_oh
98 3840 REAL, DIMENSION(klon) :: zx_mq, zx_nq, zx_oq
99 3840 REAL, DIMENSION(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat
100 3840 REAL, DIMENSION(klon) :: zx_sl, zx_coefh, zx_coefq, zx_wind
101 3840 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
1920 IF (check) WRITE(*,*)'Entree ', modname,' surface = ',nisurf
115
116
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1920 times.
1920 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
2/2
✓ Branch 0 taken 1908480 times.
✓ Branch 1 taken 1920 times.
1910400 evap = 0.
142
2/2
✓ Branch 0 taken 1908480 times.
✓ Branch 1 taken 1920 times.
1910400 fluxsens=0.
143
2/2
✓ Branch 0 taken 1908480 times.
✓ Branch 1 taken 1920 times.
1910400 fluxlat=0.
144
2/2
✓ Branch 0 taken 1908480 times.
✓ Branch 1 taken 1920 times.
1910400 dflux_s = 0.
145
2/2
✓ Branch 0 taken 1908480 times.
✓ Branch 1 taken 1920 times.
1910400 dflux_l = 0.
146
4/4
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 960 times.
✓ Branch 2 taken 954240 times.
✓ Branch 3 taken 960 times.
956160 if (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
147
4/4
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 960 times.
✓ Branch 2 taken 954240 times.
✓ Branch 3 taken 960 times.
956160 if (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
148 !
149 ! zx_qs = qsat en kg/kg
150 !****************************************************************************************
151
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
790372 DO i = 1, knon
152 788452 zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
153 IF (thermcep) THEN
154 788452 zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
155 788452 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
156 788452 zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
157 788452 zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
158 788452 zx_qs=MIN(0.5,zx_qs)
159 788452 zcor=1./(1.-retv*zx_qs)
160 788452 zx_qs=zx_qs*zcor
161 zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
162 788452 /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 788452 zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
175 788452 zx_qsat(i) = zx_qs
176 788452 zx_wind(i)=min_wind_speed+SQRT(gustiness(i)+u1lay(i)**2+v1lay(i)**2)
177 788452 zx_coefh(i) = cdragh(i) * zx_wind(i) * p1lay(i)/(RD*t1lay(i))
178 1920 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
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
790372 DO i = 1, knon
191 788452 zx_sl(i) = RLVTT
192 1920 IF (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
193 ENDDO
194
195
196
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
790372 DO i = 1, knon
197 ! Q
198 788452 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 788452 / zx_oq(i)
205 zx_nq(i) = beta(i) * zx_coefq(i) * (- fqsat * zx_dq_s_dt(i)) &
206 788452 / zx_oq(i)
207
208 ! H
209 788452 zx_oh(i) = 1. - (zx_coefh(i) * petBcoef(i) * dtime)
210 788452 zx_mh(i) = zx_coefh(i) * petAcoef(i) / zx_oh(i)
211 788452 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 788452 + 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 788452 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 788452 evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
235 788452 fluxlat(i) = - evap(i) * zx_sl(i)
236 788452 fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
237
238 ! Derives des flux dF/dTs (W m-2 K-1):
239 788452 dflux_s(i) = zx_nh(i)
240 788452 dflux_l(i) = (zx_sl(i) * zx_nq(i))
241
242 ! Nouvelle valeure de l'humidite au dessus du sol
243 788452 qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
244 788452 q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
245 788452 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
1/2
✓ Branch 0 taken 467812 times.
✗ Branch 1 not taken.
467812 if (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
264 = - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
265 467812 q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
266
3/4
✓ Branch 0 taken 467812 times.
✓ Branch 1 taken 320640 times.
✓ Branch 2 taken 467812 times.
✗ Branch 3 not taken.
788452 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
2/2
✓ Branch 0 taken 467812 times.
✓ Branch 1 taken 320640 times.
788452 if (PRESENT(lat_prec_liq)) &
272 467812 lat_prec_liq(i) = precip_rain(i) * (RLVTT - RLVTT)
273
2/2
✓ Branch 0 taken 467812 times.
✓ Branch 1 taken 320640 times.
788452 if (PRESENT(lat_prec_sol)) &
274 1920 lat_prec_sol(i) = precip_snow(i) * (RLSTT - RLVTT)
275 ENDDO
276
277
278 !**************************************************************************
279 !
280
4/4
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 960 times.
✓ Branch 2 taken 960 times.
✓ Branch 3 taken 960 times.
2880 END SUBROUTINE calcul_fluxs
281 !
282 !****************************************************************************************
283 !
284 1920 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
4/4
✓ Branch 0 taken 289280 times.
✓ Branch 1 taken 499172 times.
✓ Branch 2 taken 467812 times.
✓ Branch 3 taken 320640 times.
3300980 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
2/2
✓ Branch 0 taken 788452 times.
✓ Branch 1 taken 1920 times.
790372 DO i=1,knon
319 788452 mod_wind = min_wind_speed + SQRT(gustiness(i)+(u1(i) - u0(i))**2 + (v1(i)-v0(i))**2)
320 788452 buf = cdrag_m(i) * mod_wind * p1lay(i)/(RD*t1lay(i))
321 788452 flux_u1(i) = (AcoefU(i) - u0(i)) / (1/buf - BcoefU(i)*dtime )
322 790372 flux_v1(i) = (AcoefV(i) - v0(i)) / (1/buf - BcoefV(i)*dtime )
323 END DO
324
325 1920 END SUBROUTINE calcul_flux_wind
326 !
327 !****************************************************************************************
328 !
329 END MODULE calcul_fluxs_mod
330