GCC Code Coverage Report

Directory: ./
File: phys/cloudth_mod.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 173 603 28.7%
Branches: 74 302 24.5%
Line Branch Exec Source
1 MODULE cloudth_mod
2
3 IMPLICIT NONE
4
5 CONTAINS
6
7 ✗ SUBROUTINE cloudth(ngrid,klev,ind2, &
8 ✗ & ztv,po,zqta,fraca, &
9 & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
10 & ratqs,zqs,t)
11
12
13 IMPLICIT NONE
14
15
16 !===========================================================================
17 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
18 ! Date : 25 Mai 2010
19 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
20 !===========================================================================
21
22
23 !
24 ! $Header$
25 !
26 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
27 ! veillez � n'utiliser que des ! pour les commentaires
28 ! et � bien positionner les & des lignes de continuation
29 ! (les placer en colonne 6 et en colonne 73)
30 !
31 !
32 ! A1.0 Fundamental constants
33 REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
34 ! A1.1 Astronomical constants
35 REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
36 ! A1.1.bis Constantes concernant l'orbite de la Terre:
37 REAL R_ecc, R_peri, R_incl
38 ! A1.2 Geoide
39 REAL RA,RG,R1SA
40 ! A1.3 Radiation
41 ! REAL RSIGMA,RI0
42 REAL RSIGMA
43 ! A1.4 Thermodynamic gas phase
44 REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
45 REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
46 REAL RKAPPA,RETV, eps_w
47 ! A1.5,6 Thermodynamic liquid,solid phases
48 REAL RCW,RCS
49 ! A1.7 Thermodynamic transition of phase
50 REAL RLVTT,RLSTT,RLMLT,RTT,RATM
51 ! A1.8 Curve of saturation
52 REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
53 REAL RALPD,RBETD,RGAMD
54 !
55 COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO &
56 & ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA &
57 & ,R_ecc, R_peri, R_incl &
58 & ,RA ,RG ,R1SA &
59 & ,RSIGMA &
60 & ,R ,RMD ,RMV ,RD ,RV ,RCPD &
61 & ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 &
62 & ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w &
63 & ,RCW ,RCS &
64 & ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM &
65 & ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS &
66 & ,RALPD ,RBETD ,RGAMD
67 ! ------------------------------------------------------------------
68 !$OMP THREADPRIVATE(/YOMCST/)
69 !
70 ! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $
71 !
72 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
73 ! veillez n'utiliser que des ! pour les commentaires
74 ! et bien positionner les & des lignes de continuation
75 ! (les placer en colonne 6 et en colonne 73)
76 !
77 !* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS
78 !
79 ! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION
80 ! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR
81 ! ICE(*R_IES*).
82 ! *RVTMP2* *RVTMP2=RCPV/RCPD-1.
83 ! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.)
84 !
85 REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES
86 REAL RVTMP2, RHOH2O
87 REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
88 REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
89 LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
90 ! If FALSE, then variables set by suphel.F90
91 COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, &
92 & RVTMP2, RHOH2O, &
93 & R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, &
94 & RALFDCP,RTWAT,RTBER,RTBERCU, &
95 & RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
96 & RKOOP2, &
97 & OK_BAD_ECMWF_THERMO
98
99 !$OMP THREADPRIVATE(/YOETHF/)
100 !
101 ! $Header$
102 !
103 !
104 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
105 ! veillez n'utiliser que des ! pour les commentaires
106 ! et bien positionner les & des lignes de continuation
107 ! (les placer en colonne 6 et en colonne 73)
108 !
109 ! ------------------------------------------------------------------
110 ! This COMDECK includes the Thermodynamical functions for the cy39
111 ! ECMWF Physics package.
112 ! Consistent with YOMCST Basic physics constants, assuming the
113 ! partial pressure of water vapour is given by a first order
114 ! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants
115 ! in YOETHF
116 ! ------------------------------------------------------------------
117 REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG
118 REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl
119 LOGICAL thermcep
120 PARAMETER (thermcep=.TRUE.)
121 !
122 FOEEW ( PTARG,PDELARG ) = EXP ( &
123 & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
124 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
125 !
126 FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG &
127 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
128 !
129 qsats(ptarg) = 100.0 * 0.622 * 10.0 &
130 & ** (2.07023 - 0.00320991 * ptarg &
131 & - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
132 qsatl(ptarg) = 100.0 * 0.622 * 10.0 &
133 & ** (23.8319 - 2948.964 / ptarg &
134 & - 5.028 * LOG10(ptarg) &
135 & - 29810.16 * EXP( - 0.0699382 * ptarg) &
136 & + 25.21935 * EXP( - 2999.924 / ptarg))
137 !
138 dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
139 & +2484.896*LOG(10.)/ptarg**2 &
140 & -0.00320991*LOG(10.))
141 dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
142 & (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
143 & +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
144 & +29810.16*0.0699382*EXP(-0.0699382*ptarg))
145 integer :: iflag_thermals,nsplit_thermals
146
147 !!! nrlmd le 10/04/2012
148 integer :: iflag_trig_bl,iflag_clos_bl
149 integer :: tau_trig_shallow,tau_trig_deep
150 real :: s_trig
151 !!! fin nrlmd le 10/04/2012
152
153 real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30.
154 real :: alp_bl_k
155 real :: tau_thermals,fact_thermals_ed_dz
156 integer,parameter :: w2di_thermals=0
157 integer :: isplit
158
159 integer :: iflag_coupl,iflag_clos,iflag_wake
160 integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure
161
162 common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure
163 common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz
164 common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
165 common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
166
167 !!! nrlmd le 10/04/2012
168 common/ctherm6/iflag_trig_bl,iflag_clos_bl
169 common/ctherm7/tau_trig_shallow,tau_trig_deep
170 common/ctherm8/s_trig
171 !!! fin nrlmd le 10/04/2012
172
173 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
174 !$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
175 !
176 ! $Id: nuage.h 2945 2017-07-12 14:20:24Z jbmadeleine $
177 !
178 REAL rad_froid, rad_chau1, rad_chau2, t_glace_max, t_glace_min
179 REAL exposant_glace
180 REAL rei_min,rei_max
181 REAL tau_cld_cv,coefw_cld_cv
182
183 REAL tmax_fonte_cv
184
185 INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
186 INTEGER iflag_rain_incloud_vol
187
188 common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max, &
189 & t_glace_min,exposant_glace,rei_min,rei_max, &
190 & tau_cld_cv,coefw_cld_cv, &
191 & tmax_fonte_cv, &
192 & iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv, &
193 & iflag_rain_incloud_vol
194 !$OMP THREADPRIVATE(/nuagecom/)
195
196 INTEGER itap,ind1,ind2
197 INTEGER ngrid,klev,klon,l,ig
198
199 REAL ztv(ngrid,klev)
200 REAL po(ngrid)
201 ✗ REAL zqenv(ngrid)
202 REAL zqta(ngrid,klev)
203
204 REAL fraca(ngrid,klev+1)
205 REAL zpspsk(ngrid,klev)
206 REAL paprs(ngrid,klev+1)
207 REAL pplay(ngrid,klev)
208 REAL ztla(ngrid,klev)
209 REAL zthl(ngrid,klev)
210
211 ✗ REAL zqsatth(ngrid,klev)
212 ✗ REAL zqsatenv(ngrid,klev)
213
214
215 ✗ REAL sigma1(ngrid,klev)
216 ✗ REAL sigma2(ngrid,klev)
217 ✗ REAL qlth(ngrid,klev)
218 ✗ REAL qlenv(ngrid,klev)
219 ✗ REAL qltot(ngrid,klev)
220 ✗ REAL cth(ngrid,klev)
221 ✗ REAL cenv(ngrid,klev)
222 REAL ctot(ngrid,klev)
223 ✗ REAL rneb(ngrid,klev)
224 REAL t(ngrid,klev)
225 REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
226 REAL rdd,cppd,Lv
227 REAL alth,alenv,ath,aenv
228 REAL sth,senv,sigma1s,sigma2s,xth,xenv
229 REAL Tbef,zdelta,qsatbef,zcor
230 REAL qlbef
231 REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
232
233 REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
234 REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
235 REAL zqs(ngrid), qcloud(ngrid)
236 REAL erf
237
238
239
240
241 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
242 ! Gestion de deux versions de cloudth
243 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
244
245 ✗ IF (iflag_cloudth_vert.GE.1) THEN
246 CALL cloudth_vert(ngrid,klev,ind2, &
247 & ztv,po,zqta,fraca, &
248 & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
249 ✗ & ratqs,zqs,t)
250 ✗ RETURN
251 ENDIF
252 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
253
254
255 !-------------------------------------------------------------------------------
256 ! Initialisation des variables r?elles
257 !-------------------------------------------------------------------------------
258 ✗ sigma1(:,:)=0.
259 ✗ sigma2(:,:)=0.
260 ✗ qlth(:,:)=0.
261 ✗ qlenv(:,:)=0.
262 ✗ qltot(:,:)=0.
263 ✗ rneb(:,:)=0.
264 ✗ qcloud(:)=0.
265 ✗ cth(:,:)=0.
266 ✗ cenv(:,:)=0.
267 ✗ ctot(:,:)=0.
268 qsatmmussig1=0.
269 qsatmmussig2=0.
270 rdd=287.04
271 cppd=1005.7
272 pi=3.14159
273 Lv=2.5e6
274 sqrt2pi=sqrt(2.*pi)
275
276
277
278 !-------------------------------------------------------------------------------
279 ! Calcul de la fraction du thermique et des ?cart-types des distributions
280 !-------------------------------------------------------------------------------
281 ✗ do ind1=1,ngrid
282
283 ✗ if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
284
285 ✗ zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
286
287
288 ! zqenv(ind1)=po(ind1)
289 ✗ Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
290 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
291 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
292 ✗ qsatbef=MIN(0.5,qsatbef)
293 ✗ zcor=1./(1.-retv*qsatbef)
294 ✗ qsatbef=qsatbef*zcor
295 ✗ zqsatenv(ind1,ind2)=qsatbef
296
297
298
299
300 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
301 ✗ aenv=1./(1.+(alenv*Lv/cppd))
302 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
303
304
305
306
307 ✗ Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
308 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
309 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
310 ✗ qsatbef=MIN(0.5,qsatbef)
311 ✗ zcor=1./(1.-retv*qsatbef)
312 ✗ qsatbef=qsatbef*zcor
313 ✗ zqsatth(ind1,ind2)=qsatbef
314
315 ✗ alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)
316 ✗ ath=1./(1.+(alth*Lv/cppd))
317 ✗ sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
318
319
320
321 !------------------------------------------------------------------------------
322 ! Calcul des ?cart-types pour s
323 !------------------------------------------------------------------------------
324
325 ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
326 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
327 ! if (paprs(ind1,ind2).gt.90000) then
328 ! ratqs(ind1,ind2)=0.002
329 ! else
330 ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
331 ! endif
332 ✗ sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
333 ✗ sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
334 ! sigma1s=ratqs(ind1,ind2)*po(ind1)
335 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003
336
337 !------------------------------------------------------------------------------
338 ! Calcul de l'eau condens?e et de la couverture nuageuse
339 !------------------------------------------------------------------------------
340 sqrt2pi=sqrt(2.*pi)
341 ✗ xth=sth/(sqrt(2.)*sigma2s)
342 ✗ xenv=senv/(sqrt(2.)*sigma1s)
343 ✗ cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
344 ✗ cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
345 ✗ ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
346
347 ✗ qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
348 ✗ qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
349 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
350
351 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
352 ✗ if (ctot(ind1,ind2).lt.1.e-10) then
353 ✗ ctot(ind1,ind2)=0.
354 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
355
356 else
357
358 ctot(ind1,ind2)=ctot(ind1,ind2)
359 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
360
361 endif
362
363
364 ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
365
366
367 else ! gaussienne environnement seule
368
369 ✗ zqenv(ind1)=po(ind1)
370 ✗ Tbef=t(ind1,ind2)
371 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
372 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
373 ✗ qsatbef=MIN(0.5,qsatbef)
374 ✗ zcor=1./(1.-retv*qsatbef)
375 ✗ qsatbef=qsatbef*zcor
376 ✗ zqsatenv(ind1,ind2)=qsatbef
377
378
379 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
380 ✗ zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
381 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
382 ✗ aenv=1./(1.+(alenv*Lv/cppd))
383 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
384
385
386 ✗ sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
387
388 sqrt2pi=sqrt(2.*pi)
389 ✗ xenv=senv/(sqrt(2.)*sigma1s)
390 ✗ ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
391 ✗ qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
392
393 ✗ if (ctot(ind1,ind2).lt.1.e-3) then
394 ✗ ctot(ind1,ind2)=0.
395 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
396
397 else
398
399 ctot(ind1,ind2)=ctot(ind1,ind2)
400 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
401
402 endif
403
404
405
406
407
408
409 endif
410 enddo
411
412 return
413 ! end
414 END SUBROUTINE cloudth
415
416
417
418 !===========================================================================
419 ✗ SUBROUTINE cloudth_vert(ngrid,klev,ind2, &
420 & ztv,po,zqta,fraca, &
421 ✗ & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
422 & ratqs,zqs,t)
423
424 !===========================================================================
425 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
426 ! Date : 25 Mai 2010
427 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
428 !===========================================================================
429
430
431 ✗ USE ioipsl_getin_p_mod, ONLY : getin_p
432
433 IMPLICIT NONE
434
435 !
436 ! $Header$
437 !
438 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
439 ! veillez � n'utiliser que des ! pour les commentaires
440 ! et � bien positionner les & des lignes de continuation
441 ! (les placer en colonne 6 et en colonne 73)
442 !
443 !
444 ! A1.0 Fundamental constants
445 REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
446 ! A1.1 Astronomical constants
447 REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
448 ! A1.1.bis Constantes concernant l'orbite de la Terre:
449 REAL R_ecc, R_peri, R_incl
450 ! A1.2 Geoide
451 REAL RA,RG,R1SA
452 ! A1.3 Radiation
453 ! REAL RSIGMA,RI0
454 REAL RSIGMA
455 ! A1.4 Thermodynamic gas phase
456 REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
457 REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
458 REAL RKAPPA,RETV, eps_w
459 ! A1.5,6 Thermodynamic liquid,solid phases
460 REAL RCW,RCS
461 ! A1.7 Thermodynamic transition of phase
462 REAL RLVTT,RLSTT,RLMLT,RTT,RATM
463 ! A1.8 Curve of saturation
464 REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
465 REAL RALPD,RBETD,RGAMD
466 !
467 COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO &
468 & ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA &
469 & ,R_ecc, R_peri, R_incl &
470 & ,RA ,RG ,R1SA &
471 & ,RSIGMA &
472 & ,R ,RMD ,RMV ,RD ,RV ,RCPD &
473 & ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 &
474 & ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w &
475 & ,RCW ,RCS &
476 & ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM &
477 & ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS &
478 & ,RALPD ,RBETD ,RGAMD
479 ! ------------------------------------------------------------------
480 !$OMP THREADPRIVATE(/YOMCST/)
481 !
482 ! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $
483 !
484 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
485 ! veillez n'utiliser que des ! pour les commentaires
486 ! et bien positionner les & des lignes de continuation
487 ! (les placer en colonne 6 et en colonne 73)
488 !
489 !* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS
490 !
491 ! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION
492 ! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR
493 ! ICE(*R_IES*).
494 ! *RVTMP2* *RVTMP2=RCPV/RCPD-1.
495 ! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.)
496 !
497 REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES
498 REAL RVTMP2, RHOH2O
499 REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
500 REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
501 LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
502 ! If FALSE, then variables set by suphel.F90
503 COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, &
504 & RVTMP2, RHOH2O, &
505 & R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, &
506 & RALFDCP,RTWAT,RTBER,RTBERCU, &
507 & RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
508 & RKOOP2, &
509 & OK_BAD_ECMWF_THERMO
510
511 !$OMP THREADPRIVATE(/YOETHF/)
512 !
513 ! $Header$
514 !
515 !
516 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
517 ! veillez n'utiliser que des ! pour les commentaires
518 ! et bien positionner les & des lignes de continuation
519 ! (les placer en colonne 6 et en colonne 73)
520 !
521 ! ------------------------------------------------------------------
522 ! This COMDECK includes the Thermodynamical functions for the cy39
523 ! ECMWF Physics package.
524 ! Consistent with YOMCST Basic physics constants, assuming the
525 ! partial pressure of water vapour is given by a first order
526 ! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants
527 ! in YOETHF
528 ! ------------------------------------------------------------------
529 REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG
530 REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl
531 LOGICAL thermcep
532 PARAMETER (thermcep=.TRUE.)
533 !
534 FOEEW ( PTARG,PDELARG ) = EXP ( &
535 & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
536 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
537 !
538 FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG &
539 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
540 !
541 qsats(ptarg) = 100.0 * 0.622 * 10.0 &
542 & ** (2.07023 - 0.00320991 * ptarg &
543 & - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
544 qsatl(ptarg) = 100.0 * 0.622 * 10.0 &
545 & ** (23.8319 - 2948.964 / ptarg &
546 & - 5.028 * LOG10(ptarg) &
547 & - 29810.16 * EXP( - 0.0699382 * ptarg) &
548 & + 25.21935 * EXP( - 2999.924 / ptarg))
549 !
550 dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
551 & +2484.896*LOG(10.)/ptarg**2 &
552 & -0.00320991*LOG(10.))
553 dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
554 & (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
555 & +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
556 & +29810.16*0.0699382*EXP(-0.0699382*ptarg))
557 integer :: iflag_thermals,nsplit_thermals
558
559 !!! nrlmd le 10/04/2012
560 integer :: iflag_trig_bl,iflag_clos_bl
561 integer :: tau_trig_shallow,tau_trig_deep
562 real :: s_trig
563 !!! fin nrlmd le 10/04/2012
564
565 real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30.
566 real :: alp_bl_k
567 real :: tau_thermals,fact_thermals_ed_dz
568 integer,parameter :: w2di_thermals=0
569 integer :: isplit
570
571 integer :: iflag_coupl,iflag_clos,iflag_wake
572 integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure
573
574 common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure
575 common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz
576 common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
577 common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
578
579 !!! nrlmd le 10/04/2012
580 common/ctherm6/iflag_trig_bl,iflag_clos_bl
581 common/ctherm7/tau_trig_shallow,tau_trig_deep
582 common/ctherm8/s_trig
583 !!! fin nrlmd le 10/04/2012
584
585 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
586 !$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
587 !
588 ! $Id: nuage.h 2945 2017-07-12 14:20:24Z jbmadeleine $
589 !
590 REAL rad_froid, rad_chau1, rad_chau2, t_glace_max, t_glace_min
591 REAL exposant_glace
592 REAL rei_min,rei_max
593 REAL tau_cld_cv,coefw_cld_cv
594
595 REAL tmax_fonte_cv
596
597 INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
598 INTEGER iflag_rain_incloud_vol
599
600 common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max, &
601 & t_glace_min,exposant_glace,rei_min,rei_max, &
602 & tau_cld_cv,coefw_cld_cv, &
603 & tmax_fonte_cv, &
604 & iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv, &
605 & iflag_rain_incloud_vol
606 !$OMP THREADPRIVATE(/nuagecom/)
607
608 INTEGER itap,ind1,ind2
609 INTEGER ngrid,klev,klon,l,ig
610
611 REAL ztv(ngrid,klev)
612 REAL po(ngrid)
613 ✗ REAL zqenv(ngrid)
614 REAL zqta(ngrid,klev)
615
616 REAL fraca(ngrid,klev+1)
617 REAL zpspsk(ngrid,klev)
618 REAL paprs(ngrid,klev+1)
619 REAL pplay(ngrid,klev)
620 REAL ztla(ngrid,klev)
621 REAL zthl(ngrid,klev)
622
623 ✗ REAL zqsatth(ngrid,klev)
624 ✗ REAL zqsatenv(ngrid,klev)
625
626
627 ✗ REAL sigma1(ngrid,klev)
628 ✗ REAL sigma2(ngrid,klev)
629 ✗ REAL qlth(ngrid,klev)
630 ✗ REAL qlenv(ngrid,klev)
631 ✗ REAL qltot(ngrid,klev)
632 ✗ REAL cth(ngrid,klev)
633 ✗ REAL cenv(ngrid,klev)
634 REAL ctot(ngrid,klev)
635 ✗ REAL rneb(ngrid,klev)
636 REAL t(ngrid,klev)
637 REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
638 REAL rdd,cppd,Lv,sqrt2,sqrtpi
639 REAL alth,alenv,ath,aenv
640 REAL sth,senv,sigma1s,sigma2s,xth,xenv
641 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
642 REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
643 REAL Tbef,zdelta,qsatbef,zcor
644 REAL qlbef
645 REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
646 ! Change the width of the PDF used for vertical subgrid scale heterogeneity
647 ! (J Jouhaud, JL Dufresne, JB Madeleine)
648 REAL,SAVE :: vert_alpha
649 !$OMP THREADPRIVATE(vert_alpha)
650 LOGICAL, SAVE :: firstcall = .TRUE.
651 !$OMP THREADPRIVATE(firstcall)
652
653 REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
654 REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
655 REAL zqs(ngrid), qcloud(ngrid)
656 REAL erf
657
658 !------------------------------------------------------------------------------
659 ! Initialisation des variables r?elles
660 !------------------------------------------------------------------------------
661 ✗ sigma1(:,:)=0.
662 ✗ sigma2(:,:)=0.
663 ✗ qlth(:,:)=0.
664 ✗ qlenv(:,:)=0.
665 ✗ qltot(:,:)=0.
666 ✗ rneb(:,:)=0.
667 ✗ qcloud(:)=0.
668 ✗ cth(:,:)=0.
669 ✗ cenv(:,:)=0.
670 ✗ ctot(:,:)=0.
671 qsatmmussig1=0.
672 qsatmmussig2=0.
673 rdd=287.04
674 cppd=1005.7
675 pi=3.14159
676 Lv=2.5e6
677 sqrt2pi=sqrt(2.*pi)
678 sqrt2=sqrt(2.)
679 sqrtpi=sqrt(pi)
680
681 ✗ IF (firstcall) THEN
682 ✗ vert_alpha=0.5
683 ✗ CALL getin_p('cloudth_vert_alpha',vert_alpha)
684 ✗ WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
685 ✗ firstcall=.FALSE.
686 ENDIF
687
688 !-------------------------------------------------------------------------------
689 ! Calcul de la fraction du thermique et des ?cart-types des distributions
690 !-------------------------------------------------------------------------------
691 ✗ do ind1=1,ngrid
692
693 ✗ if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
694
695 ✗ zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
696
697
698 ! zqenv(ind1)=po(ind1)
699 ✗ Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
700 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
701 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
702 ✗ qsatbef=MIN(0.5,qsatbef)
703 ✗ zcor=1./(1.-retv*qsatbef)
704 ✗ qsatbef=qsatbef*zcor
705 ✗ zqsatenv(ind1,ind2)=qsatbef
706
707
708
709
710 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
711 ✗ aenv=1./(1.+(alenv*Lv/cppd))
712 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
713
714
715
716
717 ✗ Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
718 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
719 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
720 ✗ qsatbef=MIN(0.5,qsatbef)
721 ✗ zcor=1./(1.-retv*qsatbef)
722 ✗ qsatbef=qsatbef*zcor
723 ✗ zqsatth(ind1,ind2)=qsatbef
724
725 ✗ alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)
726 ✗ ath=1./(1.+(alth*Lv/cppd))
727 ✗ sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
728
729
730
731 !------------------------------------------------------------------------------
732 ! Calcul des ?cart-types pour s
733 !------------------------------------------------------------------------------
734
735 ✗ sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
736 ✗ sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
737 ! if (paprs(ind1,ind2).gt.90000) then
738 ! ratqs(ind1,ind2)=0.002
739 ! else
740 ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
741 ! endif
742 ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
743 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
744 ! sigma1s=ratqs(ind1,ind2)*po(ind1)
745 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003
746
747 !------------------------------------------------------------------------------
748 ! Calcul de l'eau condens?e et de la couverture nuageuse
749 !------------------------------------------------------------------------------
750 sqrt2pi=sqrt(2.*pi)
751 ✗ xth=sth/(sqrt(2.)*sigma2s)
752 ✗ xenv=senv/(sqrt(2.)*sigma1s)
753 ✗ cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
754 ✗ cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
755 ✗ ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
756
757 ✗ qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
758 ✗ qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
759 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
760
761 ✗ IF (iflag_cloudth_vert == 1) THEN
762 !-------------------------------------------------------------------------------
763 ! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
764 !-------------------------------------------------------------------------------
765 ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
766 ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
767 ✗ deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
768 ✗ deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
769 ! deltasenv=aenv*0.01*po(ind1)
770 ! deltasth=ath*0.01*zqta(ind1,ind2)
771 ✗ xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
772 ✗ xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
773 ✗ xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
774 ✗ xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
775 ✗ coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
776 ✗ coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
777
778 ✗ cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
779 ✗ cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
780 ✗ ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
781
782 ✗ IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
783 ✗ IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
784 ✗ IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
785 ✗ IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
786
787 ✗ qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
788 ! qlenv(ind1,ind2)=IntJ
789 ! print*, qlenv(ind1,ind2),'VERIF EAU'
790
791
792 ✗ IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
793 ! IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
794 ! IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
795 ✗ IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
796 ✗ IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
797 ✗ IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
798 ✗ qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
799 ! qlth(ind1,ind2)=IntJ
800 ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
801 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
802
803 ✗ ELSE IF (iflag_cloudth_vert == 2) THEN
804
805 !-------------------------------------------------------------------------------
806 ! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
807 !-------------------------------------------------------------------------------
808 ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
809 ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
810 ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
811 ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
812 ✗ deltasenv=aenv*vert_alpha*sigma1s
813 ✗ deltasth=ath*vert_alpha*sigma2s
814
815 ✗ xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
816 ✗ xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
817 ✗ xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
818 ✗ xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
819 ! coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
820 ! coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
821
822 ✗ cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
823 ✗ cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
824 ✗ ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
825
826 ✗ IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
827 ✗ IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
828 ✗ IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
829 ✗ IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
830
831 ! IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
832 ! IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
833 ! IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
834
835 ✗ qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
836 ! qlenv(ind1,ind2)=IntJ
837 ! print*, qlenv(ind1,ind2),'VERIF EAU'
838
839 ✗ IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
840 ✗ IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
841 ✗ IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
842 ✗ IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
843
844 ✗ qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
845 ! qlth(ind1,ind2)=IntJ
846 ! print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
847 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
848
849
850
851
852 ENDIF ! of if (iflag_cloudth_vert==1 or 2)
853
854 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
855
856 ✗ if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
857 ✗ ctot(ind1,ind2)=0.
858 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
859
860 else
861
862 ✗ ctot(ind1,ind2)=ctot(ind1,ind2)
863 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
864 ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
865 ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
866
867 endif
868
869
870
871 ! print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
872
873
874 else ! gaussienne environnement seule
875
876 ✗ zqenv(ind1)=po(ind1)
877 ✗ Tbef=t(ind1,ind2)
878 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
879 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
880 ✗ qsatbef=MIN(0.5,qsatbef)
881 ✗ zcor=1./(1.-retv*qsatbef)
882 ✗ qsatbef=qsatbef*zcor
883 ✗ zqsatenv(ind1,ind2)=qsatbef
884
885
886 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
887 ✗ zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
888 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
889 ✗ aenv=1./(1.+(alenv*Lv/cppd))
890 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
891
892
893 ✗ sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
894
895 sqrt2pi=sqrt(2.*pi)
896 ✗ xenv=senv/(sqrt(2.)*sigma1s)
897 ✗ ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
898 ✗ qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
899
900 ✗ if (ctot(ind1,ind2).lt.1.e-3) then
901 ✗ ctot(ind1,ind2)=0.
902 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
903
904 else
905
906 ctot(ind1,ind2)=ctot(ind1,ind2)
907 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
908
909 endif
910
911
912
913
914
915
916 endif
917 enddo
918
919 ✗ return
920 ! end
921 END SUBROUTINE cloudth_vert
922
923
924
925
926 18720 SUBROUTINE cloudth_v3(ngrid,klev,ind2, &
927 18720 & ztv,po,zqta,fraca, &
928 & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
929 & ratqs,zqs,t)
930
931
932 IMPLICIT NONE
933
934
935 !===========================================================================
936 ! Author : Arnaud Octavio Jam (LMD/CNRS)
937 ! Date : 25 Mai 2010
938 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
939 !===========================================================================
940
941
942 !
943 ! $Header$
944 !
945 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
946 ! veillez � n'utiliser que des ! pour les commentaires
947 ! et � bien positionner les & des lignes de continuation
948 ! (les placer en colonne 6 et en colonne 73)
949 !
950 !
951 ! A1.0 Fundamental constants
952 REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
953 ! A1.1 Astronomical constants
954 REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
955 ! A1.1.bis Constantes concernant l'orbite de la Terre:
956 REAL R_ecc, R_peri, R_incl
957 ! A1.2 Geoide
958 REAL RA,RG,R1SA
959 ! A1.3 Radiation
960 ! REAL RSIGMA,RI0
961 REAL RSIGMA
962 ! A1.4 Thermodynamic gas phase
963 REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
964 REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
965 REAL RKAPPA,RETV, eps_w
966 ! A1.5,6 Thermodynamic liquid,solid phases
967 REAL RCW,RCS
968 ! A1.7 Thermodynamic transition of phase
969 REAL RLVTT,RLSTT,RLMLT,RTT,RATM
970 ! A1.8 Curve of saturation
971 REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
972 REAL RALPD,RBETD,RGAMD
973 !
974 COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO &
975 & ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA &
976 & ,R_ecc, R_peri, R_incl &
977 & ,RA ,RG ,R1SA &
978 & ,RSIGMA &
979 & ,R ,RMD ,RMV ,RD ,RV ,RCPD &
980 & ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 &
981 & ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w &
982 & ,RCW ,RCS &
983 & ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM &
984 & ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS &
985 & ,RALPD ,RBETD ,RGAMD
986 ! ------------------------------------------------------------------
987 !$OMP THREADPRIVATE(/YOMCST/)
988 !
989 ! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $
990 !
991 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
992 ! veillez n'utiliser que des ! pour les commentaires
993 ! et bien positionner les & des lignes de continuation
994 ! (les placer en colonne 6 et en colonne 73)
995 !
996 !* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS
997 !
998 ! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION
999 ! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR
1000 ! ICE(*R_IES*).
1001 ! *RVTMP2* *RVTMP2=RCPV/RCPD-1.
1002 ! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.)
1003 !
1004 REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES
1005 REAL RVTMP2, RHOH2O
1006 REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
1007 REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
1008 LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
1009 ! If FALSE, then variables set by suphel.F90
1010 COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, &
1011 & RVTMP2, RHOH2O, &
1012 & R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, &
1013 & RALFDCP,RTWAT,RTBER,RTBERCU, &
1014 & RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
1015 & RKOOP2, &
1016 & OK_BAD_ECMWF_THERMO
1017
1018 !$OMP THREADPRIVATE(/YOETHF/)
1019 !
1020 ! $Header$
1021 !
1022 !
1023 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1024 ! veillez n'utiliser que des ! pour les commentaires
1025 ! et bien positionner les & des lignes de continuation
1026 ! (les placer en colonne 6 et en colonne 73)
1027 !
1028 ! ------------------------------------------------------------------
1029 ! This COMDECK includes the Thermodynamical functions for the cy39
1030 ! ECMWF Physics package.
1031 ! Consistent with YOMCST Basic physics constants, assuming the
1032 ! partial pressure of water vapour is given by a first order
1033 ! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants
1034 ! in YOETHF
1035 ! ------------------------------------------------------------------
1036 REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG
1037 REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl
1038 LOGICAL thermcep
1039 PARAMETER (thermcep=.TRUE.)
1040 !
1041 FOEEW ( PTARG,PDELARG ) = EXP ( &
1042 & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
1043 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
1044 !
1045 FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG &
1046 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
1047 !
1048 qsats(ptarg) = 100.0 * 0.622 * 10.0 &
1049 & ** (2.07023 - 0.00320991 * ptarg &
1050 & - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
1051 qsatl(ptarg) = 100.0 * 0.622 * 10.0 &
1052 & ** (23.8319 - 2948.964 / ptarg &
1053 & - 5.028 * LOG10(ptarg) &
1054 & - 29810.16 * EXP( - 0.0699382 * ptarg) &
1055 & + 25.21935 * EXP( - 2999.924 / ptarg))
1056 !
1057 dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
1058 & +2484.896*LOG(10.)/ptarg**2 &
1059 & -0.00320991*LOG(10.))
1060 dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
1061 & (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
1062 & +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
1063 & +29810.16*0.0699382*EXP(-0.0699382*ptarg))
1064 integer :: iflag_thermals,nsplit_thermals
1065
1066 !!! nrlmd le 10/04/2012
1067 integer :: iflag_trig_bl,iflag_clos_bl
1068 integer :: tau_trig_shallow,tau_trig_deep
1069 real :: s_trig
1070 !!! fin nrlmd le 10/04/2012
1071
1072 real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30.
1073 real :: alp_bl_k
1074 real :: tau_thermals,fact_thermals_ed_dz
1075 integer,parameter :: w2di_thermals=0
1076 integer :: isplit
1077
1078 integer :: iflag_coupl,iflag_clos,iflag_wake
1079 integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure
1080
1081 common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure
1082 common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz
1083 common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
1084 common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
1085
1086 !!! nrlmd le 10/04/2012
1087 common/ctherm6/iflag_trig_bl,iflag_clos_bl
1088 common/ctherm7/tau_trig_shallow,tau_trig_deep
1089 common/ctherm8/s_trig
1090 !!! fin nrlmd le 10/04/2012
1091
1092 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
1093 !$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
1094 !
1095 ! $Id: nuage.h 2945 2017-07-12 14:20:24Z jbmadeleine $
1096 !
1097 REAL rad_froid, rad_chau1, rad_chau2, t_glace_max, t_glace_min
1098 REAL exposant_glace
1099 REAL rei_min,rei_max
1100 REAL tau_cld_cv,coefw_cld_cv
1101
1102 REAL tmax_fonte_cv
1103
1104 INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
1105 INTEGER iflag_rain_incloud_vol
1106
1107 common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max, &
1108 & t_glace_min,exposant_glace,rei_min,rei_max, &
1109 & tau_cld_cv,coefw_cld_cv, &
1110 & tmax_fonte_cv, &
1111 & iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv, &
1112 & iflag_rain_incloud_vol
1113 !$OMP THREADPRIVATE(/nuagecom/)
1114
1115 INTEGER itap,ind1,ind2
1116 INTEGER ngrid,klev,klon,l,ig
1117
1118 REAL ztv(ngrid,klev)
1119 REAL po(ngrid)
1120 37440 REAL zqenv(ngrid)
1121 REAL zqta(ngrid,klev)
1122
1123 REAL fraca(ngrid,klev+1)
1124 REAL zpspsk(ngrid,klev)
1125 REAL paprs(ngrid,klev+1)
1126 REAL pplay(ngrid,klev)
1127 REAL ztla(ngrid,klev)
1128 REAL zthl(ngrid,klev)
1129
1130 37440 REAL zqsatth(ngrid,klev)
1131 37440 REAL zqsatenv(ngrid,klev)
1132
1133 37440 REAL sigma1(ngrid,klev)
1134 37440 REAL sigma2(ngrid,klev)
1135 37440 REAL qlth(ngrid,klev)
1136 37440 REAL qlenv(ngrid,klev)
1137 37440 REAL qltot(ngrid,klev)
1138 37440 REAL cth(ngrid,klev)
1139 37440 REAL cenv(ngrid,klev)
1140 REAL ctot(ngrid,klev)
1141 37440 REAL cth_vol(ngrid,klev)
1142 37440 REAL cenv_vol(ngrid,klev)
1143 REAL ctot_vol(ngrid,klev)
1144 37440 REAL rneb(ngrid,klev)
1145 REAL t(ngrid,klev)
1146 REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
1147 REAL rdd,cppd,Lv
1148 REAL alth,alenv,ath,aenv
1149 REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
1150 REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1151 REAL Tbef,zdelta,qsatbef,zcor
1152 REAL qlbef
1153 REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
1154 REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
1155 REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
1156 REAL zqs(ngrid), qcloud(ngrid)
1157 REAL erf
1158
1159
1160
1161 18720 IF (iflag_cloudth_vert.GE.1) THEN
1162 CALL cloudth_vert_v3(ngrid,klev,ind2, &
1163 & ztv,po,zqta,fraca, &
1164 & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1165 18720 & ratqs,zqs,t)
1166 18720 RETURN
1167 ENDIF
1168 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1169
1170
1171 !-------------------------------------------------------------------------------
1172 ! Initialisation des variables r?elles
1173 !-------------------------------------------------------------------------------
1174 ✗ sigma1(:,:)=0.
1175 ✗ sigma2(:,:)=0.
1176 ✗ qlth(:,:)=0.
1177 ✗ qlenv(:,:)=0.
1178 ✗ qltot(:,:)=0.
1179 ✗ rneb(:,:)=0.
1180 ✗ qcloud(:)=0.
1181 ✗ cth(:,:)=0.
1182 ✗ cenv(:,:)=0.
1183 ✗ ctot(:,:)=0.
1184 ✗ cth_vol(:,:)=0.
1185 ✗ cenv_vol(:,:)=0.
1186 ✗ ctot_vol(:,:)=0.
1187 qsatmmussig1=0.
1188 qsatmmussig2=0.
1189 rdd=287.04
1190 cppd=1005.7
1191 pi=3.14159
1192 Lv=2.5e6
1193 sqrt2pi=sqrt(2.*pi)
1194 sqrt2=sqrt(2.)
1195 sqrtpi=sqrt(pi)
1196
1197
1198 !-------------------------------------------------------------------------------
1199 ! Cloud fraction in the thermals and standard deviation of the PDFs
1200 !-------------------------------------------------------------------------------
1201 ✗ do ind1=1,ngrid
1202
1203 ✗ if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
1204
1205 ✗ zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
1206
1207 ✗ Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1208 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1209 ✗ qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1210 ✗ qsatbef=MIN(0.5,qsatbef)
1211 ✗ zcor=1./(1.-retv*qsatbef)
1212 ✗ qsatbef=qsatbef*zcor
1213 ✗ zqsatenv(ind1,ind2)=qsatbef
1214
1215
1216 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84
1217 ✗ aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
1218 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
1219
1220 !po = qt de l'environnement ET des thermique
1221 !zqenv = qt environnement
1222 !zqsatenv = qsat environnement
1223 !zthl = Tl environnement
1224
1225
1226 ✗ Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1227 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1228 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1229 ✗ qsatbef=MIN(0.5,qsatbef)
1230 ✗ zcor=1./(1.-retv*qsatbef)
1231 ✗ qsatbef=qsatbef*zcor
1232 ✗ zqsatth(ind1,ind2)=qsatbef
1233
1234 ✗ alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84
1235 ✗ ath=1./(1.+(alth*Lv/cppd)) !al, p84
1236 ✗ sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
1237
1238 !zqta = qt thermals
1239 !zqsatth = qsat thermals
1240 !ztla = Tl thermals
1241
1242 !------------------------------------------------------------------------------
1243 ! s standard deviations
1244 !------------------------------------------------------------------------------
1245
1246 ! tests
1247 ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1248 ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
1249 ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
1250 ! final option
1251 ✗ sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
1252 ✗ sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1253
1254 !------------------------------------------------------------------------------
1255 ! Condensed water and cloud cover
1256 !------------------------------------------------------------------------------
1257 ✗ xth=sth/(sqrt2*sigma2s)
1258 ✗ xenv=senv/(sqrt2*sigma1s)
1259 ✗ cth(ind1,ind2)=0.5*(1.+1.*erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
1260 ✗ cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
1261 ✗ ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1262 ✗ ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1263
1264 ✗ qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
1265 ✗ qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
1266 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1267
1268 ✗ if (ctot(ind1,ind2).lt.1.e-10) then
1269 ✗ ctot(ind1,ind2)=0.
1270 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
1271 else
1272 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1273 endif
1274
1275 else ! Environnement only, follow the if l.110
1276
1277 ✗ zqenv(ind1)=po(ind1)
1278 ✗ Tbef=t(ind1,ind2)
1279 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1280 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1281 ✗ qsatbef=MIN(0.5,qsatbef)
1282 ✗ zcor=1./(1.-retv*qsatbef)
1283 ✗ qsatbef=qsatbef*zcor
1284 ✗ zqsatenv(ind1,ind2)=qsatbef
1285
1286 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1287 ✗ zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1288 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
1289 ✗ aenv=1./(1.+(alenv*Lv/cppd))
1290 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1291
1292 ✗ sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1293
1294 ✗ xenv=senv/(sqrt2*sigma1s)
1295 ✗ ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1296 ✗ ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1297 ✗ qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
1298
1299 ✗ if (ctot(ind1,ind2).lt.1.e-3) then
1300 ✗ ctot(ind1,ind2)=0.
1301 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
1302 else
1303 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1304 endif
1305
1306
1307 endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
1308 enddo ! from the loop on ngrid l.108
1309 return
1310 ! end
1311 ✗ END SUBROUTINE cloudth_v3
1312
1313
1314
1315 !===========================================================================
1316 18720 SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, &
1317 & ztv,po,zqta,fraca, &
1318 18720 & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1319 & ratqs,zqs,t)
1320
1321 !===========================================================================
1322 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
1323 ! Date : 25 Mai 2010
1324 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
1325 !===========================================================================
1326
1327
1328
1/6
✓ Branch 0 taken 18720 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 18720 USE ioipsl_getin_p_mod, ONLY : getin_p
1329 USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
1330 & cloudth_sigmath,cloudth_sigmaenv
1331
1332 IMPLICIT NONE
1333
1334 !
1335 ! $Header$
1336 !
1337 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1338 ! veillez � n'utiliser que des ! pour les commentaires
1339 ! et � bien positionner les & des lignes de continuation
1340 ! (les placer en colonne 6 et en colonne 73)
1341 !
1342 !
1343 ! A1.0 Fundamental constants
1344 REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
1345 ! A1.1 Astronomical constants
1346 REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
1347 ! A1.1.bis Constantes concernant l'orbite de la Terre:
1348 REAL R_ecc, R_peri, R_incl
1349 ! A1.2 Geoide
1350 REAL RA,RG,R1SA
1351 ! A1.3 Radiation
1352 ! REAL RSIGMA,RI0
1353 REAL RSIGMA
1354 ! A1.4 Thermodynamic gas phase
1355 REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
1356 REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
1357 REAL RKAPPA,RETV, eps_w
1358 ! A1.5,6 Thermodynamic liquid,solid phases
1359 REAL RCW,RCS
1360 ! A1.7 Thermodynamic transition of phase
1361 REAL RLVTT,RLSTT,RLMLT,RTT,RATM
1362 ! A1.8 Curve of saturation
1363 REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
1364 REAL RALPD,RBETD,RGAMD
1365 !
1366 COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO &
1367 & ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA &
1368 & ,R_ecc, R_peri, R_incl &
1369 & ,RA ,RG ,R1SA &
1370 & ,RSIGMA &
1371 & ,R ,RMD ,RMV ,RD ,RV ,RCPD &
1372 & ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 &
1373 & ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w &
1374 & ,RCW ,RCS &
1375 & ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM &
1376 & ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS &
1377 & ,RALPD ,RBETD ,RGAMD
1378 ! ------------------------------------------------------------------
1379 !$OMP THREADPRIVATE(/YOMCST/)
1380 !
1381 ! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $
1382 !
1383 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1384 ! veillez n'utiliser que des ! pour les commentaires
1385 ! et bien positionner les & des lignes de continuation
1386 ! (les placer en colonne 6 et en colonne 73)
1387 !
1388 !* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS
1389 !
1390 ! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION
1391 ! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR
1392 ! ICE(*R_IES*).
1393 ! *RVTMP2* *RVTMP2=RCPV/RCPD-1.
1394 ! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.)
1395 !
1396 REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES
1397 REAL RVTMP2, RHOH2O
1398 REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
1399 REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
1400 LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
1401 ! If FALSE, then variables set by suphel.F90
1402 COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, &
1403 & RVTMP2, RHOH2O, &
1404 & R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, &
1405 & RALFDCP,RTWAT,RTBER,RTBERCU, &
1406 & RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
1407 & RKOOP2, &
1408 & OK_BAD_ECMWF_THERMO
1409
1410 !$OMP THREADPRIVATE(/YOETHF/)
1411 !
1412 ! $Header$
1413 !
1414 !
1415 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1416 ! veillez n'utiliser que des ! pour les commentaires
1417 ! et bien positionner les & des lignes de continuation
1418 ! (les placer en colonne 6 et en colonne 73)
1419 !
1420 ! ------------------------------------------------------------------
1421 ! This COMDECK includes the Thermodynamical functions for the cy39
1422 ! ECMWF Physics package.
1423 ! Consistent with YOMCST Basic physics constants, assuming the
1424 ! partial pressure of water vapour is given by a first order
1425 ! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants
1426 ! in YOETHF
1427 ! ------------------------------------------------------------------
1428 REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG
1429 REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl
1430 LOGICAL thermcep
1431 PARAMETER (thermcep=.TRUE.)
1432 !
1433 FOEEW ( PTARG,PDELARG ) = EXP ( &
1434 & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
1435 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
1436 !
1437 FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG &
1438 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
1439 !
1440 qsats(ptarg) = 100.0 * 0.622 * 10.0 &
1441 & ** (2.07023 - 0.00320991 * ptarg &
1442 & - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
1443 qsatl(ptarg) = 100.0 * 0.622 * 10.0 &
1444 & ** (23.8319 - 2948.964 / ptarg &
1445 & - 5.028 * LOG10(ptarg) &
1446 & - 29810.16 * EXP( - 0.0699382 * ptarg) &
1447 & + 25.21935 * EXP( - 2999.924 / ptarg))
1448 !
1449 dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
1450 & +2484.896*LOG(10.)/ptarg**2 &
1451 & -0.00320991*LOG(10.))
1452 dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
1453 & (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
1454 & +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
1455 & +29810.16*0.0699382*EXP(-0.0699382*ptarg))
1456 integer :: iflag_thermals,nsplit_thermals
1457
1458 !!! nrlmd le 10/04/2012
1459 integer :: iflag_trig_bl,iflag_clos_bl
1460 integer :: tau_trig_shallow,tau_trig_deep
1461 real :: s_trig
1462 !!! fin nrlmd le 10/04/2012
1463
1464 real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30.
1465 real :: alp_bl_k
1466 real :: tau_thermals,fact_thermals_ed_dz
1467 integer,parameter :: w2di_thermals=0
1468 integer :: isplit
1469
1470 integer :: iflag_coupl,iflag_clos,iflag_wake
1471 integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure
1472
1473 common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure
1474 common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz
1475 common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
1476 common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
1477
1478 !!! nrlmd le 10/04/2012
1479 common/ctherm6/iflag_trig_bl,iflag_clos_bl
1480 common/ctherm7/tau_trig_shallow,tau_trig_deep
1481 common/ctherm8/s_trig
1482 !!! fin nrlmd le 10/04/2012
1483
1484 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
1485 !$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
1486 !
1487 ! $Id: nuage.h 2945 2017-07-12 14:20:24Z jbmadeleine $
1488 !
1489 REAL rad_froid, rad_chau1, rad_chau2, t_glace_max, t_glace_min
1490 REAL exposant_glace
1491 REAL rei_min,rei_max
1492 REAL tau_cld_cv,coefw_cld_cv
1493
1494 REAL tmax_fonte_cv
1495
1496 INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
1497 INTEGER iflag_rain_incloud_vol
1498
1499 common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max, &
1500 & t_glace_min,exposant_glace,rei_min,rei_max, &
1501 & tau_cld_cv,coefw_cld_cv, &
1502 & tmax_fonte_cv, &
1503 & iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv, &
1504 & iflag_rain_incloud_vol
1505 !$OMP THREADPRIVATE(/nuagecom/)
1506
1507 INTEGER itap,ind1,ind2
1508 INTEGER ngrid,klev,klon,l,ig
1509
1510 REAL ztv(ngrid,klev)
1511 REAL po(ngrid)
1512 37440 REAL zqenv(ngrid)
1513 REAL zqta(ngrid,klev)
1514
1515 REAL fraca(ngrid,klev+1)
1516 REAL zpspsk(ngrid,klev)
1517 REAL paprs(ngrid,klev+1)
1518 REAL pplay(ngrid,klev)
1519 REAL ztla(ngrid,klev)
1520 REAL zthl(ngrid,klev)
1521
1522 37440 REAL zqsatth(ngrid,klev)
1523 37440 REAL zqsatenv(ngrid,klev)
1524
1525 37440 REAL sigma1(ngrid,klev)
1526 37440 REAL sigma2(ngrid,klev)
1527 37440 REAL qlth(ngrid,klev)
1528 37440 REAL qlenv(ngrid,klev)
1529 37440 REAL qltot(ngrid,klev)
1530 37440 REAL cth(ngrid,klev)
1531 37440 REAL cenv(ngrid,klev)
1532 REAL ctot(ngrid,klev)
1533 37440 REAL cth_vol(ngrid,klev)
1534 37440 REAL cenv_vol(ngrid,klev)
1535 REAL ctot_vol(ngrid,klev)
1536 37440 REAL rneb(ngrid,klev)
1537 REAL t(ngrid,klev)
1538 REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
1539 REAL rdd,cppd,Lv
1540 REAL alth,alenv,ath,aenv
1541 REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
1542 REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
1543 REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
1544 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
1545 REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
1546 REAL Tbef,zdelta,qsatbef,zcor
1547 REAL qlbef
1548 REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
1549 ! Change the width of the PDF used for vertical subgrid scale heterogeneity
1550 ! (J Jouhaud, JL Dufresne, JB Madeleine)
1551 REAL,SAVE :: vert_alpha, vert_alpha_th
1552 !$OMP THREADPRIVATE(vert_alpha, vert_alpha_th)
1553 REAL,SAVE :: sigma1s_factor=1.1
1554 REAL,SAVE :: sigma1s_power=0.6
1555 REAL,SAVE :: sigma2s_factor=0.09
1556 REAL,SAVE :: sigma2s_power=0.5
1557 REAL,SAVE :: cloudth_ratqsmin=-1.
1558 !$OMP THREADPRIVATE(sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin)
1559 INTEGER, SAVE :: iflag_cloudth_vert_noratqs=0
1560 !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)
1561
1562 LOGICAL, SAVE :: firstcall = .TRUE.
1563 !$OMP THREADPRIVATE(firstcall)
1564
1565 REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
1566 REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
1567 REAL zqs(ngrid), qcloud(ngrid)
1568 REAL erf
1569
1570 37440 REAL rhodz(ngrid,klev)
1571 37440 REAL zrho(ngrid,klev)
1572 18720 REAL dz(ngrid,klev)
1573
1574
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 DO ind1 = 1, ngrid
1575 !Layer calculation
1576 18607680 rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
1577 18607680 zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
1578 18720 dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
1579 END DO
1580
1581
1582 !------------------------------------------------------------------------------
1583 ! Initialize
1584 !------------------------------------------------------------------------------
1585
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 sigma1(:,:)=0.
1586
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 sigma2(:,:)=0.
1587
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 qlth(:,:)=0.
1588
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 qlenv(:,:)=0.
1589
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 qltot(:,:)=0.
1590
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 rneb(:,:)=0.
1591
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 qcloud(:)=0.
1592
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 cth(:,:)=0.
1593
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 cenv(:,:)=0.
1594
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 ctot(:,:)=0.
1595
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 cth_vol(:,:)=0.
1596
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 cenv_vol(:,:)=0.
1597
4/4
✓ Branch 0 taken 730080 times.
✓ Branch 1 taken 18720 times.
✓ Branch 2 taken 725699520 times.
✓ Branch 3 taken 730080 times.
 726448320 ctot_vol(:,:)=0.
