source: LMDZ5/branches/LMDZ5-dev032011/libf/phylmd/newmicro.F @ 5872

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

Un peu de ménage sur les prints

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