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

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