1598 qsatmmussig1=0.
1599 qsatmmussig2=0.
1600 rdd=287.04
1601 cppd=1005.7
1602 pi=3.14159
1603 Lv=2.5e6
1604 sqrt2pi=sqrt(2.*pi)
1605 sqrt2=sqrt(2.)
1606 sqrtpi=sqrt(pi)
1607
1608
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 18719 times.
 18720 IF (firstcall) THEN
1609 1 vert_alpha=0.5
1610 1 CALL getin_p('cloudth_vert_alpha',vert_alpha)
1611 1 WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
1612 ! The factor used for the thermal is equal to that of the environment
1613 ! if nothing is explicitly specified in the def file
1614 1 vert_alpha_th=vert_alpha
1615 1 CALL getin_p('cloudth_vert_alpha_th',vert_alpha_th)
1616 1 WRITE(*,*) 'cloudth_vert_alpha_th = ', vert_alpha_th
1617 ! Factor used in the calculation of sigma1s
1618 1 CALL getin_p('cloudth_sigma1s_factor',sigma1s_factor)
1619 1 WRITE(*,*) 'cloudth_sigma1s_factor = ', sigma1s_factor
1620 ! Power used in the calculation of sigma1s
1621 1 CALL getin_p('cloudth_sigma1s_power',sigma1s_power)
1622 1 WRITE(*,*) 'cloudth_sigma1s_power = ', sigma1s_power
1623 ! Factor used in the calculation of sigma2s
1624 1 CALL getin_p('cloudth_sigma2s_factor',sigma2s_factor)
1625 1 WRITE(*,*) 'cloudth_sigma2s_factor = ', sigma2s_factor
1626 ! Power used in the calculation of sigma2s
1627 1 CALL getin_p('cloudth_sigma2s_power',sigma2s_power)
1628 1 WRITE(*,*) 'cloudth_sigma2s_power = ', sigma2s_power
1629 ! Minimum value for the environmental air subgrid water distrib
1630 1 CALL getin_p('cloudth_ratqsmin',cloudth_ratqsmin)
1631 1 WRITE(*,*) 'cloudth_ratqsmin = ', cloudth_ratqsmin
1632 ! Remove the dependency to ratqs from the variance of the vertical PDF
1633 1 CALL getin_p('iflag_cloudth_vert_noratqs',iflag_cloudth_vert_noratqs)
1634 1 WRITE(*,*) 'iflag_cloudth_vert_noratqs = ', iflag_cloudth_vert_noratqs
1635
1636 1 firstcall=.FALSE.
1637 ENDIF
1638
1639 !-------------------------------------------------------------------------------
1640 ! Calcul de la fraction du thermique et des ecart-types des distributions
1641 !-------------------------------------------------------------------------------
1642
2/2
✓ Branch 0 taken 18607680 times.
✓ Branch 1 taken 18720 times.
 18626400 do ind1=1,ngrid
