source: LMDZ5/trunk/libf/phylmd/newmicro.F @ 1651

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

Modifications concerning the cloud scheme:

  1. In newmicro, it now possible to read a min and max effective radius of ice particles from physiq.def : rei_min and rei_max which were initially set to 3.5 and 61.29 microns in the code.

concerns : conf_phys.F90, nuage.h, newmicro.F

  1. In physiq.F, in case of combination of iflag_cldcon>=5 (A. Jam cloud scheme) and iflag_coupl=5

concerns : physiq.F and thermcell_main.F90
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.5 KB
Line 
1! $Id: newmicro.F 1525 2011-05-25 10:55:27Z fairhead $
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      REAL d_rei_dt
158
159! Abderrahmane oct 2009
160      Real reliq(klon, klev), reice(klon, klev)
161
162!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
163! FH : 2011/05/24
164!
165! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max
166! to be used for a temperature in celcius T(°C) < 0
167! rei=rei_min for T(°C) < -81.4
168!
169! Calcul de la pente de la relation entre rayon effective des cristaux
170! et la température.
171! Pour retrouver les résultats numériques de la version d'origine,
172! on impose 0.71 quand on est proche de 0.71
173
174      d_rei_dt=(rei_max-rei_min)/81.4
175      if (abs(d_rei_dt-0.71)<1.e-4) d_rei_dt=0.71
176!      print*,'d_rei_dT ',d_rei_dt,rei_min,rei_max
177!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
178c
179c Calculer l'epaisseur optique et l'emmissivite des nuages
180c
181c     IM inversion des DO
182      xflwp = 0.d0
183      xfiwp = 0.d0
184      xflwc = 0.d0
185      xfiwc = 0.d0
186
187! Initialisation
188      reliq=0.
189      reice=0.
190
191      DO k = 1, klev
192         DO i = 1, klon
193            diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
194         ENDDO
195      ENDDO
196
197      IF (ok_newmicro) THEN
198
199
200         DO k = 1, klev
201            DO i = 1, klon
202c               zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
203               zfice2(i,k) = 1.0 - (t(i,k)-t_glace_min) /
204     &                             (t_glace_max-t_glace_min)
205               zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
206c     IM Total Liquid/Ice water content                                   
207               xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
208               xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
209c     IM In-Cloud Liquid/Ice water content
210c     xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
211c     xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
212            ENDDO
213         ENDDO
214
215         IF (ok_aie) THEN
216            DO k = 1, klev
217               DO i = 1, klon
218                                ! Formula "D" of Boucher and Lohmann, Tellus, 1995
219                                !             
220                  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
221     &               log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
222                                ! Cloud droplet number concentration (CDNC) is restricted
223                                ! to be within [20, 1000 cm^3]
224                                !
225                  cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
226                                !
227                                !
228                  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
229     &               log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))
230     &               *1.e6 !-m-3
231                  cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
232               ENDDO
233            ENDDO
234            DO k = 1, klev
235               DO i = 1, klon
236!                  rad_chaud_tab(i,k) =
237!     &                 MAX(1.1e6
238!     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
239!     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.)
240                  rad_chaud_tab(i,k) =
241     &                 1.1
242     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
243     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
244                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.)
245               ENDDO           
246            ENDDO
247         ELSE
248            DO k = 1, MIN(3,klev)
249               DO i = 1, klon
250                  rad_chaud_tab(i,k) = rad_chau2
251               ENDDO           
252            ENDDO
253            DO k = MIN(3,klev)+1, klev
254               DO i = 1, klon
255                  rad_chaud_tab(i,k) = rad_chau1
256               ENDDO           
257            ENDDO
258
259         ENDIF
260         
261         DO k = 1, klev
262!            IF(.not.ok_aie) THEN
263            rad_chaud = rad_chau1
264            IF (k.LE.3) rad_chaud = rad_chau2
265!            ENDIF
266            DO i = 1, klon
267               IF (pclc(i,k) .LE. seuil_neb) THEN
268               
269c     -- effective cloud droplet radius (microns):
270               
271c     for liquid water clouds:
272                                ! For output diagnostics
273                                !
274                                ! Cloud droplet effective radius [um]
275                                !
276                                ! we multiply here with f * xl (fraction of liquid water
277                                ! clouds in the grid cell) to avoid problems in the
278                                ! averaging of the output.
279                                ! In the output of IOIPSL, derive the real cloud droplet
280                                ! effective radius as re/fl
281                                !
282                                   
283                  fl(i,k) = seuil_neb*(1.-zfice2(i,k))           
284                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
285                 
286                  rel = 0.
287                  rei = 0.
288                  pclc(i,k) = 0.0
289                  pcltau(i,k) = 0.0
290                  pclemi(i,k) = 0.0
291                  cldtaupi(i,k) = 0.0                 
292               ELSE
293
294c     -- liquid/ice cloud water paths:
295                 
296                  zflwp_var= 1000.*(1.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k)
297     &                 *diff_paprs(i,k)
298                  zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k)
299     &                 *diff_paprs(i,k)
300                 
301c     -- effective cloud droplet radius (microns):
302               
303c     for liquid water clouds:
304                                   
305                  IF (ok_aie) THEN
306                     radius =
307     &                    1.1
308     &                    *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 
309     &                    /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.)
310                     radius = MAX(radius*1e6, 5.)
311                 
312                     tc = t(i,k)-273.15
313                     rei = d_rei_dt*tc + rei_max
314                     if (tc.le.-81.4) rei = rei_min
315                     if (zflwp_var.eq.0.) radius = 1.
316                     if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
317                     cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
318     &                    + zfiwp_var * (3.448e-03  + 2.431/rei)
319
320                  ENDIF         ! ok_aie
321                                ! For output diagnostics
322                                !
323                                ! Cloud droplet effective radius [um]
324                                !
325                                ! we multiply here with f * xl (fraction of liquid water
326                                ! clouds in the grid cell) to avoid problems in the
327                                ! averaging of the output.
328                                ! In the output of IOIPSL, derive the real cloud droplet
329                                ! effective radius as re/fl
330                                !
331 
332                  fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))           
333                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
334                 
335                  rel = rad_chaud_tab(i,k)
336c     for ice clouds: as a function of the ambiant temperature
337c     [formula used by Iacobellis and Somerville (2000), with an
338c     asymptotical value of 3.5 microns at T<-81.4 C added to be
339c     consistent with observations of Heymsfield et al. 1986]:
340c  2011/05/24 : rei_min = 3.5 becomes a free parameter as well as rei_max=61.29
341                  tc = t(i,k)-273.15
342                  rei = d_rei_dt*tc + rei_max
343                  if (tc.le.-81.4) rei = rei_min
344c     -- cloud optical thickness :
345               
346c     [for liquid clouds, traditional formula,
347c     for ice clouds, Ebert & Curry (1992)]
348                 
349                 if (zflwp_var.eq.0.) rel = 1.
350                 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
351                 pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
352     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
353c     -- cloud infrared emissivity:
354               
355c     [the broadband infrared absorption coefficient is parameterized
356c     as a function of the effective cld droplet radius]
357               
358c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
359                  k_ice = k_ice0 + 1.0/rei
360                 
361                  pclemi(i,k) = 1.0
362     &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
363
364               ENDIF
365              reliq(i,k)=rel
366              reice(i,k)=rei
367!              if (i.eq.1) then
368!              print*,'Dans newmicro rel, rei :',rel, rei
369!              print*,'Dans newmicro reliq, reice :',
370!     $             reliq(i,k),reice(i,k)
371!              endif
372
373            ENDDO
374         ENDDO
375
376         DO k = 1, klev
377            DO i = 1, klon
378               xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
379               xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
380            ENDDO
381         ENDDO
382
383      ELSE
384         DO k = 1, klev
385            rad_chaud = rad_chau1
386            IF (k.LE.3) rad_chaud = rad_chau2
387            DO i = 1, klon
388                             
389               IF (pclc(i,k) .LE. seuil_neb) THEN
390
391                  pclc(i,k) = 0.0
392                  pcltau(i,k) = 0.0
393                  pclemi(i,k) = 0.0
394                  cldtaupi(i,k) = 0.0
395
396               ELSE
397
398                  zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k)
399     &                 /pclc(i,k)
400                 
401                  zfice1 = MIN(
402     &                 MAX( 1.0 - (t(i,k)-t_glace_min) /
403     &                    (t_glace_max-t_glace_min),0.0),1.0)**nexpo
404                 
405                  radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
406                  coef   = coef_chau * (1.-zfice1) + coef_froi * zfice1
407
408                  pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
409                  pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
410
411               ENDIF
412                             
413            ENDDO
414         ENDDO
415      ENDIF
416     
417      IF (.NOT.ok_aie) THEN
418         DO k = 1, klev
419            DO i = 1, klon
420               cldtaupi(i,k)=pcltau(i,k)
421            ENDDO
422         ENDDO               
423      ENDIF
424
425ccc   DO k = 1, klev
426ccc   DO i = 1, klon
427ccc   t(i,k) = t(i,k)
428ccc   pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
429ccc   lo = pclc(i,k) .GT. (2.*1.e-5)
430ccc   zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
431ccc   .          /(rg*pclc(i,k))
432ccc   zradef = 10.0 + (1.-sigs(k))*45.0
433ccc   pcltau(i,k) = 1.5 * zflwp / zradef
434ccc   zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
435ccc   zmsac = 0.13*(1.0-zfice) + 0.08*zfice
436ccc   pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
437ccc   if (.NOT.lo) pclc(i,k) = 0.0
438ccc   if (.NOT.lo) pcltau(i,k) = 0.0
439ccc   if (.NOT.lo) pclemi(i,k) = 0.0
440ccc   ENDDO
441ccc   ENDDO
442ccccc print*, 'pas de nuage dans le rayonnement'
443ccccc DO k = 1, klev
444ccccc DO i = 1, klon
445ccccc pclc(i,k) = 0.0
446ccccc pcltau(i,k) = 0.0
447ccccc pclemi(i,k) = 0.0
448ccccc ENDDO
449ccccc ENDDO
450C     
451C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
452C     
453c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
454c     initialisations
455      DO i=1,klon
456         zclear(i)=1.
457         zcloud(i)=0.
458         zcloudh(i)=0.
459         zcloudm(i)=0.
460         zcloudl(i)=0.
461         pch(i)=1.0
462         pcm(i) = 1.0
463         pcl(i) = 1.0
464         pctlwp(i) = 0.0
465      ENDDO
466C
467cIM cf CR DO k=1,klev
468      DO k = klev, 1, -1
469         DO i = 1, klon
470            pctlwp(i) = pctlwp(i)
471     &           + pqlwp(i,k)*diff_paprs(i,k)
472         ENDDO
473      ENDDO
474c     IM cf. CR
475      IF (NOVLP.EQ.1) THEN
476         DO k = klev, 1, -1
477            DO i = 1, klon
478               zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
479     &              /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
480               pct(i)=1.-zclear(i)
481               IF (paprs(i,k).LT.prmhc) THEN
482                  pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloudh(i)))
483     &                 /(1.-MIN(real(zcloudh(i), kind=8),1.-ZEPSEC))
484                  zcloudh(i)=pclc(i,k)
485               ELSE IF (paprs(i,k).GE.prmhc .AND.
486     &                 paprs(i,k).LT.prlmc) THEN
487                  pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloudm(i)))
488     &                 /(1.-MIN(real(zcloudm(i), kind=8),1.-ZEPSEC))
489                  zcloudm(i)=pclc(i,k)
490               ELSE IF (paprs(i,k).GE.prlmc) THEN
491                  pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloudl(i)))
492     &                 /(1.-MIN(real(zcloudl(i), kind=8),1.-ZEPSEC))
493                  zcloudl(i)=pclc(i,k)
494               endif
495               zcloud(i)=pclc(i,k)
496            ENDDO
497         ENDDO
498      ELSE IF (NOVLP.EQ.2) THEN
499         DO k = klev, 1, -1
500            DO i = 1, klon
501               zcloud(i)=MAX(pclc(i,k),zcloud(i))
502               pct(i)=zcloud(i)
503               IF (paprs(i,k).LT.prmhc) THEN
504                  pch(i) = MIN(pclc(i,k),pch(i))
505               ELSE IF (paprs(i,k).GE.prmhc .AND.
506     &                 paprs(i,k).LT.prlmc) THEN
507                  pcm(i) = MIN(pclc(i,k),pcm(i))
508               ELSE IF (paprs(i,k).GE.prlmc) THEN
509                  pcl(i) = MIN(pclc(i,k),pcl(i))
510               endif
511            ENDDO
512         ENDDO
513      ELSE IF (NOVLP.EQ.3) THEN
514         DO k = klev, 1, -1
515            DO i = 1, klon
516               zclear(i)=zclear(i)*(1.-pclc(i,k))
517               pct(i)=1-zclear(i)
518               IF (paprs(i,k).LT.prmhc) THEN
519                  pch(i) = pch(i)*(1.0-pclc(i,k))
520               ELSE IF (paprs(i,k).GE.prmhc .AND.
521     &                 paprs(i,k).LT.prlmc) THEN
522                  pcm(i) = pcm(i)*(1.0-pclc(i,k))
523               ELSE IF (paprs(i,k).GE.prlmc) THEN
524                  pcl(i) = pcl(i)*(1.0-pclc(i,k))
525               endif
526            ENDDO
527         ENDDO
528      ENDIF
529     
530C     
531      DO i = 1, klon
532c     IM cf. CR          pct(i)=1.-pct(i)
533         pch(i)=1.-pch(i)
534         pcm(i)=1.-pcm(i)
535         pcl(i)=1.-pcl(i)
536      ENDDO
537
538c ========================================================
539! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
540c ========================================================
541!! change by Nicolas Yan (LSCE)
542!! Cloud Droplet Number Concentration (CDNC) : 3D variable
543!! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
544!! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
545!! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
546!! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
547      IF (ok_newmicro) THEN
548         IF (ok_aie) THEN
549            DO k = 1, klev
550               DO i = 1, klon
551                  phase3d(i,k)=1-zfice2(i,k)
552                  IF (pclc(i,k) .LE. seuil_neb) THEN
553                     lcc3d(i,k)=seuil_neb*phase3d(i,k)
554                  ELSE
555                     lcc3d(i,k)=pclc(i,k)*phase3d(i,k)
556                  ENDIF
557                  scdnc(i,k)=lcc3d(i,k)*cdnc(i,k) ! m-3
558               ENDDO
559            ENDDO
560
561            DO i=1,klon
562               lcc(i)=0.
563               reffclwtop(i)=0.
564               cldncl(i)=0.
565               IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i) = 1.
566               IF(MAXIMUM) tcc(i) = 0.
567            ENDDO
568     
569
570            DO i=1,klon
571               DO k=klev-1,1,-1 !From TOA down
572
573
574            ! Test, if the cloud optical depth exceeds the necessary
575            ! threshold:
576
577             IF (pcltau(i,k).GT.thres_tau .AND. pclc(i,k).GT.thres_neb)
578     .                                                             THEN
579               ! To calculate the right Temperature at cloud top,
580               ! interpolate it between layers:
581                  t_tmp = t(i,k) +
582     .              (paprs(i,k+1)-pplay(i,k))/(pplay(i,k+1)-pplay(i,k))
583     .              * ( t(i,k+1) - t(i,k) )
584
585                  IF(MAXIMUM) THEN
586                    IF(FIRST) THEN
587                       write(*,*)'Hypothese de recouvrement: MAXIMUM'
588                       FIRST=.FALSE.
589                    ENDIF
590                    flag_max= -1.
591                    ftmp(i) = MAX(tcc(i),pclc(i,k))
592                  ENDIF
593
594                  IF(RANDOM) THEN
595                    IF(FIRST) THEN
596                       write(*,*)'Hypothese de recouvrement: RANDOM'
597                       FIRST=.FALSE.
598                    ENDIF
599                    flag_max= 1.
600                    ftmp(i) = tcc(i) * (1-pclc(i,k))
601                  ENDIF
602
603                  IF(MAXIMUM_RANDOM) THEN
604                    IF(FIRST) THEN
605                       write(*,*)'Hypothese de recouvrement: MAXIMUM_
606     .                         RANDOM'
607                       FIRST=.FALSE.
608                    ENDIF
609                    flag_max= 1.
610                    ftmp(i) = tcc(i) *
611     .              (1. - MAX(pclc(i,k),pclc(i,k+1))) /
612     .              (1. - MIN(pclc(i,k+1),1.-thres_neb))
613                  ENDIF
614c Effective radius of cloud droplet at top of cloud (m)
615                  reffclwtop(i) = reffclwtop(i) + rad_chaud_tab(i,k) *
616     .           1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max
617c CDNC at top of cloud (m-3)
618                  cldncl(i) = cldncl(i) + cdnc(i,k) * phase3d(i,k) *
619     .                 (tcc(i) - ftmp(i))*flag_max
620c Liquid Cloud Content at top of cloud
621                  lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))*
622     .                    flag_max
623c Total Cloud Content at top of cloud
624                  tcc(i)=ftmp(i)
625             
626          ENDIF ! is there a visible, not-too-small cloud? 
627          ENDDO ! loop over k
628
629          IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
630         ENDDO ! loop over i
631
632!! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC  REFFCLWS)
633            DO i = 1, klon
634               DO k = 1, klev
635                  pqlwpcon(i,k)=rnebcon(i,k)*clwcon(i,k) ! fraction eau liquide convective
636                  pqlwpstra(i,k)=pclc(i,k)*phase3d(i,k)-pqlwpcon(i,k) ! fraction eau liquide stratiforme
637                  IF (pqlwpstra(i,k) .LE. 0.0) pqlwpstra(i,k)=0.0
638! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
639                  reffclwc(i,k)=1.1
640     &                 *((pqlwpcon(i,k)*pplay(i,k)/(RD * T(i,k)))
641     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
642                  reffclwc(i,k) = MAX(reffclwc(i,k) * 1e6, 5.)
643
644! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
645                  IF ((pclc(i,k)-rnebcon(i,k)) .LE. seuil_neb) THEN ! tout sous la forme convective
646                     reffclws(i,k)=0.0
647                     lcc3dstra(i,k)= 0.0
648                  ELSE
649                     reffclws(i,k) = (pclc(i,k)*phase3d(i,k)*
650     &                               rad_chaud_tab(i,k)-
651     &                            pqlwpcon(i,k)*reffclwc(i,k))
652                     IF(reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
653                     lcc3dstra(i,k)=pqlwpstra(i,k)
654                 ENDIF
655!Convertion from um to m
656                  IF(rnebcon(i,k). LE. seuil_neb) THEN
657                    reffclwc(i,k) = reffclwc(i,k)*seuil_neb*clwcon(i,k)
658     &                              *1.0E-06
659                    lcc3dcon(i,k)= seuil_neb*clwcon(i,k)
660                  ELSE
661                    reffclwc(i,k) = reffclwc(i,k)*pqlwpcon(i,k)
662     &                              *1.0E-06
663                    lcc3dcon(i,k) = pqlwpcon(i,k)
664                  ENDIF
665
666                  reffclws(i,k) = reffclws(i,k)*1.0E-06
667
668               ENDDO !klev
669            ENDDO !klon
670
671!! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
672            DO k = 1, klev
673               DO i = 1, klon
674                   pdel(i,k) = paprs(i,k)-paprs(i,k+1)
675                   zrho(i,k)=pplay(i,k)/t(i,k)/RD                  ! kg/m3
676                   dh(i,k)=pdel(i,k)/(gravit*zrho(i,k)) ! hauteur de chaque boite (m)
677               ENDDO
678            ENDDO
679c
680            DO i = 1, klon
681               cldnvi(i)=0.
682               lcc_integrat(i)=0.
683               height(i)=0.
684               DO k = 1, klev
685                  cldnvi(i)=cldnvi(i)+cdnc(i,k)*lcc3d(i,k)*dh(i,k)
686                  lcc_integrat(i)=lcc_integrat(i)+lcc3d(i,k)*dh(i,k)
687                  height(i)=height(i)+dh(i,k)
688               ENDDO ! klev
689               lcc_integrat(i)=lcc_integrat(i)/height(i)
690               IF (lcc_integrat(i) .LE. 1.0E-03) THEN
691                  cldnvi(i)=cldnvi(i)*lcc(i)/seuil_neb
692               ELSE
693                  cldnvi(i)=cldnvi(i)*lcc(i)/lcc_integrat(i)
694               ENDIF
695            ENDDO ! klon
696           
697            DO i = 1, klon
698               DO k = 1, klev
699                IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0
700                IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
701                IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0
702                IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0
703                IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0
704                IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
705               ENDDO
706               IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0
707               IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0
708               IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0
709               IF (lcc(i) .LE. 0.0) lcc(i)=0.0
710            ENDDO
711
712         ENDIF !ok_aie
713      ENDIF !ok newmicro
714c
715C
716      RETURN
717      END
Note: See TracBrowser for help on using the repository browser.