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

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

Additions to aerosol outputs for CMIP5 exercise

Additions aux sorties aérosols pour l'exercice CMIP5

Michael, Anne
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.3 KB
Line 
1! $Id: newmicro.F 1337 2010-04-02 11:31:05Z 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      USE phys_local_var_mod, only: scdnc,cldncl,reffclwtop,lcc,
14     .                              reffclws,reffclwc,cldnvi,lcc3d,
15     .                              lcc3dcon,lcc3dstra
16      USE phys_state_var_mod, only: rnebcon,clwcon
17      IMPLICIT none
18c======================================================================
19c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
20c Objet: Calculer epaisseur optique et emmissivite des nuages
21c======================================================================
22c Arguments:
23c t-------input-R-temperature
24c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
25c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
26c
27c ok_aie--input-L-apply aerosol indirect effect or not
28c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
29c mass_solu_aero_pi--input-R-dito, pre-industrial value
30c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
31c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
32c     
33c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
34c                   needed for the diagnostics of the aerosol indirect
35c                   radiative forcing (see radlwsw)
36c re------output-R-Cloud droplet effective radius multiplied by fl [um]
37c fl------output-R-Denominator to re, introduced to avoid problems in
38c                  the averaging of the output. fl is the fraction of liquid
39c                  water clouds within a grid cell           
40c pcltau--output-R-epaisseur optique des nuages
41c pclemi--output-R-emissivite des nuages (0 a 1)
42c======================================================================
43C
44#include "YOMCST.h"
45c
46cym#include "dimensions.h"
47cym#include "dimphy.h"
48#include "nuage.h"
49cIM cf. CR: include pour NOVLP et ZEPSEC
50#include "radepsi.h"
51#include "radopt.h"
52c choix de l'hypothese de recouvrememnt nuageuse
53      LOGICAL RANDOM,MAXIMUM_RANDOM,MAXIMUM,FIRST
54      parameter (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.)
55c Hypoyhese de recouvrement : MAXIMUM_RANDOM
56      INTEGER flag_max
57      REAL phase3d(klon, klev),dh(klon, klev),pdel(klon, klev),
58     .     zrho(klon, klev)
59      REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
60      REAL thres_tau,thres_neb
61      PARAMETER (thres_tau=0.3, thres_neb=0.001)
62      REAL t_tmp
63      REAL gravit
64      PARAMETER (gravit=9.80616)  !m/s2
65      REAL pqlwpcon(klon, klev), pqlwpstra(klon, klev)
66c
67      REAL paprs(klon,klev+1), pplay(klon,klev)
68      REAL t(klon,klev)
69c
70      REAL pclc(klon,klev)
71      REAL pqlwp(klon,klev)
72      REAL pcltau(klon,klev), pclemi(klon,klev)
73c
74      REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
75c
76      LOGICAL lo
77c
78      REAL cetahb, cetamb
79      PARAMETER (cetahb = 0.45, cetamb = 0.80)
80C
81      INTEGER i, k
82cIM: 091003   REAL zflwp, zradef, zfice, zmsac
83      REAL zflwp(klon), zradef, zfice, zmsac
84cIM: 091003 rajout
85      REAL xflwp(klon), xfiwp(klon)
86      REAL xflwc(klon,klev), xfiwc(klon,klev)
87c
88      REAL radius, rad_chaud
89cc      PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
90ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
91c sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
92      REAL coef, coef_froi, coef_chau
93      PARAMETER (coef_chau=0.13, coef_froi=0.09)
94      REAL seuil_neb
95      PARAMETER (seuil_neb=0.001)
96      INTEGER nexpo ! exponentiel pour glace/eau
97      PARAMETER (nexpo=6)
98ccc      PARAMETER (nexpo=1)
99
100c -- sb:
101      logical ok_newmicro
102c     parameter (ok_newmicro=.FALSE.)
103cIM: 091003   real rel, tc, rei, zfiwp
104      real rel, tc, rei, zfiwp(klon)
105      real k_liq, k_ice0, k_ice, DF
106      parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g
107      parameter (DF=1.66) ! diffusivity factor
108c sb --
109cjq for the aerosol indirect effect
110cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
111cjq     
112      LOGICAL ok_aie            ! Apply AIE or not?
113      LOGICAL ok_a1lwpdep       ! a1 LWP dependent?
114     
115      REAL mass_solu_aero(klon, klev)    ! total mass concentration for all soluble aerosols [ug m-3]
116      REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value)
117      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
118      REAL re(klon, klev)       ! cloud droplet effective radius [um]
119      REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
120      REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
121     
122      REAL fl(klon, klev)       ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
123     
124      REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
125     
126      REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
127cjq-end   
128cIM cf. CR:parametres supplementaires
129      REAL zclear(klon)
130      REAL zcloud(klon)
131
132c **************************
133c *                        *
134c * DEBUT PARTIE OPTIMISEE *
135c *                        *
136c **************************
137
138      REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev)
139      REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var
140
141! Abderrahmane oct 2009
142      Real reliq(klon, klev), reice(klon, klev)
143
144c
145c Calculer l'epaisseur optique et l'emmissivite des nuages
146c
147c     IM inversion des DO
148      xflwp = 0.d0
149      xfiwp = 0.d0
150      xflwc = 0.d0
151      xfiwc = 0.d0
152
153! Initialisation
154      reliq=0.
155      reice=0.
156
157      DO k = 1, klev
158         DO i = 1, klon
159            diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
160         ENDDO
161      ENDDO
162
163      IF (ok_newmicro) THEN
164
165
166         DO k = 1, klev
167            DO i = 1, klon
168c               zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
169               zfice2(i,k) = 1.0 - (t(i,k)-t_glace_min) /
170     &                             (t_glace_max-t_glace_min)
171               zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
172c     IM Total Liquid/Ice water content                                   
173               xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
174               xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
175c     IM In-Cloud Liquid/Ice water content
176c     xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
177c     xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
178            ENDDO
179         ENDDO
180
181         IF (ok_aie) THEN
182            DO k = 1, klev
183               DO i = 1, klon
184                                ! Formula "D" of Boucher and Lohmann, Tellus, 1995
185                                !             
186                  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
187     &               log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
188                                ! Cloud droplet number concentration (CDNC) is restricted
189                                ! to be within [20, 1000 cm^3]
190                                !
191                  cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
192                                !
193                                !
194                  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
195     &               log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))
196     &               *1.e6 !-m-3
197                  cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
198               ENDDO
199            ENDDO
200            DO k = 1, klev
201               DO i = 1, klon
202!                  rad_chaud_tab(i,k) =
203!     &                 MAX(1.1e6
204!     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
205!     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.)
206                  rad_chaud_tab(i,k) =
207     &                 1.1
208     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
209     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
210                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.)
211               ENDDO           
212            ENDDO
213         ELSE
214            DO k = 1, MIN(3,klev)
215               DO i = 1, klon
216                  rad_chaud_tab(i,k) = rad_chau2
217               ENDDO           
218            ENDDO
219            DO k = MIN(3,klev)+1, klev
220               DO i = 1, klon
221                  rad_chaud_tab(i,k) = rad_chau1
222               ENDDO           
223            ENDDO
224
225         ENDIF
226         
227         DO k = 1, klev
228!            IF(.not.ok_aie) THEN
229            rad_chaud = rad_chau1
230            IF (k.LE.3) rad_chaud = rad_chau2
231!            ENDIF
232            DO i = 1, klon
233               IF (pclc(i,k) .LE. seuil_neb) THEN
234               
235c     -- effective cloud droplet radius (microns):
236               
237c     for liquid water clouds:
238                                ! For output diagnostics
239                                !
240                                ! Cloud droplet effective radius [um]
241                                !
242                                ! we multiply here with f * xl (fraction of liquid water
243                                ! clouds in the grid cell) to avoid problems in the
244                                ! averaging of the output.
245                                ! In the output of IOIPSL, derive the real cloud droplet
246                                ! effective radius as re/fl
247                                !
248                                   
249                  fl(i,k) = seuil_neb*(1.-zfice2(i,k))           
250                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
251                 
252                  rel = 0.
253                  rei = 0.
254                  pclc(i,k) = 0.0
255                  pcltau(i,k) = 0.0
256                  pclemi(i,k) = 0.0
257                  cldtaupi(i,k) = 0.0                 
258               ELSE
259
260c     -- liquid/ice cloud water paths:
261                 
262                  zflwp_var= 1000.*(1.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k)
263     &                 *diff_paprs(i,k)
264                  zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k)
265     &                 *diff_paprs(i,k)
266                 
267c     -- effective cloud droplet radius (microns):
268               
269c     for liquid water clouds:
270                                   
271                  IF (ok_aie) THEN
272                     radius =
273     &                    1.1
274     &                    *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
275     &                    /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.)
276                     radius = MAX(radius*1e6, 5.)
277                 
278                     tc = t(i,k)-273.15
279                     rei = 0.71*tc + 61.29
280                     if (tc.le.-81.4) rei = 3.5
281                     if (zflwp_var.eq.0.) radius = 1.
282                     if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
283                     cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
284     &                    + zfiwp_var * (3.448e-03  + 2.431/rei)
285
286                  ENDIF         ! ok_aie
287                                ! For output diagnostics
288                                !
289                                ! Cloud droplet effective radius [um]
290                                !
291                                ! we multiply here with f * xl (fraction of liquid water
292                                ! clouds in the grid cell) to avoid problems in the
293                                ! averaging of the output.
294                                ! In the output of IOIPSL, derive the real cloud droplet
295                                ! effective radius as re/fl
296                                !
297 
298                  fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))           
299                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
300                 
301                  rel = rad_chaud_tab(i,k)
302c     for ice clouds: as a function of the ambiant temperature
303c     [formula used by Iacobellis and Somerville (2000), with an
304c     asymptotical value of 3.5 microns at T<-81.4 C added to be
305c     consistent with observations of Heymsfield et al. 1986]:
306                  tc = t(i,k)-273.15
307                  rei = 0.71*tc + 61.29
308                  if (tc.le.-81.4) rei = 3.5
309c     -- cloud optical thickness :
310               
311c     [for liquid clouds, traditional formula,
312c     for ice clouds, Ebert & Curry (1992)]
313                 
314                 if (zflwp_var.eq.0.) rel = 1.
315                 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
316                 pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
317     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
318c     -- cloud infrared emissivity:
319               
320c     [the broadband infrared absorption coefficient is parameterized
321c     as a function of the effective cld droplet radius]
322               
323c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
324                  k_ice = k_ice0 + 1.0/rei
325                 
326                  pclemi(i,k) = 1.0
327     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
328
329               ENDIF
330              reliq(i,k)=rel
331              reice(i,k)=rei
332!              if (i.eq.1) then
333!              print*,'Dans newmicro rel, rei :',rel, rei
334!              print*,'Dans newmicro reliq, reice :',
335!     $             reliq(i,k),reice(i,k)
336!              endif
337
338            ENDDO
339         ENDDO
340
341         DO k = 1, klev
342            DO i = 1, klon
343               xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
344               xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
345            ENDDO
346         ENDDO
347
348      ELSE
349         DO k = 1, klev
350            rad_chaud = rad_chau1
351            IF (k.LE.3) rad_chaud = rad_chau2
352            DO i = 1, klon
353                             
354               IF (pclc(i,k) .LE. seuil_neb) THEN
355
356                  pclc(i,k) = 0.0
357                  pcltau(i,k) = 0.0
358                  pclemi(i,k) = 0.0
359                  cldtaupi(i,k) = 0.0
360
361               ELSE
362
363                  zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k)
364     &                 /pclc(i,k)
365                 
366                  zfice1 = MIN(
367     &                 MAX( 1.0 - (t(i,k)-t_glace_min) /
368     &                    (t_glace_max-t_glace_min),0.0),1.0)**nexpo
369                 
370                  radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
371                  coef   = coef_chau * (1.-zfice1) + coef_froi * zfice1
372
373                  pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
374                  pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
375
376               ENDIF
377                             
378            ENDDO
379         ENDDO
380      ENDIF
381     
382      IF (.NOT.ok_aie) THEN
383         DO k = 1, klev
384            DO i = 1, klon
385               cldtaupi(i,k)=pcltau(i,k)
386            ENDDO
387         ENDDO               
388      ENDIF
389
390ccc   DO k = 1, klev
391ccc   DO i = 1, klon
392ccc   t(i,k) = t(i,k)
393ccc   pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
394ccc   lo = pclc(i,k) .GT. (2.*1.e-5)
395ccc   zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
396ccc   .          /(rg*pclc(i,k))
397ccc   zradef = 10.0 + (1.-sigs(k))*45.0
398ccc   pcltau(i,k) = 1.5 * zflwp / zradef
399ccc   zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
400ccc   zmsac = 0.13*(1.0-zfice) + 0.08*zfice
401ccc   pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
402ccc   if (.NOT.lo) pclc(i,k) = 0.0
403ccc   if (.NOT.lo) pcltau(i,k) = 0.0
404ccc   if (.NOT.lo) pclemi(i,k) = 0.0
405ccc   ENDDO
406ccc   ENDDO
407ccccc print*, 'pas de nuage dans le rayonnement'
408ccccc DO k = 1, klev
409ccccc DO i = 1, klon
410ccccc pclc(i,k) = 0.0
411ccccc pcltau(i,k) = 0.0
412ccccc pclemi(i,k) = 0.0
413ccccc ENDDO
414ccccc ENDDO
415C     
416C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
417C     
418c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
419c     initialisations
420      DO i=1,klon
421         zclear(i)=1.
422         zcloud(i)=0.
423         pch(i)=1.0
424         pcm(i) = 1.0
425         pcl(i) = 1.0
426         pctlwp(i) = 0.0
427      ENDDO
428C
429cIM cf CR DO k=1,klev
430      DO k = klev, 1, -1
431         DO i = 1, klon
432            pctlwp(i) = pctlwp(i)
433     &           + pqlwp(i,k)*diff_paprs(i,k)
434         ENDDO
435      ENDDO
436c     IM cf. CR
437      IF (NOVLP.EQ.1) THEN
438         DO k = klev, 1, -1
439            DO i = 1, klon
440               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
441     &              /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
442               pct(i)=1.-zclear(i)
443               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
444                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
445     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
446               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
447     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
448                  pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
449     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
450               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
451                  pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
452     &                 /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
453               endif
454               zcloud(i)=pclc(i,k)
455            ENDDO
456         ENDDO
457      ELSE IF (NOVLP.EQ.2) THEN
458         DO k = klev, 1, -1
459            DO i = 1, klon
460               zcloud(i)=MAX(pclc(i,k),zcloud(i))
461               pct(i)=zcloud(i)
462               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
463                  pch(i) = MIN(pclc(i,k),pch(i))
464               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
465     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
466                  pcm(i) = MIN(pclc(i,k),pcm(i))
467               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
468                  pcl(i) = MIN(pclc(i,k),pcl(i))
469               endif
470            ENDDO
471         ENDDO
472      ELSE IF (NOVLP.EQ.3) THEN
473         DO k = klev, 1, -1
474            DO i = 1, klon
475               zclear(i)=zclear(i)*(1.-pclc(i,k))
476               pct(i)=1-zclear(i)
477               IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
478                  pch(i) = pch(i)*(1.0-pclc(i,k))
479               ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
480     &                 pplay(i,k).LE.cetamb*paprs(i,1)) THEN
481                  pcm(i) = pcm(i)*(1.0-pclc(i,k))
482               ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
483                  pcl(i) = pcl(i)*(1.0-pclc(i,k))
484               endif
485            ENDDO
486         ENDDO
487      ENDIF
488     
489C     
490      DO i = 1, klon
491c     IM cf. CR          pct(i)=1.-pct(i)
492         pch(i)=1.-pch(i)
493         pcm(i)=1.-pcm(i)
494         pcl(i)=1.-pcl(i)
495      ENDDO
496
497c ========================================================
498! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
499c ========================================================
500!! change by Nicolas Yan (LSCE)
501!! Cloud Droplet Number Concentration (CDNC) : 3D variable
502!! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
503!! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
504!! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
505!! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
506      IF (ok_newmicro) THEN
507         IF (ok_aie) THEN
508            DO k = 1, klev
509               DO i = 1, klon
510                  phase3d(i,k)=1-zfice2(i,k)
511                  IF (pclc(i,k) .LE. seuil_neb) THEN
512                     lcc3d(i,k)=seuil_neb*phase3d(i,k)
513                  ELSE
514                     lcc3d(i,k)=pclc(i,k)*phase3d(i,k)
515                  ENDIF
516                  scdnc(i,k)=lcc3d(i,k)*cdnc(i,k) ! m-3
517               ENDDO
518            ENDDO
519
520            DO i=1,klon
521               lcc(i)=0.
522               reffclwtop(i)=0.
523               cldncl(i)=0.
524               IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i) = 1.
525               IF(MAXIMUM) tcc(i) = 0.
526            ENDDO
527     
528            FIRST=.TRUE.
529
530            DO i=1,klon
531               DO k=klev-1,1,-1 !From TOA down
532
533
534            ! Test, if the cloud optical depth exceeds the necessary
535            ! threshold:
536
537             IF (pcltau(i,k).GT.thres_tau .AND. pclc(i,k).GT.thres_neb)
538     .                                                             THEN
539               ! To calculate the right Temperature at cloud top,
540               ! interpolate it between layers:
541                  t_tmp = t(i,k) +
542     .              (paprs(i,k+1)-pplay(i,k))/(pplay(i,k+1)-pplay(i,k))
543     .              * ( t(i,k+1) - t(i,k) )
544
545                  IF(MAXIMUM) THEN
546                    IF(FIRST) THEN
547                       write(*,*)'Hypothese de recouvrement: MAXIMUM'
548                       FIRST=.FALSE.
549                    ENDIF
550                    flag_max= -1.
551                    ftmp(i) = MAX(tcc(i),pclc(i,k))
552                  ENDIF
553
554                  IF(RANDOM) THEN
555                    IF(FIRST) THEN
556                       write(*,*)'Hypothese de recouvrement: RANDOM'
557                       FIRST=.FALSE.
558                    ENDIF
559                    flag_max= 1.
560                    ftmp(i) = tcc(i) * (1-pclc(i,k))
561                  ENDIF
562
563                  IF(MAXIMUM_RANDOM) THEN
564                    IF(FIRST) THEN
565                       write(*,*)'Hypothese de recouvrement: MAXIMUM_
566     .                         RANDOM'
567                       FIRST=.FALSE.
568                    ENDIF
569                    flag_max= 1.
570                    ftmp(i) = tcc(i) *
571     .              (1. - MAX(pclc(i,k),pclc(i,k+1))) /
572     .              (1. - MIN(pclc(i,k+1),1.-thres_neb))
573                  ENDIF
574c Effective radius of cloud droplet at top of cloud (m)
575                  reffclwtop(i) = reffclwtop(i) + rad_chaud_tab(i,k) *
576     .           1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max
577c CDNC at top of cloud (m-3)
578                  cldncl(i) = cldncl(i) + cdnc(i,k) * phase3d(i,k) *
579     .                 (tcc(i) - ftmp(i))*flag_max
580c Liquid Cloud Content at top of cloud
581                  lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))*
582     .                    flag_max
583c Total Cloud Content at top of cloud
584                  tcc(i)=ftmp(i)
585             
586          ENDIF ! is there a visible, not-too-small cloud? 
587          ENDDO ! loop over k
588
589          IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
590         ENDDO ! loop over i
591
592!! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC  REFFCLWS)
593            DO i = 1, klon
594               DO k = 1, klev
595                  pqlwpcon(i,k)=rnebcon(i,k)*clwcon(i,k) ! fraction eau liquide convective
596                  pqlwpstra(i,k)=pclc(i,k)*phase3d(i,k)-pqlwpcon(i,k) ! fraction eau liquide stratiforme
597                  IF (pqlwpstra(i,k) .LE. 0.0) pqlwpstra(i,k)=0.0
598! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
599                  reffclwc(i,k)=1.1
600     &                 *((pqlwpcon(i,k)*pplay(i,k)/(RD * T(i,k)))
601     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
602                  reffclwc(i,k) = MAX(reffclwc(i,k) * 1e6, 5.)
603
604! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
605                  IF ((pclc(i,k)-rnebcon(i,k)) .LE. seuil_neb) THEN ! tout sous la forme convective
606                     reffclws(i,k)=0.0
607                     lcc3dstra(i,k)= 0.0
608                  ELSE
609                     reffclws(i,k) = (pclc(i,k)*phase3d(i,k)*
610     &                               rad_chaud_tab(i,k)-
611     &                            pqlwpcon(i,k)*reffclwc(i,k))
612                     IF(reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
613                     lcc3dstra(i,k)=pqlwpstra(i,k)
614                 ENDIF
615!Convertion from um to m
616                  IF(rnebcon(i,k). LE. seuil_neb) THEN
617                    reffclwc(i,k) = reffclwc(i,k)*seuil_neb*clwcon(i,k)
618     &                              *1.0E-06
619                    lcc3dcon(i,k)= seuil_neb*clwcon(i,k)
620                  ELSE
621                    reffclwc(i,k) = reffclwc(i,k)*pqlwpcon(i,k)
622     &                              *1.0E-06
623                    lcc3dcon(i,k) = pqlwpcon(i,k)
624                  ENDIF
625
626                  reffclws(i,k) = reffclws(i,k)*1.0E-06
627
628               ENDDO !klev
629            ENDDO !klon
630
631!! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
632            DO k = 1, klev
633               DO i = 1, klon
634                   pdel(i,k) = paprs(i,k)-paprs(i,k+1)
635                   zrho(i,k)=pplay(i,k)/t(i,k)/RD                  ! kg/m3
636                   dh(i,k)=pdel(i,k)/(gravit*zrho(i,k)) ! hauteur de chaque boite (m)
637               ENDDO
638            ENDDO
639c
640            DO i = 1, klon
641               cldnvi(i)=0.
642               lcc_integrat(i)=0.
643               height(i)=0.
644               DO k = 1, klev
645                  cldnvi(i)=cldnvi(i)+cdnc(i,k)*lcc3d(i,k)*dh(i,k)
646                  lcc_integrat(i)=lcc_integrat(i)+lcc3d(i,k)*dh(i,k)
647                  height(i)=height(i)+dh(i,k)
648               ENDDO ! klev
649               lcc_integrat(i)=lcc_integrat(i)/height(i)
650               IF (lcc_integrat(i) .LE. 1.0E-03) THEN
651                  cldnvi(i)=cldnvi(i)*lcc(i)/seuil_neb
652               ELSE
653                  cldnvi(i)=cldnvi(i)*lcc(i)/lcc_integrat(i)
654               ENDIF
655            ENDDO ! klon
656           
657            DO i = 1, klon
658               DO k = 1, klev
659                IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0
660                IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
661                IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0
662                IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0
663                IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0
664                IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
665               ENDDO
666               IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0
667               IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0
668               IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0
669               IF (lcc(i) .LE. 0.0) lcc(i)=0.0
670            ENDDO
671
672         ENDIF !ok_aie
673      ENDIF !ok newmicro
674c
675C
676      RETURN
677      END
Note: See TracBrowser for help on using the repository browser.