1643
1644
4/4
✓ Branch 0 taken 9678123 times.
✓ Branch 1 taken 8929557 times.
✓ Branch 2 taken 1304630 times.
✓ Branch 3 taken 8373493 times.
 18607680 if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
1645
1646 1304630 zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
1647
1648
1649 1304630 Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
1650 1304630 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1651 1304630 qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1652 1304630 qsatbef=MIN(0.5,qsatbef)
1653 1304630 zcor=1./(1.-retv*qsatbef)
1654 1304630 qsatbef=qsatbef*zcor
1655 1304630 zqsatenv(ind1,ind2)=qsatbef
1656
1657
1658 1304630 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84
1659 1304630 aenv=1./(1.+(alenv*Lv/cppd)) !al, p84
1660 1304630 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84
1661
1662 !zqenv = qt environnement
1663 !zqsatenv = qsat environnement
1664 !zthl = Tl environnement
1665
1666
1667 1304630 Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
1668 1304630 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1669 1304630 qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1670 1304630 qsatbef=MIN(0.5,qsatbef)
1671 1304630 zcor=1./(1.-retv*qsatbef)
1672 1304630 qsatbef=qsatbef*zcor
1673 1304630 zqsatth(ind1,ind2)=qsatbef
1674
1675 1304630 alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84
1676 1304630 ath=1./(1.+(alth*Lv/cppd)) !al, p84
1677 1304630 sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84
1678
1679
1680 !zqta = qt thermals
1681 !zqsatth = qsat thermals
1682 !ztla = Tl thermals
1683
1684 !------------------------------------------------------------------------------
1685 ! s standard deviation
1686 !------------------------------------------------------------------------------
1687
1688 sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
1689 1304630 & (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
1690 ! sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
1691 1304630 IF (cloudth_ratqsmin>0.) THEN
1692 ✗ sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
1693 ELSE
1694 1304630 sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
1695 ENDIF
1696 1304630 sigma1s = sigma1s_fraca + sigma1s_ratqs
1697 1304630 sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
1698 ! tests
1699 ! sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
1700 ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
1701 ! sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
1702 ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
1703 ! if (paprs(ind1,ind2).gt.90000) then
1704 ! ratqs(ind1,ind2)=0.002
1705 ! else
1706 ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
1707 ! endif
1708 ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1709 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1710 ! sigma1s=ratqs(ind1,ind2)*po(ind1)
1711 ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003
1712
1713 1304630 IF (iflag_cloudth_vert == 1) THEN
1714 !-------------------------------------------------------------------------------
1715 ! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
1716 !-------------------------------------------------------------------------------
1717
1718 ✗ deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1719 ✗ deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1720
1721 ✗ xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
1722 ✗ xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
1723 ✗ xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
1724 ✗ xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
1725 ✗ coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
1726 ✗ coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
1727
1728 ✗ cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
1729 ✗ cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2))
1730 ✗ ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1731
1732 ! Environment
1733 ✗ IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
1734 ✗ IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
1735 ✗ IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
1736 ✗ IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
1737
1738 ✗ qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1739
1740 ! Thermal
1741 ✗ IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
1742 ✗ IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
1743 ✗ IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
1744 ✗ IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
1745 ✗ qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1746 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1747
1748
1/2
✓ Branch 0 taken 1304630 times.
✗ Branch 1 not taken.
 1304630 ELSE IF (iflag_cloudth_vert >= 3) THEN
