/[lmdze]/trunk/libf/phylmd/nuage.f
ViewVC logotype

Annotation of /trunk/libf/phylmd/nuage.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (hide annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 4 months ago) by guez
File size: 14033 byte(s)
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

Variables of module "yoethf" were defined in procedure "suphec". Moved
these definitions to a new procedure "yoethf" in module "yoethf_m".

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

  ViewVC Help
Powered by ViewVC 1.1.21