source: lmdz_wrf/trunk/WRFV3/lmdz/nuage_mod.F90 @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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