1749
1/2
✓ Branch 0 taken 1304630 times.
✗ Branch 1 not taken.
 1304630 IF (iflag_cloudth_vert < 5) THEN
1750 !-------------------------------------------------------------------------------
1751 ! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
1752 !-------------------------------------------------------------------------------
1753 ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
1754 ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
1755 ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
1756 ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
1757
1/2
✓ Branch 0 taken 1304630 times.
✗ Branch 1 not taken.
 1304630 IF (iflag_cloudth_vert == 3) THEN
1758 1304630 deltasenv=aenv*vert_alpha*sigma1s
1759 1304630 deltasth=ath*vert_alpha_th*sigma2s
1760 ELSE IF (iflag_cloudth_vert == 4) THEN
1761 ✗ IF (iflag_cloudth_vert_noratqs == 1) THEN
1762 ✗ deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
1763 ✗ deltasth=vert_alpha_th*sigma2s
1764 ELSE
1765 ✗ deltasenv=vert_alpha*sigma1s
1766 ✗ deltasth=vert_alpha_th*sigma2s
1767 ENDIF
1768 ENDIF
1769
1770 1304630 xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
1771 1304630 xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
1772 1304630 exp_xenv1 = exp(-1.*xenv1**2)
1773 1304630 exp_xenv2 = exp(-1.*xenv2**2)
1774 1304630 xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
1775 1304630 xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
1776 1304630 exp_xth1 = exp(-1.*xth1**2)
1777 1304630 exp_xth2 = exp(-1.*xth2**2)
1778
1779 !CF_surfacique
1780 1304630 cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
1781 1304630 cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
1782 1304630 ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1783
1784
1785 !CF_volumique & eau condense
1786 !environnement
1787 1304630 IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
1788 1304630 IntJ_CF=0.5*(1.-1.*erf(xenv2))
1789
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1304630 times.
 1304630 if (deltasenv .lt. 1.e-10) then
