source: LMDZ4/branches/V3_test/libf/phylmd/nuage.F @ 703

Last change on this file since 703 was 703, checked in by (none), 19 years ago

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