source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/newmicro.F @ 5192

Last change on this file since 5192 was 634, checked in by Laurent Fairhead, 20 years ago

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