1790 ✗ qlenv(ind1,ind2)=IntJ
1791 ✗ cenv_vol(ind1,ind2)=IntJ_CF
1792 else
1793 1304630 IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
1794 1304630 IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
1795 1304630 IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
1796 1304630 IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
1797 1304630 IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
1798 1304630 qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1799 1304630 cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1800 endif
1801
1802 !thermique
1803 1304630 IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
1804 1304630 IntJ_CF=0.5*(1.-1.*erf(xth2))
1805
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1304630 times.
 1304630 if (deltasth .lt. 1.e-10) then
1806 ✗ qlth(ind1,ind2)=IntJ
1807 ✗ cth_vol(ind1,ind2)=IntJ_CF
1808 else
1809 1304630 IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
1810 1304630 IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
1811 1304630 IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
1812 1304630 IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
1813 1304630 IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
1814 1304630 qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
1815 1304630 cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
1816 endif
1817
1818 1304630 qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1819 1304630 ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1820
1821 ✗ ELSE IF (iflag_cloudth_vert == 5) THEN
1822 ✗ sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1) !Environment
1823 ✗ sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2) !Thermals
1824 !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
1825 !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
1826 ✗ xth=sth/(sqrt(2.)*sigma2s)
1827 ✗ xenv=senv/(sqrt(2.)*sigma1s)
1828
1829 !Volumique
1830 ✗ cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
1831 ✗ cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1832 ✗ ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
1833 !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)
1834
1835 ✗ qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
1836 ✗ qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2))
1837 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
1838
1839 !Surfacique
1840 !Neggers
1841 !beta=0.0044
1842 !inverse_rho=1.+beta*dz(ind1,ind2)
1843 !print *,'jeanjean : beta=',beta
1844 !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
1845 !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
1846 !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
1847
1848 !Brooks
1849 a_Brooks=0.6694
1850 b_Brooks=0.1882
1851 A_Maj_Brooks=0.1635 !-- sans shear
1852 !A_Maj_Brooks=0.17 !-- ARM LES
1853 !A_Maj_Brooks=0.18 !-- RICO LES
1854 !A_Maj_Brooks=0.19 !-- BOMEX LES
1855 Dx_Brooks=200000.
1856 ✗ f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
1857 !print *,'jeanjean_f=',f_Brooks
1858
1859 ✗ cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
1860 ✗ cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
1861 ✗ ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
1862 !print *,'JJ_ctot_1',ctot(ind1,ind2)
1863
1864
1865
1866
1867
1868 ENDIF ! of if (iflag_cloudth_vert<5)
1869 ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
1870
1871 ! if (ctot(ind1,ind2).lt.1.e-10) then
1872
4/4
✓ Branch 0 taken 629590 times.
✓ Branch 1 taken 675040 times.
✓ Branch 2 taken 62464 times.
✓ Branch 3 taken 567126 times.
 1304630 if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
