source: LMDZ4/branches/IPSL-CM4_IPCC_patches/libf/phylmd/newmicro.F @ 5036

Last change on this file since 5036 was 540, checked in by (none), 20 years ago

This commit was manufactured by cvs2svn to create branch
'IPSL-CM4_IPCC_patches'.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.4 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"
44      REAL paprs(klon,klev+1), pplay(klon,klev)
45      REAL t(klon,klev)
46c
47      REAL pclc(klon,klev)
48      REAL pqlwp(klon,klev)
49      REAL pcltau(klon,klev), pclemi(klon,klev)
50c
51      REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
52c
53      LOGICAL lo
54c
55      REAL cetahb, cetamb
56      PARAMETER (cetahb = 0.45, cetamb = 0.80)
57C
58      INTEGER i, k
59cIM: 091003   REAL zflwp, zradef, zfice, zmsac
60      REAL zflwp(klon), zradef, zfice, zmsac
61cIM: 091003 rajout
62      REAL xflwp(klon), xfiwp(klon)
63      REAL xflwc(klon,klev), xfiwc(klon,klev)
64c
65      REAL radius, rad_chaud
66cc      PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
67ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
68c sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
69      REAL coef, coef_froi, coef_chau
70      PARAMETER (coef_chau=0.13, coef_froi=0.09)
71      REAL seuil_neb, t_glace
72      PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
73      INTEGER nexpo ! exponentiel pour glace/eau
74      PARAMETER (nexpo=6)
75ccc      PARAMETER (nexpo=1)
76
77c -- sb:
78      logical ok_newmicro
79c     parameter (ok_newmicro=.FALSE.)
80cIM: 091003   real rel, tc, rei, zfiwp
81      real rel, tc, rei, zfiwp(klon)
82      real k_liq, k_ice0, k_ice, DF
83      parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g
84      parameter (DF=1.66) ! diffusivity factor
85c sb --
86cjq for the aerosol indirect effect
87cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
88cjq     
89      LOGICAL ok_aie            ! Apply AIE or not?
90      LOGICAL ok_a1lwpdep       ! a1 LWP dependent?
91     
92      REAL sulfate(klon, klev)  ! sulfate aerosol mass concentration [ug m-3]
93      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
94      REAL re(klon, klev)       ! cloud droplet effective radius [um]
95      REAL sulfate_pi(klon, klev)  ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
96      REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
97      REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
98     
99      REAL fl(klon, klev)       ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
100     
101      REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
102     
103      REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
104cjq-end   
105c
106c Calculer l'epaisseur optique et l'emmissivite des nuages
107c
108cIM inversion des DO
109      DO i = 1, klon
110       xflwp(i)=0.
111       xfiwp(i)=0.
112      DO k = 1, klev
113c
114       xflwc(i,k)=0.
115       xfiwc(i,k)=0.
116c
117         rad_chaud = rad_chau1
118         IF (k.LE.3) rad_chaud = rad_chau2
119         pclc(i,k) = MAX(pclc(i,k), seuil_neb)
120         zflwp(i) = 1000.*pqlwp(i,k)/RG/pclc(i,k)
121     .          *(paprs(i,k)-paprs(i,k+1))
122         zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
123         zfice = MIN(MAX(zfice,0.0),1.0)
124         zfice = zfice**nexpo
125         radius = rad_chaud * (1.-zfice) + rad_froid * zfice
126         coef = coef_chau * (1.-zfice) + coef_froi * zfice
127         pcltau(i,k) = 3.0/2.0 * zflwp(i) / radius
128         pclemi(i,k) = 1.0 - EXP( - coef * zflwp(i))
129
130         if (ok_newmicro) then
131
132c -- liquid/ice cloud water paths:
133
134         zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
135         zfice = MIN(MAX(zfice,0.0),1.0)
136
137         zflwp(i) = 1000.*(1.-zfice)*pqlwp(i,k)/pclc(i,k)
138     :          *(paprs(i,k)-paprs(i,k+1))/RG
139         zfiwp(i) = 1000.*zfice*pqlwp(i,k)/pclc(i,k)
140     :          *(paprs(i,k)-paprs(i,k+1))/RG
141
142         xflwp(i) = xflwp(i)+ (1.-zfice)*pqlwp(i,k)
143     :          *(paprs(i,k)-paprs(i,k+1))/RG
144         xfiwp(i) = xfiwp(i)+ zfice*pqlwp(i,k)
145     :          *(paprs(i,k)-paprs(i,k+1))/RG
146
147cIM Total Liquid/Ice water content
148         xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)
149         xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)
150cIM In-Cloud Liquid/Ice water content
151c        xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
152c        xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
153
154c -- effective cloud droplet radius (microns):
155
156c for liquid water clouds:
157         IF (ok_aie) THEN
158            ! Formula "D" of Boucher and Lohmann, Tellus, 1995
159            !             
160            cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
161     .           log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
162            ! Cloud droplet number concentration (CDNC) is restricted
163            ! to be within [20, 1000 cm^3]
164            !
165            cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
166            !
167            !
168            cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
169     .           log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
170            cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
171            !           
172            !
173            ! air density: pplay(i,k) / (RD * zT(i,k))
174            ! factor 1.1: derive effective radius from volume-mean radius
175            ! factor 1000 is the water density
176            ! _chaud means that this is the CDR for liquid water clouds
177            !
178            rad_chaud =
179     .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) ) 
180     .               / (4./3. * RPI * 1000. * cdnc(i,k)) )**(1./3.)
181            !
182            ! Convert to um. CDR shall be at least 3 um.
183            !
184c           rad_chaud = MAX(rad_chaud*1.e6, 3.)
185            rad_chaud = MAX(rad_chaud*1.e6, 5.)
186           
187            ! Pre-industrial cloud opt thickness
188            !
189            ! "radius" is calculated as rad_chaud above (plus the
190            ! ice cloud contribution) but using cdnc_pi instead of
191            ! cdnc.
192            radius =
193     .           1.1 * ( (pqlwp(i,k) * pplay(i,k) / (RD * T(i,k)) ) 
194     .               / (4./3. * RPI * 1000. * cdnc_pi(i,k)) )**(1./3.)
[534]195            radius = MAX(radius*1.e6, 5.)
[524]196           
197            tc = t(i,k)-273.15
198            rei = 0.71*tc + 61.29
199            if (tc.le.-81.4) rei = 3.5
200            if (zflwp(i).eq.0.) radius = 1.
201            if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1.
202            cldtaupi(i,k) = 3.0/2.0 * zflwp(i) / radius
203     .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
204         ENDIF                  ! ok_aie
205         ! For output diagnostics
206         !
207         ! Cloud droplet effective radius [um]
208         !
209         ! we multiply here with f * xl (fraction of liquid water
210         ! clouds in the grid cell) to avoid problems in the
211         ! averaging of the output.
212         ! In the output of IOIPSL, derive the real cloud droplet
213         ! effective radius as re/fl
214         !
215         fl(i,k) = pclc(i,k)*(1.-zfice)           
216         re(i,k) = rad_chaud*fl(i,k)
217           
218c-jq end         
219         
220         rel = rad_chaud
221c for ice clouds: as a function of the ambiant temperature
222c [formula used by Iacobellis and Somerville (2000), with an
223c asymptotical value of 3.5 microns at T<-81.4 C added to be
224c consistent with observations of Heymsfield et al. 1986]:
225         tc = t(i,k)-273.15
226         rei = 0.71*tc + 61.29
227         if (tc.le.-81.4) rei = 3.5
228
229c -- cloud optical thickness :
230
231c [for liquid clouds, traditional formula,
232c  for ice clouds, Ebert & Curry (1992)]
233
234         if (zflwp(i).eq.0.) rel = 1.
235         if (zfiwp(i).eq.0. .or. rei.le.0.) rei = 1.
236         pcltau(i,k) = 3.0/2.0 * ( zflwp(i)/rel )
237     .             + zfiwp(i) * (3.448e-03  + 2.431/rei)
238
239c -- cloud infrared emissivity:
240
241c [the broadband infrared absorption coefficient is parameterized
242c  as a function of the effective cld droplet radius]
243
244c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
245         k_ice = k_ice0 + 1.0/rei
246
247         pclemi(i,k) = 1.0
248     .      - EXP( - coef_chau*zflwp(i) - DF*k_ice*zfiwp(i) )
249
250         endif ! ok_newmicro
251
252         lo = (pclc(i,k) .LE. seuil_neb)
253         IF (lo) pclc(i,k) = 0.0
254         IF (lo) pcltau(i,k) = 0.0
255         IF (lo) pclemi(i,k) = 0.0
256         
257         IF (lo) cldtaupi(i,k) = 0.0
258         IF (.NOT.ok_aie) cldtaupi(i,k)=pcltau(i,k)           
259      ENDDO
260      ENDDO
261ccc      DO k = 1, klev
262ccc      DO i = 1, klon
263ccc         t(i,k) = t(i,k)
264ccc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
265ccc         lo = pclc(i,k) .GT. (2.*1.e-5)
266ccc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
267ccc     .          /(rg*pclc(i,k))
268ccc         zradef = 10.0 + (1.-sigs(k))*45.0
269ccc         pcltau(i,k) = 1.5 * zflwp / zradef
270ccc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
271ccc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
272ccc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
273ccc         if (.NOT.lo) pclc(i,k) = 0.0
274ccc         if (.NOT.lo) pcltau(i,k) = 0.0
275ccc         if (.NOT.lo) pclemi(i,k) = 0.0
276ccc      ENDDO
277ccc      ENDDO
278cccccc      print*, 'pas de nuage dans le rayonnement'
279cccccc      DO k = 1, klev
280cccccc      DO i = 1, klon
281cccccc         pclc(i,k) = 0.0
282cccccc         pcltau(i,k) = 0.0
283cccccc         pclemi(i,k) = 0.0
284cccccc      ENDDO
285cccccc      ENDDO
286C
287C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
288C
289      DO i = 1, klon
290         pct(i)=1.0
291         pch(i)=1.0
292         pcm(i) = 1.0
293         pcl(i) = 1.0
294         pctlwp(i) = 0.0
295      ENDDO
296C
297      DO k = klev, 1, -1
298      DO i = 1, klon
299         pctlwp(i) = pctlwp(i)
300     .             + pqlwp(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
301         pct(i) = pct(i)*(1.0-pclc(i,k))
302         if (pplay(i,k).LE.cetahb*paprs(i,1))
303     .      pch(i) = pch(i)*(1.0-pclc(i,k))
304         if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
305     .       pplay(i,k).LE.cetamb*paprs(i,1))
306     .      pcm(i) = pcm(i)*(1.0-pclc(i,k))
307         if (pplay(i,k).GT.cetamb*paprs(i,1))
308     .      pcl(i) = pcl(i)*(1.0-pclc(i,k))
309      ENDDO
310      ENDDO
311C
312      DO i = 1, klon
313         pct(i)=1.-pct(i)
314         pch(i)=1.-pch(i)
315         pcm(i)=1.-pcm(i)
316         pcl(i)=1.-pcl(i)
317      ENDDO
318C
319      RETURN
320      END
Note: See TracBrowser for help on using the repository browser.