source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F @ 1105

Last change on this file since 1105 was 1092, checked in by yann meurdesoif, 15 years ago

Optimisation Othman Bouizi : newmicro

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