source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/newmicro.F @ 1522

Last change on this file since 1522 was 1522, checked in by idelkadi, 14 years ago

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