1873 737504 ctot(ind1,ind2)=0.
1874 737504 ctot_vol(ind1,ind2)=0.
1875 737504 qcloud(ind1)=zqsatenv(ind1,ind2)
1876
1877 else
1878
1879 567126 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
1880 ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
1881 ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
1882
1883 endif
1884
1885 else ! gaussienne environnement seule
1886
1887 17303050 zqenv(ind1)=po(ind1)
1888 17303050 Tbef=t(ind1,ind2)
1889 17303050 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
1890 17303050 qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
1891 17303050 qsatbef=MIN(0.5,qsatbef)
1892 17303050 zcor=1./(1.-retv*qsatbef)
1893 17303050 qsatbef=qsatbef*zcor
1894 17303050 zqsatenv(ind1,ind2)=qsatbef
1895
1896
1897 ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
1898 17303050 zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
1899 17303050 alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
1900 17303050 aenv=1./(1.+(alenv*Lv/cppd))
1901 17303050 senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
1902 sth=0.
1903
1904
1905 17303050 sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
1906 sigma2s=0.
1907
1908 sqrt2pi=sqrt(2.*pi)
1909 17303050 xenv=senv/(sqrt(2.)*sigma1s)
1910 17303050 ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
1911 17303050 ctot_vol(ind1,ind2)=ctot(ind1,ind2)
1912 17303050 qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
1913
1914 17303050 if (ctot(ind1,ind2).lt.1.e-3) then
1915 13770146 ctot(ind1,ind2)=0.
1916 13770146 qcloud(ind1)=zqsatenv(ind1,ind2)
1917
1918 else
1919
1920 ! ctot(ind1,ind2)=ctot(ind1,ind2)
1921 3532904 qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
1922
1923 endif
1924
1925
1926
1927
1928 endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
1929 ! Outputs used to check the PDFs
1930 18607680 cloudth_senv(ind1,ind2) = senv
1931 18607680 cloudth_sth(ind1,ind2) = sth
1932 18607680 cloudth_sigmaenv(ind1,ind2) = sigma1s
1933 18626400 cloudth_sigmath(ind1,ind2) = sigma2s
1934
1935 enddo ! from the loop on ngrid l.333
1936 18720 return
1937 ! end
1938 END SUBROUTINE cloudth_vert_v3
1939 !
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951 ✗ SUBROUTINE cloudth_v6(ngrid,klev,ind2, &
1952 & ztv,po,zqta,fraca, &
1953 ✗ & qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
1954 & ratqs,zqs,T)
1955
1956
1957
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1304630 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1304630 times.
✓ Branch 4 taken 13770146 times.
✓ Branch 5 taken 3532904 times.
 38519990 USE ioipsl_getin_p_mod, ONLY : getin_p
