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