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

Last change on this file since 1286 was 1286, checked in by Laurent Fairhead, 14 years ago

Les parametres definissant la transition eau glace dans les nuages sont
maintenant lus dans physiq.def. Les valeurs par defaut donnent les memes
resultats que precedemment. JLD
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.3 KB
Line 
1! $Id: newmicro.F 1286 2009-12-17 13:47:10Z fairhead $
2!     
3      SUBROUTINE newmicro (paprs, pplay,ok_newmicro,
4     .                  t, pqlwp, pclc, pcltau, pclemi,
5     .                  pch, pcl, pcm, pct, pctlwp,
6     s                  xflwp, xfiwp, xflwc, xfiwc,
7     e                  ok_aie,
8     e                  mass_solu_aero, mass_solu_aero_pi,
9     e                  bl95_b0, bl95_b1,
10     s                  cldtaupi, re, fl, reliq, reice)
11
12      USE dimphy
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 mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
25c mass_solu_aero_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
42cym#include "dimensions.h"
43cym#include "dimphy.h"
44#include "nuage.h"
45cIM cf. CR: include pour NOVLP et ZEPSEC
46#include "radepsi.h"
47#include "radopt.h"
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
76      PARAMETER (seuil_neb=0.001)
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 mass_solu_aero(klon, klev)    ! total mass concentration for all soluble aerosols [ug m-3]
97      REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value)
98      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
99      REAL re(klon, klev)       ! cloud droplet effective radius [um]
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   
109cIM cf. CR:parametres supplementaires
110      REAL zclear(klon)
111      REAL zcloud(klon)
112
113c **************************
114c *                        *
115c * DEBUT PARTIE OPTIMISEE *
116c *                        *
117c **************************
118
119      REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev)
120      REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var
121
122! Abderrahmane oct 2009
123      Real reliq(klon, klev), reice(klon, klev)
124
125c
126c Calculer l'epaisseur optique et l'emmissivite des nuages
127c
128c     IM inversion des DO
129      xflwp = 0.d0
130      xfiwp = 0.d0
131      xflwc = 0.d0
132      xfiwc = 0.d0
133
134      DO k = 1, klev
135         DO i = 1, klon
136            diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
137         ENDDO
138      ENDDO
139
140      IF (ok_newmicro) THEN
141
142
143         DO k = 1, klev
144            DO i = 1, klon
145c               zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
146               zfice2(i,k) = 1.0 - (t(i,k)-t_glace_min) /
147     &                             (t_glace_max-t_glace_min)
148               zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
149c     IM Total Liquid/Ice water content                                   
150               xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
151               xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
152c     IM 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            ENDDO
156         ENDDO
157
158         IF (ok_aie) THEN
159            DO k = 1, klev
160               DO i = 1, klon
161                                ! Formula "D" of Boucher and Lohmann, Tellus, 1995
162                                !             
163                  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
164     &               log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
165                                ! Cloud droplet number concentration (CDNC) is restricted
166                                ! to be within [20, 1000 cm^3]
167                                !
168                  cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
169                                !
170                                !
171                  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
172     &               log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
173                  cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
174               ENDDO
175            ENDDO
176            DO k = 1, klev
177               DO i = 1, klon
178!                  rad_chaud_tab(i,k) =
179!     &                 MAX(1.1e6
180!     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
181!     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.)
182                  rad_chaud_tab(i,k) =
183     &                 1.1
184     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
185     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
186                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.)
187               ENDDO           
188            ENDDO
189         ELSE
190            DO k = 1, MIN(3,klev)
191               DO i = 1, klon
192                  rad_chaud_tab(i,k) = rad_chau2
193               ENDDO           
194            ENDDO
195            DO k = MIN(3,klev)+1, klev
196               DO i = 1, klon
197                  rad_chaud_tab(i,k) = rad_chau1
198               ENDDO           
199            ENDDO
200
201         ENDIF
202         
203         DO k = 1, klev
204!            IF(.not.ok_aie) THEN
205            rad_chaud = rad_chau1
206            IF (k.LE.3) rad_chaud = rad_chau2
207!            ENDIF
208            DO i = 1, klon
209               IF (pclc(i,k) .LE. seuil_neb) THEN
210               
211c     -- effective cloud droplet radius (microns):
212               
213c     for liquid water clouds:
214                                ! For output diagnostics
215                                !
216                                ! Cloud droplet effective radius [um]
217                                !
218                                ! we multiply here with f * xl (fraction of liquid water
219                                ! clouds in the grid cell) to avoid problems in the
220                                ! averaging of the output.
221                                ! In the output of IOIPSL, derive the real cloud droplet
222                                ! effective radius as re/fl
223                                !
224                                   
225                  fl(i,k) = seuil_neb*(1.-zfice2(i,k))           
226                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
227                 
228                  rel = 0.
229                  rei = 0.
230                  pclc(i,k) = 0.0
231                  pcltau(i,k) = 0.0
232                  pclemi(i,k) = 0.0
233                  cldtaupi(i,k) = 0.0                 
234               ELSE
235
236c     -- liquid/ice cloud water paths:
237                 
238                  zflwp_var= 1000.*(1.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k)
239     &                 *diff_paprs(i,k)
240                  zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k)
241     &                 *diff_paprs(i,k)
242                 
243c     -- effective cloud droplet radius (microns):
244               
245c     for liquid water clouds:
246                                   
247                  IF (ok_aie) THEN
248                     radius =
249     &                    1.1
250     &                    *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
251     &                    /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.)
252                     radius = MAX(radius*1e6, 5.)
253                 
254                     tc = t(i,k)-273.15
255                     rei = 0.71*tc + 61.29
256                     if (tc.le.-81.4) rei = 3.5
257                     if (zflwp_var.eq.0.) radius = 1.
258                     if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
259                     cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
260     &                    + zfiwp_var * (3.448e-03  + 2.431/rei)
261
262                  ENDIF         ! ok_aie
263                                ! For output diagnostics
264                                !
265                                ! Cloud droplet effective radius [um]
266                                !
267                                ! we multiply here with f * xl (fraction of liquid water
268                                ! clouds in the grid cell) to avoid problems in the
269                                ! averaging of the output.
270                                ! In the output of IOIPSL, derive the real cloud droplet
271                                ! effective radius as re/fl
272                                !
273 
274                  fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))           
275                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
276                 
277                  rel = rad_chaud_tab(i,k)
278c     for ice clouds: as a function of the ambiant temperature
279c     [formula used by Iacobellis and Somerville (2000), with an
280c     asymptotical value of 3.5 microns at T<-81.4 C added to be
281c     consistent with observations of Heymsfield et al. 1986]:
282                  tc = t(i,k)-273.15
283                  rei = 0.71*tc + 61.29
284                  if (tc.le.-81.4) rei = 3.5
285c     -- cloud optical thickness :
286               
287c     [for liquid clouds, traditional formula,
288c     for ice clouds, Ebert & Curry (1992)]
289                 
290                 if (zflwp_var.eq.0.) rel = 1.
291                 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
292                 pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
293     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
294c     -- cloud infrared emissivity:
295               
296c     [the broadband infrared absorption coefficient is parameterized
297c     as a function of the effective cld droplet radius]
298               
299c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
300                  k_ice = k_ice0 + 1.0/rei
301                 
302                  pclemi(i,k) = 1.0
303     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
304
305               ENDIF
306              reliq(i,k)=rel
307              reice(i,k)=rei
308!              if (i.eq.1) then
309!              print*,'Dans newmicro rel, rei :',rel, rei
310!              print*,'Dans newmicro reliq, reice :',
311!     $             reliq(i,k),reice(i,k)
312!              endif
313
314            ENDDO
315         ENDDO
316
317         DO k = 1, klev
318            DO i = 1, klon
319               xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
320               xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
321            ENDDO
322         ENDDO
323
324      ELSE
325         DO k = 1, klev
326            rad_chaud = rad_chau1
327            IF (k.LE.3) rad_chaud = rad_chau2
328            DO i = 1, klon
329                             
330               IF (pclc(i,k) .LE. seuil_neb) THEN
331
332                  pclc(i,k) = 0.0
333                  pcltau(i,k) = 0.0
334                  pclemi(i,k) = 0.0
335                  cldtaupi(i,k) = 0.0
336
337               ELSE
338
339                  zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k)
340     &                 /pclc(i,k)
341                 
342                  zfice1 = MIN(
343     &                 MAX( 1.0 - (t(i,k)-t_glace_min) /
344     &                    (t_glace_max-t_glace_min),0.0),1.0)**nexpo
345                 
346                  radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
347                  coef   = coef_chau * (1.-zfice1) + coef_froi * zfice1
348
349                  pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
350                  pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
351
352               ENDIF
353                             
354            ENDDO
355         ENDDO
356      ENDIF
357     
358      IF (.NOT.ok_aie) THEN
359         DO k = 1, klev
360            DO i = 1, klon
361               cldtaupi(i,k)=pcltau(i,k)
362            ENDDO
363         ENDDO               
364      ENDIF
365
366ccc   DO k = 1, klev
367ccc   DO i = 1, klon
368ccc   t(i,k) = t(i,k)
369ccc   pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
370ccc   lo = pclc(i,k) .GT. (2.*1.e-5)
371ccc   zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
372ccc   .          /(rg*pclc(i,k))
373ccc   zradef = 10.0 + (1.-sigs(k))*45.0
374ccc   pcltau(i,k) = 1.5 * zflwp / zradef
375ccc   zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
376ccc   zmsac = 0.13*(1.0-zfice) + 0.08*zfice
377ccc   pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
378ccc   if (.NOT.lo) pclc(i,k) = 0.0
379ccc   if (.NOT.lo) pcltau(i,k) = 0.0
380ccc   if (.NOT.lo) pclemi(i,k) = 0.0
381ccc   ENDDO
382ccc   ENDDO
383ccccc print*, 'pas de nuage dans le rayonnement'
384ccccc DO k = 1, klev
385ccccc DO i = 1, klon
386ccccc pclc(i,k) = 0.0
387ccccc pcltau(i,k) = 0.0
388ccccc pclemi(i,k) = 0.0
389ccccc ENDDO
390ccccc ENDDO
391C     
392C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
393C     
394c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
395c     initialisations
396      DO i=1,klon
397         zclear(i)=1.
398         zcloud(i)=0.
399         pch(i)=1.0
400         pcm(i) = 1.0
401         pcl(i) = 1.0
402         pctlwp(i) = 0.0
403      ENDDO
404C
405cIM cf CR DO k=1,klev
406      DO k = klev, 1, -1
407         DO i = 1, klon
408            pctlwp(i) = pctlwp(i)
409     &           + pqlwp(i,k)*diff_paprs(i,k)
410         ENDDO
411      ENDDO
412c     IM cf. CR
413      IF (NOVLP.EQ.1) THEN
414         DO k = klev, 1, -1
415            DO i = 1, klon
416               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
417     &              /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
418               pct(i)=1.-zclear(i)
419               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
420                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
421     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
422               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
423     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
424                  pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
425     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
426               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
427                  pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
428     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
429               endif
430               zcloud(i)=pclc(i,k)
431            ENDDO
432         ENDDO
433      ELSE IF (NOVLP.EQ.2) THEN
434         DO k = klev, 1, -1
435            DO i = 1, klon
436               zcloud(i)=MAX(pclc(i,k),zcloud(i))
437               pct(i)=zcloud(i)
438               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
439                  pch(i) = MIN(pclc(i,k),pch(i))
440               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
441     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
442                  pcm(i) = MIN(pclc(i,k),pcm(i))
443               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
444                  pcl(i) = MIN(pclc(i,k),pcl(i))
445               endif
446            ENDDO
447         ENDDO
448      ELSE IF (NOVLP.EQ.3) THEN
449         DO k = klev, 1, -1
450            DO i = 1, klon
451               zclear(i)=zclear(i)*(1.-pclc(i,k))
452               pct(i)=1-zclear(i)
453               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
454                  pch(i) = pch(i)*(1.0-pclc(i,k))
455               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
456     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
457                  pcm(i) = pcm(i)*(1.0-pclc(i,k))
458               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
459                  pcl(i) = pcl(i)*(1.0-pclc(i,k))
460               endif
461            ENDDO
462         ENDDO
463      ENDIF
464     
465C     
466      DO i = 1, klon
467c     IM cf. CR          pct(i)=1.-pct(i)
468         pch(i)=1.-pch(i)
469         pcm(i)=1.-pcm(i)
470         pcl(i)=1.-pcl(i)
471      ENDDO
472     
473C
474      RETURN
475      END
Note: See TracBrowser for help on using the repository browser.