GCC Code Coverage Report


Directory: ./
File: phys/nuage.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 110 0.0%
Branches: 0 82 0.0%

Line Branch Exec Source
1 ! $Id: nuage.F90 2346 2015-08-21 15:13:46Z emillour $
2
3 ✗ SUBROUTINE nuage(paprs, pplay, t, pqlwp, pclc, pcltau, pclemi, pch, pcl, pcm, &
4 ✗ pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, &
5 cldtaupi, re, fl)
6 USE dimphy
7 USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
8 IMPLICIT NONE
9 ! ======================================================================
10 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
11 ! Objet: Calculer epaisseur optique et emmissivite des nuages
12 ! ======================================================================
13 ! Arguments:
14 ! t-------input-R-temperature
15 ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
16 ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
17 ! ok_aie--input-L-apply aerosol indirect effect or not
18 ! mass_solu_aero-----input-R-total mass concentration for all soluble
19 ! aerosols[ug/m^3]
20 ! mass_solu_aero_pi--input-R-dito, pre-industrial value
21 ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
22 ! bl95_b1-input-R-a parameter, may be varied for tests ( -"- )
23
24 ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
25 ! needed for the diagnostics of the aerosol indirect
26 ! radiative forcing (see radlwsw)
27 ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
28 ! fl------output-R-Denominator to re, introduced to avoid problems in
29 ! the averaging of the output. fl is the fraction of liquid
30 ! water clouds within a grid cell
31
32 ! pcltau--output-R-epaisseur optique des nuages
33 ! pclemi--output-R-emissivite des nuages (0 a 1)
34 ! ======================================================================
35
36 include "YOMCST.h"
37 include "nuage.h" ! JBM 3/14
38
39 REAL paprs(klon, klev+1), pplay(klon, klev)
40 REAL t(klon, klev)
41
42 REAL pclc(klon, klev)
43 REAL pqlwp(klon, klev)
44 REAL pcltau(klon, klev), pclemi(klon, klev)
45
46 REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
47
48 LOGICAL lo
49
50 REAL cetahb, cetamb
51 PARAMETER (cetahb=0.45, cetamb=0.80)
52
53 INTEGER i, k
54 ✗ REAL zflwp, zradef, zfice(klon), zmsac
55
56 REAL radius, rad_chaud
57 ! JBM (3/14) parameters already defined in nuage.h:
58 ! REAL rad_froid, rad_chau1, rad_chau2
59 ! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
60 ! cc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
61 ! sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
62 REAL coef, coef_froi, coef_chau
63 PARAMETER (coef_chau=0.13, coef_froi=0.09)
64 REAL seuil_neb
65 PARAMETER (seuil_neb=0.001)
66 ! JBM (3/14) nexpo is replaced by exposant_glace
67 ! REAL nexpo ! exponentiel pour glace/eau
68 ! PARAMETER (nexpo=6.)
69 REAL, PARAMETER :: t_glace_min_old = 258.
70 INTEGER, PARAMETER :: exposant_glace_old = 6
71
72
73 ! jq for the aerosol indirect effect
74 ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
75 ! jq
76 LOGICAL ok_aie ! Apply AIE or not?
77
78 REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols[ug m-3]
79 REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
80 ✗ REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
81 REAL re(klon, klev) ! cloud droplet effective radius [um]
82 ✗ REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
83 REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
84
85 REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds
86 ! within the grid cell)
87
88 REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
89
90 REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
91 ! jq-end
92
93 ! cc PARAMETER (nexpo=1)
94
95 ! Calculer l'epaisseur optique et l'emmissivite des nuages
96
97 ✗ DO k = 1, klev
98 ✗ IF (iflag_t_glace.EQ.0) THEN
99 ✗ DO i = 1, klon
100 ✗ zfice(i) = 1.0 - (t(i,k)-t_glace_min_old)/(273.13-t_glace_min_old)
101 ✗ zfice(i) = min(max(zfice(i),0.0), 1.0)
102 ✗ zfice(i) = zfice(i)**exposant_glace_old
103 ENDDO
104 ELSE ! of IF (iflag_t_glace.EQ.0)
105 ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
106 ! zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
107 ! t_glace_max, exposant_glace)
108 ✗ CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
109 ENDIF
110
111 ✗ DO i = 1, klon
112 ✗ rad_chaud = rad_chau1
113 ✗ IF (k<=3) rad_chaud = rad_chau2
114
115 ✗ pclc(i, k) = max(pclc(i,k), seuil_neb)
116 ✗ zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))
117
118 ✗ IF (ok_aie) THEN
119 ! Formula "D" of Boucher and Lohmann, Tellus, 1995
120 !
121 cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
122 ✗ 1.E-4))/log(10.))*1.E6 !-m-3
123 ! Cloud droplet number concentration (CDNC) is restricted
124 ! to be within [20, 1000 cm^3]
125 !
126 ✗ cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
127 cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
128 ✗ 1.E-4))/log(10.))*1.E6 !-m-3
129 ✗ cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
130 !
131 !
132 ! air density: pplay(i,k) / (RD * zT(i,k))
133 ! factor 1.1: derive effective radius from volume-mean radius
134 ! factor 1000 is the water density
135 ! _chaud means that this is the CDR for liquid water clouds
136 !
137 rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &
138 ✗ *cdnc(i,k)))**(1./3.)
139 !
140 ! Convert to um. CDR shall be at least 3 um.
141 !
142 ✗ rad_chaud = max(rad_chaud*1.E6, 3.)
143
144 ! For output diagnostics
145 !
146 ! Cloud droplet effective radius [um]
147 !
148 ! we multiply here with f * xl (fraction of liquid water
149 ! clouds in the grid cell) to avoid problems in the
150 ! averaging of the output.
151 ! In the output of IOIPSL, derive the real cloud droplet
152 ! effective radius as re/fl
153 !
154 ✗ fl(i, k) = pclc(i, k)*(1.-zfice(i))
155 ✗ re(i, k) = rad_chaud*fl(i, k)
156
157 ! Pre-industrial cloud opt thickness
158 !
159 ! "radius" is calculated as rad_chaud above (plus the
160 ! ice cloud contribution) but using cdnc_pi instead of
161 ! cdnc.
162 radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &
163 ✗ 1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i)
164 ✗ cldtaupi(i, k) = 3.0/2.0*zflwp/radius
165 END IF ! ok_aie
166
167 ✗ radius = rad_chaud*(1.-zfice(i)) + rad_froid*zfice(i)
168 ✗ coef = coef_chau*(1.-zfice(i)) + coef_froi*zfice(i)
169 ✗ pcltau(i, k) = 3.0/2.0*zflwp/radius
170 ✗ pclemi(i, k) = 1.0 - exp(-coef*zflwp)
171 lo = (pclc(i,k)<=seuil_neb)
172 ✗ IF (lo) pclc(i, k) = 0.0
173 ✗ IF (lo) pcltau(i, k) = 0.0
174 ✗ IF (lo) pclemi(i, k) = 0.0
175
176 ✗ IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k)
177 END DO
178 END DO
179 ! cc DO k = 1, klev
180 ! cc DO i = 1, klon
181 ! cc t(i,k) = t(i,k)
182 ! cc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
183 ! cc lo = pclc(i,k) .GT. (2.*1.e-5)
184 ! cc zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
185 ! cc . /(rg*pclc(i,k))
186 ! cc zradef = 10.0 + (1.-sigs(k))*45.0
187 ! cc pcltau(i,k) = 1.5 * zflwp / zradef
188 ! cc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
189 ! cc zmsac = 0.13*(1.0-zfice) + 0.08*zfice
190 ! cc pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
191 ! cc if (.NOT.lo) pclc(i,k) = 0.0
192 ! cc if (.NOT.lo) pcltau(i,k) = 0.0
193 ! cc if (.NOT.lo) pclemi(i,k) = 0.0
194 ! cc ENDDO
195 ! cc ENDDO
196 ! ccccc print*, 'pas de nuage dans le rayonnement'
197 ! ccccc DO k = 1, klev
198 ! ccccc DO i = 1, klon
199 ! ccccc pclc(i,k) = 0.0
200 ! ccccc pcltau(i,k) = 0.0
201 ! ccccc pclemi(i,k) = 0.0
202 ! ccccc ENDDO
203 ! ccccc ENDDO
204
205 ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
206
207 ✗ DO i = 1, klon
208 ✗ pct(i) = 1.0
209 ✗ pch(i) = 1.0
210 ✗ pcm(i) = 1.0
211 ✗ pcl(i) = 1.0
212 ✗ pctlwp(i) = 0.0
213 END DO
214
215 ✗ DO k = klev, 1, -1
216 ✗ DO i = 1, klon
217 ✗ pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
218 ✗ pct(i) = pct(i)*(1.0-pclc(i,k))
219 ✗ IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))
220 ✗ IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &
221 ✗ pcm(i) = pcm(i)*(1.0-pclc(i,k))
222 ✗ IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))
223 END DO
224 END DO
225
226 ✗ DO i = 1, klon
227 ✗ pct(i) = 1. - pct(i)
228 ✗ pch(i) = 1. - pch(i)
229 ✗ pcm(i) = 1. - pcm(i)
230 ✗ pcl(i) = 1. - pcl(i)
231 END DO
232
233 ✗ RETURN
234 END SUBROUTINE nuage
235 ✗ SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
236 USE dimphy
237 IMPLICIT NONE
238
239 ! Laurent Li (LMD/CNRS), le 12 octobre 1998
240 ! (adaptation du code ECMWF)
241
242 ! Dans certains cas, le schema pronostique des nuages n'est
243 ! pas suffisament performant. On a donc besoin de diagnostiquer
244 ! ces nuages. Je dois avouer que c'est une frustration.
245
246 include "YOMCST.h"
247
248 ! Arguments d'entree:
249 REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
250 REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
251 REAL t(klon, klev) ! temperature (K)
252 REAL q(klon, klev) ! humidite specifique (Kg/Kg)
253 REAL rain(klon) ! pluie convective (kg/m2/s)
254 REAL snow(klon) ! neige convective (kg/m2/s)
255 INTEGER ktop(klon) ! sommet de la convection
256 INTEGER kbot(klon) ! bas de la convection
257
258 ! Arguments de sortie:
259 REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
260 REAL dialiq(klon, klev) ! eau liquide nuageuse
261
262 ! Constantes ajustables:
263 REAL canva, canvb, canvh
264 PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
265 REAL cca, ccb, ccc
266 PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
267 REAL ccfct, ccscal
268 PARAMETER (ccfct=0.400)
269 PARAMETER (ccscal=1.0E+11)
270 REAL cetahb, cetamb
271 PARAMETER (cetahb=0.45, cetamb=0.80)
272 REAL cclwmr
273 PARAMETER (cclwmr=1.E-04)
274 REAL zepscr
275 PARAMETER (zepscr=1.0E-10)
276
277 ! Variables locales:
278 INTEGER i, k
279 ✗ REAL zcc(klon)
280
281 ! Initialisation:
282
283 ✗ DO k = 1, klev
284 ✗ DO i = 1, klon
285 ✗ diafra(i, k) = 0.0
286 ✗ dialiq(i, k) = 0.0
287 END DO
288 END DO
289
290 ✗ DO i = 1, klon ! Calculer la fraction nuageuse
291 ✗ zcc(i) = 0.0
292 ✗ IF ((rain(i)+snow(i))>0.) THEN
293 ✗ zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
294 ✗ zcc(i) = min(ccc, max(0.0,zcc(i)))
295 END IF
296 END DO
297
298 ✗ DO i = 1, klon ! pour traiter les enclumes
299 ✗ diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
300 ✗ IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
301 1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
302 ✗ i)*ccfct,canva*(zcc(i)-canvb)))
303 ✗ dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
304 END DO
305
306 ✗ DO k = 1, klev ! nuages convectifs (sauf enclumes)
307 ✗ DO i = 1, klon
308 ✗ IF (k<ktop(i) .AND. k>=kbot(i)) THEN
309 ✗ diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
310 ✗ dialiq(i, k) = cclwmr*diafra(i, k)
311 END IF
312 END DO
313 END DO
314
315 ✗ RETURN
316 END SUBROUTINE diagcld1
317 ✗ SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
318 USE dimphy
319 IMPLICIT NONE
320
321 include "YOMCST.h"
322
323 ! Arguments d'entree:
324 REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
325 REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
326 REAL t(klon, klev) ! temperature (K)
327 REAL q(klon, klev) ! humidite specifique (Kg/Kg)
328
329 ! Arguments de sortie:
330 REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
331 REAL dialiq(klon, klev) ! eau liquide nuageuse
332
333 REAL cetamb
334 PARAMETER (cetamb=0.80)
335 REAL cloia, cloib, cloic, cloid
336 PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
337 ! cc PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
338 REAL rgammas
339 PARAMETER (rgammas=0.05)
340 REAL crhl
341 PARAMETER (crhl=0.15)
342 ! cc PARAMETER (CRHL=0.70)
343 REAL t_coup
344 PARAMETER (t_coup=234.0)
345
346 ! Variables locales:
347 ✗ INTEGER i, k, kb, invb(klon)
348 ✗ REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
349 REAL zdelta, zcor
350
351 ! Fonctions thermodynamiques:
352 include "YOETHF.h"
353 include "FCTTRE.h"
354
355 ! Initialisation:
356
357 ✗ DO k = 1, klev
358 ✗ DO i = 1, klon
359 ✗ diafra(i, k) = 0.0
360 ✗ dialiq(i, k) = 0.0
361 END DO
362 END DO
363
364 ✗ DO i = 1, klon
365 ✗ invb(i) = klev
366 ✗ zdthmin(i) = 0.0
367 END DO
368
369 ✗ DO k = 2, klev/2 - 1
370 ✗ DO i = 1, klon
371 zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
372 ✗ rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
373 ✗ zdthdp = zdthdp*cloia
374 ✗ IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
375 ✗ zdthmin(i) = zdthdp
376 ✗ invb(i) = k
377 END IF
378 END DO
379 END DO
380
381 ✗ DO i = 1, klon
382 ✗ kb = invb(i)
383 IF (thermcep) THEN
384 ✗ zdelta = max(0., sign(1.,rtt-t(i,kb)))
385 ✗ zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
386 ✗ zqs = min(0.5, zqs)
387 ✗ zcor = 1./(1.-retv*zqs)
388 ✗ zqs = zqs*zcor
389 ELSE
390 IF (t(i,kb)<t_coup) THEN
391 zqs = qsats(t(i,kb))/pplay(i, kb)
392 ELSE
393 zqs = qsatl(t(i,kb))/pplay(i, kb)
394 END IF
395 END IF
396 ✗ zcll = cloib*zdthmin(i) + cloic
397 ✗ zcll = min(1.0, max(0.0,zcll))
398 ✗ zrhb = q(i, kb)/zqs
399 ✗ IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
400 ✗ zcll = min(1.0, max(0.0,zcll))
401 ✗ diafra(i, kb) = max(diafra(i,kb), zcll)
402 ✗ dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
403 END DO
404
405 ✗ RETURN
406 END SUBROUTINE diagcld2
407