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

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

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.2 KB
Line 
1! $Id: newmicro.F 1279 2009-12-10 09:02:56Z 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, 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 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
145               zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
146               zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
147c     IM Total Liquid/Ice water content                                   
148               xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
149               xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
150c     IM 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            ENDDO
154         ENDDO
155
156         IF (ok_aie) THEN
157            DO k = 1, klev
158               DO i = 1, klon
159                                ! Formula "D" of Boucher and Lohmann, Tellus, 1995
160                                !             
161                  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
162     &               log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
163                                ! Cloud droplet number concentration (CDNC) is restricted
164                                ! to be within [20, 1000 cm^3]
165                                !
166                  cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
167                                !
168                                !
169                  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
170     &               log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
171                  cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
172               ENDDO
173            ENDDO
174            DO k = 1, klev
175               DO i = 1, klon
176!                  rad_chaud_tab(i,k) =
177!     &                 MAX(1.1e6
178!     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
179!     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.)
180                  rad_chaud_tab(i,k) =
181     &                 1.1
182     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
183     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
184                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.)
185               ENDDO           
186            ENDDO
187         ELSE
188            DO k = 1, MIN(3,klev)
189               DO i = 1, klon
190                  rad_chaud_tab(i,k) = rad_chau2
191               ENDDO           
192            ENDDO
193            DO k = MIN(3,klev)+1, klev
194               DO i = 1, klon
195                  rad_chaud_tab(i,k) = rad_chau1
196               ENDDO           
197            ENDDO
198
199         ENDIF
200         
201         DO k = 1, klev
202!            IF(.not.ok_aie) THEN
203            rad_chaud = rad_chau1
204            IF (k.LE.3) rad_chaud = rad_chau2
205!            ENDIF
206            DO i = 1, klon
207               IF (pclc(i,k) .LE. seuil_neb) THEN
208               
209c     -- effective cloud droplet radius (microns):
210               
211c     for liquid water clouds:
212                                ! For output diagnostics
213                                !
214                                ! Cloud droplet effective radius [um]
215                                !
216                                ! we multiply here with f * xl (fraction of liquid water
217                                ! clouds in the grid cell) to avoid problems in the
218                                ! averaging of the output.
219                                ! In the output of IOIPSL, derive the real cloud droplet
220                                ! effective radius as re/fl
221                                !
222                                   
223                  fl(i,k) = seuil_neb*(1.-zfice2(i,k))           
224                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
225                 
226                  rel = 0.
227                  rei = 0.
228                  pclc(i,k) = 0.0
229                  pcltau(i,k) = 0.0
230                  pclemi(i,k) = 0.0
231                  cldtaupi(i,k) = 0.0                 
232               ELSE
233
234c     -- liquid/ice cloud water paths:
235                 
236                  zflwp_var= 1000.*(1.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k)
237     &                 *diff_paprs(i,k)
238                  zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k)
239     &                 *diff_paprs(i,k)
240                 
241c     -- effective cloud droplet radius (microns):
242               
243c     for liquid water clouds:
244                                   
245                  IF (ok_aie) THEN
246                     radius =
247     &                    1.1
248     &                    *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
249     &                    /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.)
250                     radius = MAX(radius*1e6, 5.)
251                 
252                     tc = t(i,k)-273.15
253                     rei = 0.71*tc + 61.29
254                     if (tc.le.-81.4) rei = 3.5
255                     if (zflwp_var.eq.0.) radius = 1.
256                     if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
257                     cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
258     &                    + zfiwp_var * (3.448e-03  + 2.431/rei)
259
260                  ENDIF         ! ok_aie
261                                ! For output diagnostics
262                                !
263                                ! Cloud droplet effective radius [um]
264                                !
265                                ! we multiply here with f * xl (fraction of liquid water
266                                ! clouds in the grid cell) to avoid problems in the
267                                ! averaging of the output.
268                                ! In the output of IOIPSL, derive the real cloud droplet
269                                ! effective radius as re/fl
270                                !
271 
272                  fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))           
273                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
274                 
275                  rel = rad_chaud_tab(i,k)
276c     for ice clouds: as a function of the ambiant temperature
277c     [formula used by Iacobellis and Somerville (2000), with an
278c     asymptotical value of 3.5 microns at T<-81.4 C added to be
279c     consistent with observations of Heymsfield et al. 1986]:
280                  tc = t(i,k)-273.15
281                  rei = 0.71*tc + 61.29
282                  if (tc.le.-81.4) rei = 3.5
283c     -- cloud optical thickness :
284               
285c     [for liquid clouds, traditional formula,
286c     for ice clouds, Ebert & Curry (1992)]
287                 
288                 if (zflwp_var.eq.0.) rel = 1.
289                 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
290                 pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
291     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
292c     -- cloud infrared emissivity:
293               
294c     [the broadband infrared absorption coefficient is parameterized
295c     as a function of the effective cld droplet radius]
296               
297c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
298                  k_ice = k_ice0 + 1.0/rei
299                 
300                  pclemi(i,k) = 1.0
301     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
302
303               ENDIF
304              reliq(i,k)=rel
305              reice(i,k)=rei
306!              if (i.eq.1) then
307!              print*,'Dans newmicro rel, rei :',rel, rei
308!              print*,'Dans newmicro reliq, reice :',
309!     $             reliq(i,k),reice(i,k)
310!              endif
311
312            ENDDO
313         ENDDO
314
315         DO k = 1, klev
316            DO i = 1, klon
317               xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
318               xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
319            ENDDO
320         ENDDO
321
322      ELSE
323         DO k = 1, klev
324            rad_chaud = rad_chau1
325            IF (k.LE.3) rad_chaud = rad_chau2
326            DO i = 1, klon
327                             
328               IF (pclc(i,k) .LE. seuil_neb) THEN
329
330                  pclc(i,k) = 0.0
331                  pcltau(i,k) = 0.0
332                  pclemi(i,k) = 0.0
333                  cldtaupi(i,k) = 0.0
334
335               ELSE
336
337                  zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k)
338     &                 /pclc(i,k)
339                 
340                  zfice1 = MIN(
341     &                 MAX( 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
342     &                 ,0.0),1.0)**nexpo
343                 
344                  radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
345                  coef   = coef_chau * (1.-zfice1) + coef_froi * zfice1
346
347                  pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
348                  pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
349
350               ENDIF
351                             
352            ENDDO
353         ENDDO
354      ENDIF
355     
356      IF (.NOT.ok_aie) THEN
357         DO k = 1, klev
358            DO i = 1, klon
359               cldtaupi(i,k)=pcltau(i,k)
360            ENDDO
361         ENDDO               
362      ENDIF
363
364ccc   DO k = 1, klev
365ccc   DO i = 1, klon
366ccc   t(i,k) = t(i,k)
367ccc   pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
368ccc   lo = pclc(i,k) .GT. (2.*1.e-5)
369ccc   zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
370ccc   .          /(rg*pclc(i,k))
371ccc   zradef = 10.0 + (1.-sigs(k))*45.0
372ccc   pcltau(i,k) = 1.5 * zflwp / zradef
373ccc   zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
374ccc   zmsac = 0.13*(1.0-zfice) + 0.08*zfice
375ccc   pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
376ccc   if (.NOT.lo) pclc(i,k) = 0.0
377ccc   if (.NOT.lo) pcltau(i,k) = 0.0
378ccc   if (.NOT.lo) pclemi(i,k) = 0.0
379ccc   ENDDO
380ccc   ENDDO
381ccccc print*, 'pas de nuage dans le rayonnement'
382ccccc DO k = 1, klev
383ccccc DO i = 1, klon
384ccccc pclc(i,k) = 0.0
385ccccc pcltau(i,k) = 0.0
386ccccc pclemi(i,k) = 0.0
387ccccc ENDDO
388ccccc ENDDO
389C     
390C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
391C     
392c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
393c     initialisations
394      DO i=1,klon
395         zclear(i)=1.
396         zcloud(i)=0.
397         pch(i)=1.0
398         pcm(i) = 1.0
399         pcl(i) = 1.0
400         pctlwp(i) = 0.0
401      ENDDO
402C
403cIM cf CR DO k=1,klev
404      DO k = klev, 1, -1
405         DO i = 1, klon
406            pctlwp(i) = pctlwp(i)
407     &           + pqlwp(i,k)*diff_paprs(i,k)
408         ENDDO
409      ENDDO
410c     IM cf. CR
411      IF (NOVLP.EQ.1) THEN
412         DO k = klev, 1, -1
413            DO i = 1, klon
414               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
415     &              /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
416               pct(i)=1.-zclear(i)
417               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
418                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
419     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
420               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
421     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
422                  pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
423     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
424               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
425                  pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
426     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
427               endif
428               zcloud(i)=pclc(i,k)
429            ENDDO
430         ENDDO
431      ELSE IF (NOVLP.EQ.2) THEN
432         DO k = klev, 1, -1
433            DO i = 1, klon
434               zcloud(i)=MAX(pclc(i,k),zcloud(i))
435               pct(i)=zcloud(i)
436               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
437                  pch(i) = MIN(pclc(i,k),pch(i))
438               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
439     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
440                  pcm(i) = MIN(pclc(i,k),pcm(i))
441               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
442                  pcl(i) = MIN(pclc(i,k),pcl(i))
443               endif
444            ENDDO
445         ENDDO
446      ELSE IF (NOVLP.EQ.3) THEN
447         DO k = klev, 1, -1
448            DO i = 1, klon
449               zclear(i)=zclear(i)*(1.-pclc(i,k))
450               pct(i)=1-zclear(i)
451               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
452                  pch(i) = pch(i)*(1.0-pclc(i,k))
453               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
454     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
455                  pcm(i) = pcm(i)*(1.0-pclc(i,k))
456               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
457                  pcl(i) = pcl(i)*(1.0-pclc(i,k))
458               endif
459            ENDDO
460         ENDDO
461      ENDIF
462     
463C     
464      DO i = 1, klon
465c     IM cf. CR          pct(i)=1.-pct(i)
466         pch(i)=1.-pch(i)
467         pcm(i)=1.-pcm(i)
468         pcl(i)=1.-pcl(i)
469      ENDDO
470     
471C
472      RETURN
473      END
Note: See TracBrowser for help on using the repository browser.