source: lmdz_wrf/WRFV3/lmdz/nuage.F90 @ 1

Last change on this file since 1 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!
3      SUBROUTINE nuage (paprs, pplay,                                                &
4       &                  t, pqlwp, pclc, pcltau, pclemi,                            &
5       &                  pch, pcl, pcm, pct, pctlwp,                                &
6       &                  ok_aie,                                                    &
7       &                  mass_solu_aero, mass_solu_aero_pi,                         &
8       &                  bl95_b0, bl95_b1,                                          &
9       &                  cldtaupi, re, fl)
10      USE dimphy
11      IMPLICIT none
12!c======================================================================
13!c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
14!c Objet: Calculer epaisseur optique et emmissivite des nuages
15!c======================================================================
16!c Arguments:
17!c t-------input-R-temperature
18!c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
19!c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
20!c ok_aie--input-L-apply aerosol indirect effect or not
21!c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
22!c mass_solu_aero_pi--input-R-dito, pre-industrial value
23!c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
24!c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
25!c     
26!c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
27!c                   needed for the diagnostics of the aerosol indirect
28!c                   radiative forcing (see radlwsw)
29!c re------output-R-Cloud droplet effective radius multiplied by fl [um]
30!c fl------output-R-Denominator to re, introduced to avoid problems in
31!c                  the averaging of the output. fl is the fraction of liquid
32!c                  water clouds within a grid cell     
33!c
34!c pcltau--output-R-epaisseur optique des nuages
35!c pclemi--output-R-emissivite des nuages (0 a 1)
36!c======================================================================
37!C
38#include "YOMCST.h"
39!c
40!cym#include "dimensions.h"
41!cym#include "dimphy.h"
42      REAL paprs(klon,klev+1), pplay(klon,klev)
43      REAL t(klon,klev)
44!c
45      REAL pclc(klon,klev)
46      REAL pqlwp(klon,klev)
47      REAL pcltau(klon,klev), pclemi(klon,klev)
48!c
49      REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
50!c
51      LOGICAL lo
52!c
53      REAL cetahb, cetamb
54      PARAMETER (cetahb = 0.45, cetamb = 0.80)
55!C
56      INTEGER i, k
57      REAL zflwp, zradef, zfice, zmsac
58!c
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)
61!ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
62!c 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     
70!cjq for the aerosol indirect effect
71!cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
72!cjq     
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
87!cjq-end     
88     
89!ccc      PARAMETER (nexpo=1)
90!c
91!c Calculer l'epaisseur optique et l'emmissivite des nuages
92!c
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
169!ccc      DO k = 1, klev
170!ccc      DO i = 1, klon
171!ccc         t(i,k) = t(i,k)
172!ccc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
173!ccc         lo = pclc(i,k) .GT. (2.*1.e-5)
174!ccc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
175!ccc     .          /(rg*pclc(i,k))
176!ccc         zradef = 10.0 + (1.-sigs(k))*45.0
177!ccc         pcltau(i,k) = 1.5 * zflwp / zradef
178!ccc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
179!ccc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
180!ccc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
181!ccc         if (.NOT.lo) pclc(i,k) = 0.0
182!ccc         if (.NOT.lo) pcltau(i,k) = 0.0
183!ccc         if (.NOT.lo) pclemi(i,k) = 0.0
184!ccc      ENDDO
185!ccc      ENDDO
186!cccccc      print*, 'pas de nuage dans le rayonnement'
187!cccccc      DO k = 1, klev
188!cccccc      DO i = 1, klon
189!cccccc         pclc(i,k) = 0.0
190!cccccc         pcltau(i,k) = 0.0
191!cccccc         pclemi(i,k) = 0.0
192!cccccc      ENDDO
193!cccccc      ENDDO
194!C
195!C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
196!C
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
204!C
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
219!C
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
226!C
227      RETURN
228      END SUBROUTINE nuage
229
230      SUBROUTINE diagcld1(paprs,pplay,rain,snow,kbot,ktop,                           &
231       &                   diafra,dialiq)
232      use dimphy
233      IMPLICIT none
234!c
235!c Laurent Li (LMD/CNRS), le 12 octobre 1998
236!c                        (adaptation du code ECMWF)
237!c
238!c Dans certains cas, le schema pronostique des nuages n'est
239!c pas suffisament performant. On a donc besoin de diagnostiquer
240!c ces nuages. Je dois avouer que c'est une frustration.
241!c
242!cym#include "dimensions.h"
243!cym#include "dimphy.h"
244#include "YOMCST.h"
245!c
246!c Arguments d'entree:
247      REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
248      REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
249      REAL t(klon,klev) ! temperature (K)
250      REAL q(klon,klev) ! humidite specifique (Kg/Kg)
251      REAL rain(klon) ! pluie convective (kg/m2/s)
252      REAL snow(klon) ! neige convective (kg/m2/s)
253      INTEGER ktop(klon) ! sommet de la convection
254      INTEGER kbot(klon) ! bas de la convection
255!c
256!c Arguments de sortie:
257      REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
258      REAL dialiq(klon,klev) ! eau liquide nuageuse
259!c
260!c Constantes ajustables:
261      REAL CANVA, CANVB, CANVH
262      PARAMETER (CANVA=2.0, CANVB=0.3, CANVH=0.4)
263      REAL CCA, CCB, CCC
264      PARAMETER (CCA=0.125, CCB=1.5, CCC=0.8)
265      REAL CCFCT, CCSCAL
266      PARAMETER (CCFCT=0.400)
267      PARAMETER (CCSCAL=1.0E+11)
268      REAL CETAHB, CETAMB
269      PARAMETER (CETAHB=0.45, CETAMB=0.80)
270      REAL CCLWMR
271      PARAMETER (CCLWMR=1.E-04)
272      REAL ZEPSCR
273      PARAMETER (ZEPSCR=1.0E-10)
274!c
275!c Variables locales:
276      INTEGER i, k
277      REAL zcc(klon)
278!c
279!c Initialisation:
280!c
281      DO k = 1, klev
282      DO i = 1, klon
283         diafra(i,k) = 0.0
284         dialiq(i,k) = 0.0
285      ENDDO
286      ENDDO
287!c
288      DO i = 1, klon ! Calculer la fraction nuageuse
289      zcc(i) = 0.0
290      IF((rain(i)+snow(i)).GT.0.) THEN
291         zcc(i)= CCA * LOG(MAX(ZEPSCR,(rain(i)+snow(i))*CCSCAL))-CCB
292         zcc(i)= MIN(CCC,MAX(0.0,zcc(i)))
293      ENDIF
294      ENDDO
295!c
296      DO i = 1, klon ! pour traiter les enclumes
297      diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),zcc(i)*CCFCT)
298      IF ((zcc(i).GE.CANVH) .AND.                                                    &
299       &    (pplay(i,ktop(i)).LE.CETAHB*paprs(i,1)))                                 &
300       & diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),                                  &
301       &                         MAX(zcc(i)*CCFCT,CANVA*(zcc(i)-CANVB)))
302      dialiq(i,ktop(i))=CCLWMR*diafra(i,ktop(i))
303      ENDDO
304!c
305      DO k = 1, klev ! nuages convectifs (sauf enclumes)
306      DO i = 1, klon
307      IF (k.LT.ktop(i) .AND. k.GE.kbot(i)) THEN
308         diafra(i,k)=MAX(diafra(i,k),zcc(i)*CCFCT)
309         dialiq(i,k)=CCLWMR*diafra(i,k)
310      ENDIF
311      ENDDO
312      ENDDO
313!c
314      RETURN
315      END SUBROUTINE diagcld1
316
317      SUBROUTINE diagcld2(paprs,pplay,t,q, diafra,dialiq)
318      use dimphy
319      IMPLICIT none
320!c
321!cym#include "dimensions.h"
322!cym#include "dimphy.h"
323#include "YOMCST.h"
324!c
325!c Arguments d'entree:
326      REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
327      REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
328      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 Fonctions thermodynamiques:
354#include "YOETHF.h"
355#include "FCTTRE.h"
356!c
357!c Initialisation:
358!c
359      DO k = 1, klev
360      DO i = 1, klon
361         diafra(i,k) = 0.0
362         dialiq(i,k) = 0.0
363      ENDDO
364      ENDDO
365!c
366      DO i = 1, klon
367         invb(i) = klev
368         zdthmin(i)=0.0
369      ENDDO
370
371      DO k = 2, klev/2-1
372      DO i = 1, klon
373         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))                        &
374       &          - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
375         zdthdp = zdthdp * CLOIA
376         IF (pplay(i,k).GT.CETAMB*paprs(i,1) .AND.                                   &
377       &       zdthdp.LT.zdthmin(i) ) THEN
378            zdthmin(i) = zdthdp
379            invb(i) = k
380         ENDIF
381      ENDDO
382      ENDDO
383
384      DO i = 1, klon
385         kb=invb(i)
386         IF (thermcep) THEN
387           zdelta=MAX(0.,SIGN(1.,RTT-t(i,kb)))
388           zqs= R2ES*FOEEW(t(i,kb),zdelta)/pplay(i,kb)
389           zqs=MIN(0.5,zqs)
390           zcor=1./(1.-RETV*zqs)
391           zqs=zqs*zcor
392         ELSE
393           IF (t(i,kb) .LT. t_coup) THEN
394              zqs = qsats(t(i,kb)) / pplay(i,kb)
395           ELSE
396              zqs = qsatl(t(i,kb)) / pplay(i,kb)
397           ENDIF
398         ENDIF
399         zcll = CLOIB * zdthmin(i) + CLOIC
400         zcll = MIN(1.0,MAX(0.0,zcll))
401         zrhb= q(i,kb)/zqs
402         IF (zcll.GT.0.0.AND.zrhb.LT.CRHL)                                           &
403       &   zcll=zcll*(1.-(CRHL-zrhb)*CLOID)
404         zcll=MIN(1.0,MAX(0.0,zcll))
405         diafra(i,kb) = MAX(diafra(i,kb),zcll)
406         dialiq(i,kb)= diafra(i,kb) * RGAMMAS*zqs
407      ENDDO
408!c
409      RETURN
410      END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.