source: LMDZ.3.3/branches/rel-LF/libf/phylmd/newmicro.F @ 5456

Last change on this file since 5456 was 535, checked in by lmdzadmin, 21 years ago

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