1958 USE phys_output_var_mod, ONLY : cloudth_sth,cloudth_senv, &
1959 & cloudth_sigmath,cloudth_sigmaenv
1960
1961 IMPLICIT NONE
1962
1963 !
1964 ! $Header$
1965 !
1966 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1967 ! veillez � n'utiliser que des ! pour les commentaires
1968 ! et � bien positionner les & des lignes de continuation
1969 ! (les placer en colonne 6 et en colonne 73)
1970 !
1971 !
1972 ! A1.0 Fundamental constants
1973 REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
1974 ! A1.1 Astronomical constants
1975 REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
1976 ! A1.1.bis Constantes concernant l'orbite de la Terre:
1977 REAL R_ecc, R_peri, R_incl
1978 ! A1.2 Geoide
1979 REAL RA,RG,R1SA
1980 ! A1.3 Radiation
1981 ! REAL RSIGMA,RI0
1982 REAL RSIGMA
1983 ! A1.4 Thermodynamic gas phase
1984 REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
1985 REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
1986 REAL RKAPPA,RETV, eps_w
1987 ! A1.5,6 Thermodynamic liquid,solid phases
1988 REAL RCW,RCS
1989 ! A1.7 Thermodynamic transition of phase
1990 REAL RLVTT,RLSTT,RLMLT,RTT,RATM
1991 ! A1.8 Curve of saturation
1992 REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
1993 REAL RALPD,RBETD,RGAMD
1994 !
1995 COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO &
1996 & ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA &
1997 & ,R_ecc, R_peri, R_incl &
1998 & ,RA ,RG ,R1SA &
1999 & ,RSIGMA &
2000 & ,R ,RMD ,RMV ,RD ,RV ,RCPD &
2001 & ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 &
2002 & ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w &
2003 & ,RCW ,RCS &
2004 & ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM &
2005 & ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS &
2006 & ,RALPD ,RBETD ,RGAMD
2007 ! ------------------------------------------------------------------
2008 !$OMP THREADPRIVATE(/YOMCST/)
2009 !
2010 ! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $
2011 !
2012 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
2013 ! veillez n'utiliser que des ! pour les commentaires
2014 ! et bien positionner les & des lignes de continuation
2015 ! (les placer en colonne 6 et en colonne 73)
2016 !
2017 !* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS
2018 !
2019 ! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION
2020 ! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR
2021 ! ICE(*R_IES*).
2022 ! *RVTMP2* *RVTMP2=RCPV/RCPD-1.
2023 ! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.)
2024 !
2025 REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES
2026 REAL RVTMP2, RHOH2O
2027 REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU
2028 REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2
2029 LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90
2030 ! If FALSE, then variables set by suphel.F90
2031 COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, &
2032 & RVTMP2, RHOH2O, &
2033 & R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, &
2034 & RALFDCP,RTWAT,RTBER,RTBERCU, &
2035 & RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
2036 & RKOOP2, &
2037 & OK_BAD_ECMWF_THERMO
2038
2039 !$OMP THREADPRIVATE(/YOETHF/)
2040 !
2041 ! $Header$
2042 !
2043 !
2044 ! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
2045 ! veillez n'utiliser que des ! pour les commentaires
2046 ! et bien positionner les & des lignes de continuation
2047 ! (les placer en colonne 6 et en colonne 73)
2048 !
2049 ! ------------------------------------------------------------------
2050 ! This COMDECK includes the Thermodynamical functions for the cy39
2051 ! ECMWF Physics package.
2052 ! Consistent with YOMCST Basic physics constants, assuming the
2053 ! partial pressure of water vapour is given by a first order
2054 ! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants
2055 ! in YOETHF
2056 ! ------------------------------------------------------------------
2057 REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG
2058 REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl
2059 LOGICAL thermcep
2060 PARAMETER (thermcep=.TRUE.)
2061 !
2062 FOEEW ( PTARG,PDELARG ) = EXP ( &
2063 & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
2064 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
2065 !
2066 FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG &
2067 & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2
2068 !
2069 qsats(ptarg) = 100.0 * 0.622 * 10.0 &
2070 & ** (2.07023 - 0.00320991 * ptarg &
2071 & - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
2072 qsatl(ptarg) = 100.0 * 0.622 * 10.0 &
2073 & ** (23.8319 - 2948.964 / ptarg &
2074 & - 5.028 * LOG10(ptarg) &
2075 & - 29810.16 * EXP( - 0.0699382 * ptarg) &
2076 & + 25.21935 * EXP( - 2999.924 / ptarg))
2077 !
2078 dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
2079 & +2484.896*LOG(10.)/ptarg**2 &
2080 & -0.00320991*LOG(10.))
2081 dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
2082 & (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
2083 & +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
2084 & +29810.16*0.0699382*EXP(-0.0699382*ptarg))
2085 integer :: iflag_thermals,nsplit_thermals
2086
2087 !!! nrlmd le 10/04/2012
2088 integer :: iflag_trig_bl,iflag_clos_bl
2089 integer :: tau_trig_shallow,tau_trig_deep
2090 real :: s_trig
2091 !!! fin nrlmd le 10/04/2012
2092
2093 real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30.
2094 real :: alp_bl_k
2095 real :: tau_thermals,fact_thermals_ed_dz
2096 integer,parameter :: w2di_thermals=0
2097 integer :: isplit
2098
2099 integer :: iflag_coupl,iflag_clos,iflag_wake
2100 integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure
2101
2102 common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure
2103 common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz
2104 common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
2105 common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
2106
2107 !!! nrlmd le 10/04/2012
2108 common/ctherm6/iflag_trig_bl,iflag_clos_bl
2109 common/ctherm7/tau_trig_shallow,tau_trig_deep
2110 common/ctherm8/s_trig
2111 !!! fin nrlmd le 10/04/2012
2112
2113 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/)
2114 !$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/)
2115 !
2116 ! $Id: nuage.h 2945 2017-07-12 14:20:24Z jbmadeleine $
2117 !
2118 REAL rad_froid, rad_chau1, rad_chau2, t_glace_max, t_glace_min
2119 REAL exposant_glace
2120 REAL rei_min,rei_max
2121 REAL tau_cld_cv,coefw_cld_cv
2122
2123 REAL tmax_fonte_cv
2124
2125 INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv
2126 INTEGER iflag_rain_incloud_vol
2127
2128 common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max, &
2129 & t_glace_min,exposant_glace,rei_min,rei_max, &
2130 & tau_cld_cv,coefw_cld_cv, &
2131 & tmax_fonte_cv, &
2132 & iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv, &
2133 & iflag_rain_incloud_vol
2134 !$OMP THREADPRIVATE(/nuagecom/)
2135
2136
2137 !Domain variables
2138 INTEGER ngrid !indice Max lat-lon
2139 INTEGER klev !indice Max alt
2140 INTEGER ind1 !indice in [1:ngrid]
2141 INTEGER ind2 !indice in [1:klev]
2142 !thermal plume fraction
2143 REAL fraca(ngrid,klev+1) !thermal plumes fraction in the gridbox
2144 !temperatures
2145 REAL T(ngrid,klev) !temperature
2146 REAL zpspsk(ngrid,klev) !factor (p/p0)**kappa (used for potential variables)
2147 REAL ztv(ngrid,klev) !potential temperature (voir thermcell_env.F90)
2148 REAL ztla(ngrid,klev) !liquid temperature in the thermals (Tl_th)
2149 REAL zthl(ngrid,klev) !liquid temperature in the environment (Tl_env)
2150 !pressure
2151 REAL paprs(ngrid,klev+1) !pressure at the interface of levels
2152 REAL pplay(ngrid,klev) !pressure at the middle of the level
2153 !humidity
2154 REAL ratqs(ngrid,klev) !width of the total water subgrid-scale distribution
2155 REAL po(ngrid) !total water (qt)
2156 ✗ REAL zqenv(ngrid) !total water in the environment (qt_env)
2157 REAL zqta(ngrid,klev) !total water in the thermals (qt_th)
2158 ✗ REAL zqsatth(ngrid,klev) !water saturation level in the thermals (q_sat_th)
2159 ✗ REAL zqsatenv(ngrid,klev) !water saturation level in the environment (q_sat_env)
2160 ✗ REAL qlth(ngrid,klev) !condensed water in the thermals
2161 ✗ REAL qlenv(ngrid,klev) !condensed water in the environment
2162 ✗ REAL qltot(ngrid,klev) !condensed water in the gridbox
2163 !cloud fractions
2164 ✗ REAL cth_vol(ngrid,klev) !cloud fraction by volume in the thermals
2165 ✗ REAL cenv_vol(ngrid,klev) !cloud fraction by volume in the environment
2166 REAL ctot_vol(ngrid,klev) !cloud fraction by volume in the gridbox
2167 ✗ REAL cth_surf(ngrid,klev) !cloud fraction by surface in the thermals
2168 ✗ REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment
2169 REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
2170 !PDF of saturation deficit variables
2171 REAL rdd,cppd,Lv
2172 REAL Tbef,zdelta,qsatbef,zcor
2173 REAL alth,alenv,ath,aenv
2174 REAL sth,senv !saturation deficits in the thermals and environment
2175 REAL sigma_env,sigma_th !standard deviations of the biGaussian PDF
2176 !cloud fraction variables
2177 REAL xth,xenv
2178 REAL inverse_rho,beta !Neggers et al. (2011) method
2179 REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
2180 !Incloud total water variables
2181 REAL zqs(ngrid) !q_sat
2182 REAL qcloud(ngrid) !eau totale dans le nuage
2183 !Some arithmetic variables
2184 REAL erf,pi,sqrt2,sqrt2pi
2185 !Depth of the layer
2186 ✗ REAL dz(ngrid,klev) !epaisseur de la couche en metre
2187 ✗ REAL rhodz(ngrid,klev)
2188 ✗ REAL zrho(ngrid,klev)
2189 ✗ DO ind1 = 1, ngrid
2190 ✗ rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
2191 ✗ zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd ![kg/m3]
2192 ✗ dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) ![m]
2193 END DO
2194
2195 !------------------------------------------------------------------------------
2196 ! Initialization
2197 !------------------------------------------------------------------------------
2198 ✗ qlth(:,:)=0.
2199 ✗ qlenv(:,:)=0.
2200 ✗ qltot(:,:)=0.
2201 ✗ cth_vol(:,:)=0.
2202 ✗ cenv_vol(:,:)=0.
2203 ✗ ctot_vol(:,:)=0.
2204 ✗ cth_surf(:,:)=0.
2205 ✗ cenv_surf(:,:)=0.
2206 ✗ ctot_surf(:,:)=0.
2207 ✗ qcloud(:)=0.
2208 rdd=287.04
2209 cppd=1005.7
2210 pi=3.14159
2211 Lv=2.5e6
2212 sqrt2=sqrt(2.)
2213 sqrt2pi=sqrt(2.*pi)
2214
2215
2216 ✗ DO ind1=1,ngrid
2217 !-------------------------------------------------------------------------------
2218 !Both thermal and environment in the gridbox
2219 !-------------------------------------------------------------------------------
2220 ✗ IF ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) THEN
2221 !--------------------------------------------
2222 !calcul de qsat_env
2223 !--------------------------------------------
2224 ✗ Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
2225 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2226 ✗ qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2227 ✗ qsatbef=MIN(0.5,qsatbef)
2228 ✗ zcor=1./(1.-retv*qsatbef)
2229 ✗ qsatbef=qsatbef*zcor
2230 ✗ zqsatenv(ind1,ind2)=qsatbef
2231 !--------------------------------------------
2232 !calcul de s_env
2233 !--------------------------------------------
2234 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam
2235 ✗ aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
2236 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
2237 !--------------------------------------------
2238 !calcul de qsat_th
2239 !--------------------------------------------
2240 ✗ Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
2241 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2242 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2243 ✗ qsatbef=MIN(0.5,qsatbef)
2244 ✗ zcor=1./(1.-retv*qsatbef)
2245 ✗ qsatbef=qsatbef*zcor
2246 ✗ zqsatth(ind1,ind2)=qsatbef
2247 !--------------------------------------------
2248 !calcul de s_th
2249 !--------------------------------------------
2250 ✗ alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2) !qsl, p84 these Arnaud Jam
2251 ✗ ath=1./(1.+(alth*Lv/cppd)) !al, p84 these Arnaud Jam
2252 ✗ sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) !s, p84 these Arnaud Jam
2253 !--------------------------------------------
2254 !calcul standard deviations bi-Gaussian PDF
2255 !--------------------------------------------
2256 ✗ sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)
2257 ✗ sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5)+ratqs(ind1,ind2)*po(ind1)
2258 ✗ xth=sth/(sqrt2*sigma_th)
2259 ✗ xenv=senv/(sqrt2*sigma_env)
2260 !--------------------------------------------
2261 !Cloud fraction by volume CF_vol
2262 !--------------------------------------------
2263 ✗ cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
2264 ✗ cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
2265 ✗ ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
2266 !--------------------------------------------
2267 !Condensed water qc
2268 !--------------------------------------------
2269 ✗ qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
2270 ✗ qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2))
2271 ✗ qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
2272 !--------------------------------------------
2273 !Cloud fraction by surface CF_surf
2274 !--------------------------------------------
2275 !Method Neggers et al. (2011) : ok for cumulus clouds only
2276 !beta=0.0044 (Jouhaud et al.2018)
2277 !inverse_rho=1.+beta*dz(ind1,ind2)
2278 !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
2279 !Method Brooks et al. (2005) : ok for all types of clouds
2280 a_Brooks=0.6694
2281 b_Brooks=0.1882
2282 A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
2283 Dx_Brooks=200000. !-- si l'on considere des mailles de 200km de cote
2284 ✗ f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
2285 ✗ ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
2286 !--------------------------------------------
2287 !Incloud Condensed water qcloud
2288 !--------------------------------------------
2289 ✗ if (ctot_surf(ind1,ind2) .lt. 1.e-10) then
2290 ✗ ctot_vol(ind1,ind2)=0.
2291 ✗ ctot_surf(ind1,ind2)=0.
2292 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
2293 else
2294 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
2295 endif
2296
2297
2298
2299 !-------------------------------------------------------------------------------
2300 !Environment only in the gridbox
2301 !-------------------------------------------------------------------------------
2302 ELSE
2303 !--------------------------------------------
2304 !calcul de qsat_env
2305 !--------------------------------------------
2306 ✗ Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
2307 ✗ zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
2308 ✗ qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
2309 ✗ qsatbef=MIN(0.5,qsatbef)
2310 ✗ zcor=1./(1.-retv*qsatbef)
2311 ✗ qsatbef=qsatbef*zcor
2312 ✗ zqsatenv(ind1,ind2)=qsatbef
2313 !--------------------------------------------
2314 !calcul de s_env
2315 !--------------------------------------------
2316 ✗ alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) !qsl, p84 these Arnaud Jam
2317 ✗ aenv=1./(1.+(alenv*Lv/cppd)) !al, p84 these Arnaud Jam
2318 ✗ senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) !s, p84 these Arnaud Jam
2319 !--------------------------------------------
2320 !calcul standard deviations Gaussian PDF
2321 !--------------------------------------------
2322 ✗ zqenv(ind1)=po(ind1)
2323 ✗ sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
2324 ✗ xenv=senv/(sqrt2*sigma_env)
2325 !--------------------------------------------
2326 !Cloud fraction by volume CF_vol
2327 !--------------------------------------------
2328 ✗ ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
2329 !--------------------------------------------
2330 !Condensed water qc
2331 !--------------------------------------------
2332 ✗ qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
2333 !--------------------------------------------
2334 !Cloud fraction by surface CF_surf
2335 !--------------------------------------------
2336 !Method Neggers et al. (2011) : ok for cumulus clouds only
2337 !beta=0.0044 (Jouhaud et al.2018)
2338 !inverse_rho=1.+beta*dz(ind1,ind2)
2339 !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
2340 !Method Brooks et al. (2005) : ok for all types of clouds
2341 a_Brooks=0.6694
2342 b_Brooks=0.1882
2343 A_Maj_Brooks=0.1635 !-- sans dependence au shear
2344 Dx_Brooks=200000.
2345 ✗ f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
2346 ✗ ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
2347 !--------------------------------------------
2348 !Incloud Condensed water qcloud
2349 !--------------------------------------------
2350 ✗ if (ctot_surf(ind1,ind2) .lt. 1.e-8) then
2351 ✗ ctot_vol(ind1,ind2)=0.
2352 ✗ ctot_surf(ind1,ind2)=0.
2353 ✗ qcloud(ind1)=zqsatenv(ind1,ind2)
2354 else
2355 ✗ qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
2356 endif
2357
2358
2359 END IF ! From the separation (thermal/envrionnement) et (environnement only)
2360
2361 ! Outputs used to check the PDFs
2362 ✗ cloudth_senv(ind1,ind2) = senv
2363 ✗ cloudth_sth(ind1,ind2) = sth
2364 ✗ cloudth_sigmaenv(ind1,ind2) = sigma_env
2365 ✗ cloudth_sigmath(ind1,ind2) = sigma_th
2366
2367 END DO ! From the loop on ngrid
2368 ✗ return
2369
2370 END SUBROUTINE cloudth_v6
2371 END MODULE cloudth_mod
2372
2373
2374
2375
2376