source: LMDZ4/branches/V3_test/libf/phylmd/newmicro.F @ 789

Last change on this file since 789 was 747, checked in by Laurent Fairhead, 18 years ago

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