source: LMDZ4/branches/LMDZ4V5.0-LF/libf/phylmd/newmicro.F @ 2124

Last change on this file since 2124 was 1306, checked in by Laurent Fairhead, 15 years ago

Too many characters on line and at the wrong place

Ligne trop longue au mauvais endroit

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