source: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/nuage.F @ 4126

Last change on this file since 4126 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.7 KB
Line 
1! $Id: nuage.F 1279 2009-12-10 09:02:56Z fairhead $
2!
3      SUBROUTINE nuage (paprs, pplay,
4     .                  t, pqlwp, pclc, pcltau, pclemi,
5     .                  pch, pcl, pcm, pct, pctlwp,
6     e                  ok_aie,
7     e                  mass_solu_aero, mass_solu_aero_pi,
8     e                  bl95_b0, bl95_b1,
9     s                  cldtaupi, re, fl)
10      USE dimphy
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 mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
22c mass_solu_aero_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
40cym#include "dimensions.h"
41cym#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 mass_solu_aero(klon, klev)    ! total mass concentration for all soluble aerosols[ug m-3]
76      REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
77      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
78      REAL re(klon, klev)       ! cloud droplet effective radius [um]
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(mass_solu_aero(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(mass_solu_aero_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      use dimphy
232      IMPLICIT none
233c
234c Laurent Li (LMD/CNRS), le 12 octobre 1998
235c                        (adaptation du code ECMWF)
236c
237c Dans certains cas, le schema pronostique des nuages n'est
238c pas suffisament performant. On a donc besoin de diagnostiquer
239c ces nuages. Je dois avouer que c'est une frustration.
240c
241cym#include "dimensions.h"
242cym#include "dimphy.h"
243#include "YOMCST.h"
244c
245c Arguments d'entree:
246      REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
247      REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
248      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
254c
255c Arguments de sortie:
256      REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
257      REAL dialiq(klon,klev) ! eau liquide nuageuse
258c
259c 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)
273c
274c Variables locales:
275      INTEGER i, k
276      REAL zcc(klon)
277c
278c Initialisation:
279c
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
286c
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
294c
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
303c
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
312c
313      RETURN
314      END
315      SUBROUTINE diagcld2(paprs,pplay,t,q, diafra,dialiq)
316      use dimphy
317      IMPLICIT none
318c
319cym#include "dimensions.h"
320cym#include "dimphy.h"
321#include "YOMCST.h"
322c
323c 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)
328c
329c Arguments de sortie:
330      REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
331      REAL dialiq(klon,klev) ! eau liquide nuageuse
332c
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)
337ccc      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)
342ccc      PARAMETER (CRHL=0.70)
343      REAL t_coup
344      PARAMETER (t_coup=234.0)
345c
346c Variables locales:
347      INTEGER i, k, kb, invb(klon)
348      REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
349      REAL zdelta, zcor
350c
351c Fonctions thermodynamiques:
352#include "YOETHF.h"
353#include "FCTTRE.h"
354c
355c Initialisation:
356c
357      DO k = 1, klev
358      DO i = 1, klon
359         diafra(i,k) = 0.0
360         dialiq(i,k) = 0.0
361      ENDDO
362      ENDDO
363c
364      DO i = 1, klon
365         invb(i) = klev
366         zdthmin(i)=0.0
367      ENDDO
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).GT.CETAMB*paprs(i,1) .AND.
375     .       zdthdp.LT.zdthmin(i) ) THEN
376            zdthmin(i) = zdthdp
377            invb(i) = k
378         ENDIF
379      ENDDO
380      ENDDO
381
382      DO i = 1, klon
383         kb=invb(i)
384         IF (thermcep) THEN
385           zdelta=MAX(0.,SIGN(1.,RTT-t(i,kb)))
386           zqs= R2ES*FOEEW(t(i,kb),zdelta)/pplay(i,kb)
387           zqs=MIN(0.5,zqs)
388           zcor=1./(1.-RETV*zqs)
389           zqs=zqs*zcor
390         ELSE
391           IF (t(i,kb) .LT. t_coup) THEN
392              zqs = qsats(t(i,kb)) / pplay(i,kb)
393           ELSE
394              zqs = qsatl(t(i,kb)) / pplay(i,kb)
395           ENDIF
396         ENDIF
397         zcll = CLOIB * zdthmin(i) + CLOIC
398         zcll = MIN(1.0,MAX(0.0,zcll))
399         zrhb= q(i,kb)/zqs
400         IF (zcll.GT.0.0.AND.zrhb.LT.CRHL)
401     .   zcll=zcll*(1.-(CRHL-zrhb)*CLOID)
402         zcll=MIN(1.0,MAX(0.0,zcll))
403         diafra(i,kb) = MAX(diafra(i,kb),zcll)
404         dialiq(i,kb)= diafra(i,kb) * RGAMMAS*zqs
405      ENDDO
406c
407      RETURN
408      END
Note: See TracBrowser for help on using the repository browser.