source: LMDZ4/trunk/libf/phylmd/newmicro.F @ 758

Last change on this file since 758 was 685, checked in by lmdzadmin, 19 years ago

Prise en compte du type de recouvrement cf. CR
IM
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.3 KB
RevLine 
[524]1!
2! $Header$
3!
4      SUBROUTINE newmicro (paprs, pplay,ok_newmicro,
5     .                  t, pqlwp, pclc, pcltau, pclemi,
6     .                  pch, pcl, pcm, pct, pctlwp,
7     s                  xflwp, xfiwp, xflwc, xfiwc,
8     e                  ok_aie,
9     e                  sulfate, sulfate_pi,
10     e                  bl95_b0, bl95_b1,
11     s                  cldtaupi, re, fl)
12      IMPLICIT none
13c======================================================================
14c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
15c Objet: Calculer epaisseur optique et emmissivite des nuages
16c======================================================================
17c Arguments:
18c t-------input-R-temperature
19c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
20c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
21c
22c ok_aie--input-L-apply aerosol indirect effect or not
23c sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
24c sulfate_pi-input-R-dito, pre-industrial value
25c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
26c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
27c     
28c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
29c                   needed for the diagnostics of the aerosol indirect
30c                   radiative forcing (see radlwsw)
31c re------output-R-Cloud droplet effective radius multiplied by fl [um]
32c fl------output-R-Denominator to re, introduced to avoid problems in
33c                  the averaging of the output. fl is the fraction of liquid
34c                  water clouds within a grid cell           
35c pcltau--output-R-epaisseur optique des nuages
36c pclemi--output-R-emissivite des nuages (0 a 1)
37c======================================================================
38C
39#include "YOMCST.h"
40c
41#include "dimensions.h"
42#include "dimphy.h"
43#include "nuage.h"
[685]44cIM cf. CR: include pour NOVLP et ZEPSEC
45#include "radepsi.h"
46#include "radopt.h"
[524]47      REAL paprs(klon,klev+1), pplay(klon,klev)
48      REAL t(klon,klev)
49c
50      REAL pclc(klon,klev)
51      REAL pqlwp(klon,klev)
52      REAL pcltau(klon,klev), pclemi(klon,klev)
53c
54      REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
55c
56      LOGICAL lo
57c
58      REAL cetahb, cetamb
59      PARAMETER (cetahb = 0.45, cetamb = 0.80)
60C
61      INTEGER i, k
62cIM: 091003   REAL zflwp, zradef, zfice, zmsac
63      REAL zflwp(klon), zradef, zfice, zmsac
64cIM: 091003 rajout
65      REAL xflwp(klon), xfiwp(klon)
66      REAL xflwc(klon,klev), xfiwc(klon,klev)
67c
68      REAL radius, rad_chaud
69cc      PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
70ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
71c sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
72      REAL coef, coef_froi, coef_chau
73      PARAMETER (coef_chau=0.13, coef_froi=0.09)
74      REAL seuil_neb, t_glace
75      PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
76      INTEGER nexpo ! exponentiel pour glace/eau
77      PARAMETER (nexpo=6)
78ccc      PARAMETER (nexpo=1)
79
80c -- sb:
81      logical ok_newmicro
82c     parameter (ok_newmicro=.FALSE.)
83cIM: 091003   real rel, tc, rei, zfiwp
84      real rel, tc, rei, zfiwp(klon)
85      real k_liq, k_ice0, k_ice, DF
86      parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g
87      parameter (DF=1.66) ! diffusivity factor
88c sb --
89cjq for the aerosol indirect effect
90cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
91cjq     
92      LOGICAL ok_aie            ! Apply AIE or not?
93      LOGICAL ok_a1lwpdep       ! a1 LWP dependent?
94     
95      REAL sulfate(klon, klev)  ! sulfate aerosol mass concentration [ug m-3]
96      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
97      REAL re(klon, klev)       ! cloud droplet effective radius [um]
98      REAL sulfate_pi(klon, klev)  ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
99      REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
100      REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
101     
102      REAL fl(klon, klev)       ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
103     
104      REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
105     
106      REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
107cjq-end   
[685]108cIM cf. CR:parametres supplementaires
109      REAL zclear(klon)
110      REAL zcloud(klon)
[524]111c
112c Calculer l'epaisseur optique et l'emmissivite des nuages
113c
114cIM inversion des DO
115      DO i = 1, klon
116       xflwp(i)=0.
117       xfiwp(i)=0.
118      DO k = 1, klev
119c
120       xflwc(i,k)=0.
121       xfiwc(i,k)=0.
122c
123         rad_chaud = rad_chau1
124         IF (k.LE.3) rad_chaud = rad_chau2
125         pclc(i,k) = MAX(pclc(i,k), seuil_neb)
126         zflwp(i) = 1000.*pqlwp(i,k)/RG/pclc(i,k)
127     .          *(paprs(i,k)-paprs(i,k+1))
128         zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
129         zfice = MIN(MAX(zfice,0.0),1.0)
130         zfice = zfice**nexpo
131         radius = rad_chaud * (1.-zfice) + rad_froid * zfice
132         coef = coef_chau * (1.-zfice) + coef_froi * zfice
133         pcltau(i,k) = 3.0/2.0 * zflwp(i) / radius
134         pclemi(i,k) = 1.0 - EXP( - coef * zflwp(i))
135
136         if (ok_newmicro) then
137
138c -- liquid/ice cloud water paths:
139
140         zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
141         zfice = MIN(MAX(zfice,0.0),1.0)
142
143         zflwp(i) = 1000.*(1.-zfice)*pqlwp(i,k)/pclc(i,k)
144     :          *(paprs(i,k)-paprs(i,k+1))/RG
145         zfiwp(i) = 1000.*zfice*pqlwp(i,k)/pclc(i,k)
146     :          *(paprs(i,k)-paprs(i,k+1))/RG
147
148         xflwp(i) = xflwp(i)+ (1.-zfice)*pqlwp(i,k)
149     :          *(paprs(i,k)-paprs(i,k+1))/RG
150         xfiwp(i) = xfiwp(i)+ zfice*pqlwp(i,k)
151     :          *(paprs(i,k)-paprs(i,k+1))/RG
152
153cIM Total Liquid/Ice water content
154         xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)
155         xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)
156cIM In-Cloud Liquid/Ice water content
157c        xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
158c        xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
159
160c -- effective cloud droplet radius (microns):
161
162c for liquid water clouds:
163         IF (ok_aie) THEN
164            ! Formula "D" of Boucher and Lohmann, Tellus, 1995
165            !             
166            cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
167     .           log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
168            ! Cloud droplet number concentration (CDNC) is restricted
169            ! to be within [20, 1000 cm^3]
170            !
171            cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
172            !
173            !
174            cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
175     .           log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
176            cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
177            !           
178            !
179            ! air density: pplay(i,k) / (RD * zT(i,k))
180            ! factor 1.1: derive effective radius from volume-mean radius
181            ! factor 1000 is the water density
182            ! _chaud means that this is the CDR for liquid water clouds
183            !
184            rad_chaud =
185     .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) ) 
186     .               / (4./3. * RPI * 1000. * cdnc(i,k)) )**(1./3.)
187            !
188            ! Convert to um. CDR shall be at least 3 um.
189            !
190c           rad_chaud = MAX(rad_chaud*1.e6, 3.)
191            rad_chaud = MAX(rad_chaud*1.e6, 5.)
192           
193            ! Pre-industrial cloud opt thickness
194            !
195            ! "radius" is calculated as rad_chaud above (plus the
196            ! ice cloud contribution) but using cdnc_pi instead of
197            ! cdnc.
198            radius =
199     .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) ) 
200     .               / (4./3. * RPI * 1000. * cdnc_pi(i,k)) )**(1./3.)
[534]201            radius = MAX(radius*1.e6, 5.)
[524]202           
203            tc = t(i,k)-273.15
204            rei = 0.71*tc + 61.29
205            if (tc.le.-81.4) rei = 3.5
206            if (zflwp(i).eq.0.) radius = 1.
207            if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1.
208            cldtaupi(i,k) = 3.0/2.0 * zflwp(i) / radius
209     .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
210         ENDIF                  ! ok_aie
211         ! For output diagnostics
212         !
213         ! Cloud droplet effective radius [um]
214         !
215         ! we multiply here with f * xl (fraction of liquid water
216         ! clouds in the grid cell) to avoid problems in the
217         ! averaging of the output.
218         ! In the output of IOIPSL, derive the real cloud droplet
219         ! effective radius as re/fl
220         !
221         fl(i,k) = pclc(i,k)*(1.-zfice)           
222         re(i,k) = rad_chaud*fl(i,k)
223           
224c-jq end         
225         
226         rel = rad_chaud
227c for ice clouds: as a function of the ambiant temperature
228c [formula used by Iacobellis and Somerville (2000), with an
229c asymptotical value of 3.5 microns at T<-81.4 C added to be
230c consistent with observations of Heymsfield et al. 1986]:
231         tc = t(i,k)-273.15
232         rei = 0.71*tc + 61.29
233         if (tc.le.-81.4) rei = 3.5
234
235c -- cloud optical thickness :
236
237c [for liquid clouds, traditional formula,
238c  for ice clouds, Ebert & Curry (1992)]
239
240         if (zflwp(i).eq.0.) rel = 1.
241         if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1.
242         pcltau(i,k) = 3.0/2.0 * ( zflwp(i)/rel )
243     .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
244
245c -- cloud infrared emissivity:
246
247c [the broadband infrared absorption coefficient is parameterized
248c  as a function of the effective cld droplet radius]
249
250c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
251         k_ice = k_ice0 + 1.0/rei
252
253         pclemi(i,k) = 1.0
254     .      - EXP( - coef_chau*zflwp(i) - DF*k_ice*zfiwp(i) )
255
256         endif ! ok_newmicro
257
258         lo = (pclc(i,k) .LE. seuil_neb)
259         IF (lo) pclc(i,k) = 0.0
260         IF (lo) pcltau(i,k) = 0.0
261         IF (lo) pclemi(i,k) = 0.0
262         
263         IF (lo) cldtaupi(i,k) = 0.0
264         IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)           
265      ENDDO
266      ENDDO
267ccc      DO k = 1, klev
268ccc      DO i = 1, klon
269ccc         t(i,k) = t(i,k)
270ccc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
271ccc         lo = pclc(i,k) .GT. (2.*1.e-5)
272ccc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
273ccc     .          /(rg*pclc(i,k))
274ccc         zradef = 10.0 + (1.-sigs(k))*45.0
275ccc         pcltau(i,k) = 1.5 * zflwp / zradef
276ccc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
277ccc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
278ccc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
279ccc         if (.NOT.lo) pclc(i,k) = 0.0
280ccc         if (.NOT.lo) pcltau(i,k) = 0.0
281ccc         if (.NOT.lo) pclemi(i,k) = 0.0
282ccc      ENDDO
283ccc      ENDDO
284cccccc      print*, 'pas de nuage dans le rayonnement'
285cccccc      DO k = 1, klev
286cccccc      DO i = 1, klon
287cccccc         pclc(i,k) = 0.0
288cccccc         pcltau(i,k) = 0.0
289cccccc         pclemi(i,k) = 0.0
290cccccc      ENDDO
291cccccc      ENDDO
292C
293C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
294C
[685]295cIM cf. CR:test: calcul prenant ou non en compte le recouvrement
296cinitialisations
297      DO i=1,klon
298         zclear(i)=1.
299         zcloud(i)=0.
[524]300         pch(i)=1.0
301         pcm(i) = 1.0
302         pcl(i) = 1.0
303         pctlwp(i) = 0.0
304      ENDDO
305C
[685]306cIM cf CR DO k=1,klev
[524]307      DO k = klev, 1, -1
308      DO i = 1, klon
309         pctlwp(i) = pctlwp(i)
310     .             + pqlwp(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
[685]311cIM cf. CR
312            IF (NOVLP.EQ.1) THEN
313               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
314     s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
315               pct(i)=1.-zclear(i)
316               if (pplay(i,k).LE.cetahb*paprs(i,1)) then
317                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
318     s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
319               else if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
320     .                  pplay(i,k).LE.cetamb*paprs(i,1)) then
321                  pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
322     s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
323               else if (pplay(i,k).GT.cetamb*paprs(i,1)) then
324                  pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
325     s                            /(1.-MIN(zcloud(i),1.-ZEPSEC))
326               endif
327               zcloud(i)=pclc(i,k)
328            ELSE IF (NOVLP.EQ.2) THEN
329               zcloud(i)=MAX(pclc(i,k),zcloud(i))
330               pct(i)=zcloud(i)
331               if (pplay(i,k).LE.cetahb*paprs(i,1)) then
332                  pch(i) = MIN(pclc(i,k),pch(i))
333               else if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
334     .                  pplay(i,k).LE.cetamb*paprs(i,1)) then
335                  pcm(i) = MIN(pclc(i,k),pcm(i))
336               else if (pplay(i,k).GT.cetamb*paprs(i,1)) then
337                  pcl(i) = MIN(pclc(i,k),pcl(i))
338               endif
339            ELSE IF (NOVLP.EQ.3) THEN
340               zclear(i)=zclear(i)*(1.-pclc(i,k))
341               pct(i)=1-zclear(i)
342               if (pplay(i,k).LE.cetahb*paprs(i,1)) then
343               pch(i) = pch(i)*(1.0-pclc(i,k))
344               else if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
345     .                  pplay(i,k).LE.cetamb*paprs(i,1)) then
346               pcm(i) = pcm(i)*(1.0-pclc(i,k))
347               else if (pplay(i,k).GT.cetamb*paprs(i,1)) then
348               pcl(i) = pcl(i)*(1.0-pclc(i,k))
349               endif
350            ENDIF
351         ENDDO
[524]352      ENDDO
353C
354      DO i = 1, klon
[685]355cIM cf. CR          pct(i)=1.-pct(i)
[524]356         pch(i)=1.-pch(i)
357         pcm(i)=1.-pcm(i)
358         pcl(i)=1.-pcl(i)
359      ENDDO
360C
361      RETURN
362      END
Note: See TracBrowser for help on using the repository browser.