source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/nuage.F90 @ 5233

Last change on this file since 5233 was 4727, checked in by idelkadi, 13 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.7 KB
Line 
1! $Id: nuage.F90 4727 2023-10-19 14:02:57Z abarral $
2
3SUBROUTINE nuage(paprs, pplay, t, pqlwp,picefra, pclc, pcltau, pclemi, pch, pcl, pcm, &
4    pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, distcltop, &
5    temp_cltop, cldtaupi, re, fl)
6  USE dimphy
7  USE lmdz_lscp_tools, only: icefrac_lscp
8  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
9  USE lmdz_lscp_ini, only : iflag_t_glace
10  USE phys_local_var_mod, ONLY: ptconv
11  IMPLICIT NONE
12  ! ======================================================================
13  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
14  ! Objet: Calculer epaisseur optique et emmissivite des nuages
15  ! ======================================================================
16  ! Arguments:
17  ! t-------input-R-temperature
18  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
19  ! picefra--inout-R-fraction de glace dans les nuages (-)
20  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
21  ! ok_aie--input-L-apply aerosol indirect effect or not
22  ! mass_solu_aero-----input-R-total mass concentration for all soluble
23  ! aerosols[ug/m^3]
24  ! mass_solu_aero_pi--input-R-dito, pre-industrial value
25  ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
26  ! bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
27
28  ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
29  ! needed for the diagnostics of the aerosol indirect
30  ! radiative forcing (see radlwsw)
31  ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
32  ! fl------output-R-Denominator to re, introduced to avoid problems in
33  ! the averaging of the output. fl is the fraction of liquid
34  ! water clouds within a grid cell
35
36  ! pcltau--output-R-epaisseur optique des nuages
37  ! pclemi--output-R-emissivite des nuages (0 a 1)
38  ! ======================================================================
39
40  include "YOMCST.h"
41  include "nuage.h" ! JBM 3/14
42  include "clesphys.h"
43
44  REAL paprs(klon, klev+1), pplay(klon, klev)
45  REAL t(klon, klev)
46
47  REAL pclc(klon, klev)
48  REAL pqlwp(klon, klev), picefra(klon,klev)
49  REAL pcltau(klon, klev), pclemi(klon, klev)
50
51  REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
52  REAL distcltop(klon,klev)
53  REAL temp_cltop(klon,klev)
54  LOGICAL lo
55
56  REAL cetahb, cetamb
57  PARAMETER (cetahb=0.45, cetamb=0.80)
58
59  INTEGER i, k
60  REAL zflwp, zradef, zfice(klon), zmsac
61
62  REAL radius, rad_chaud
63! JBM (3/14) parameters already defined in nuage.h:
64! REAL rad_froid, rad_chau1, rad_chau2
65! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
66  ! cc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
67  ! sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
68  REAL coef, coef_froi, coef_chau
69  PARAMETER (coef_chau=0.13, coef_froi=0.09)
70  REAL seuil_neb
71  PARAMETER (seuil_neb=0.001)
72! JBM (3/14) nexpo is replaced by exposant_glace
73! REAL nexpo ! exponentiel pour glace/eau
74! PARAMETER (nexpo=6.)
75  REAL, PARAMETER :: t_glace_min_old = 258.
76  INTEGER, PARAMETER :: exposant_glace_old = 6
77
78
79  ! jq for the aerosol indirect effect
80  ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
81  ! jq
82  LOGICAL ok_aie ! Apply AIE or not?
83
84  REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols[ug m-3]
85  REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
86  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
87  REAL re(klon, klev) ! cloud droplet effective radius [um]
88  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
89  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
90
91  REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds
92  ! within the grid cell)
93
94  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
95
96  REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
97  REAl dzfice(klon)
98  ! jq-end
99
100  ! cc      PARAMETER (nexpo=1)
101
102  ! Calculer l'epaisseur optique et l'emmissivite des nuages
103
104  DO k = 1, klev
105     IF (iflag_t_glace.EQ.0) THEN
106       DO i = 1, klon
107        zfice(i) = 1.0 - (t(i,k)-t_glace_min_old)/(273.13-t_glace_min_old)
108        zfice(i) = min(max(zfice(i),0.0), 1.0)
109        zfice(i) = zfice(i)**exposant_glace_old
110       ENDDO
111     ELSE ! of IF (iflag_t_glace.EQ.0)
112! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
113!       zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
114!                           t_glace_max, exposant_glace)
115        IF (ok_new_lscp) THEN
116            CALL icefrac_lscp(klon,t(:,k),iflag_ice_thermo,distcltop(:,k),temp_cltop(:,k),zfice(:),dzfice(:))
117        ELSE
118            CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
119
120        ENDIF
121
122    IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
123    ! EV: take the ice fraction directly from the lscp code
124    ! consistent only for non convective grid points
125    ! critical for mixed phase clouds
126        DO i=1,klon
127        IF (.NOT. ptconv(i,k)) THEN
128           zfice(i)=picefra(i,k)
129        ENDIF
130        ENDDO
131    ENDIF
132
133
134     ENDIF
135
136    DO i = 1, klon
137      rad_chaud = rad_chau1
138      IF (k<=3) rad_chaud = rad_chau2
139
140      pclc(i, k) = max(pclc(i,k), seuil_neb)
141      zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))
142
143      IF (ok_aie) THEN
144          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
145          !
146        cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
147          1.E-4))/log(10.))*1.E6 !-m-3
148          ! Cloud droplet number concentration (CDNC) is restricted
149          ! to be within [20, 1000 cm^3]
150          !
151        cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
152        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
153          1.E-4))/log(10.))*1.E6 !-m-3
154        cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
155          !
156          !
157          ! air density: pplay(i,k) / (RD * zT(i,k))
158          ! factor 1.1: derive effective radius from volume-mean radius
159          ! factor 1000 is the water density
160          ! _chaud means that this is the CDR for liquid water clouds
161          !
162        rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &
163          *cdnc(i,k)))**(1./3.)
164          !
165          ! Convert to um. CDR shall be at least 3 um.
166          !
167        rad_chaud = max(rad_chaud*1.E6, 3.)
168
169          ! For output diagnostics
170          !
171          ! Cloud droplet effective radius [um]
172          !
173          ! we multiply here with f * xl (fraction of liquid water
174          ! clouds in the grid cell) to avoid problems in the
175          ! averaging of the output.
176          ! In the output of IOIPSL, derive the real cloud droplet
177          ! effective radius as re/fl
178          !
179        fl(i, k) = pclc(i, k)*(1.-zfice(i))
180        re(i, k) = rad_chaud*fl(i, k)
181
182          ! Pre-industrial cloud opt thickness
183          !
184          ! "radius" is calculated as rad_chaud above (plus the
185          ! ice cloud contribution) but using cdnc_pi instead of
186          ! cdnc.
187        radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &
188          1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i)
189        cldtaupi(i, k) = 3.0/2.0*zflwp/radius
190      END IF ! ok_aie
191
192      radius = rad_chaud*(1.-zfice(i)) + rad_froid*zfice(i)
193      coef = coef_chau*(1.-zfice(i)) + coef_froi*zfice(i)
194      pcltau(i, k) = 3.0/2.0*zflwp/radius
195      pclemi(i, k) = 1.0 - exp(-coef*zflwp)
196      lo = (pclc(i,k)<=seuil_neb)
197      IF (lo) pclc(i, k) = 0.0
198      IF (lo) pcltau(i, k) = 0.0
199      IF (lo) pclemi(i, k) = 0.0
200
201      IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k)
202    END DO
203  END DO
204  ! cc      DO k = 1, klev
205  ! cc      DO i = 1, klon
206  ! cc         t(i,k) = t(i,k)
207  ! cc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
208  ! cc         lo = pclc(i,k) .GT. (2.*1.e-5)
209  ! cc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
210  ! cc     .          /(rg*pclc(i,k))
211  ! cc         zradef = 10.0 + (1.-sigs(k))*45.0
212  ! cc         pcltau(i,k) = 1.5 * zflwp / zradef
213  ! cc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
214  ! cc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
215  ! cc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
216  ! cc         if (.NOT.lo) pclc(i,k) = 0.0
217  ! cc         if (.NOT.lo) pcltau(i,k) = 0.0
218  ! cc         if (.NOT.lo) pclemi(i,k) = 0.0
219  ! cc      ENDDO
220  ! cc      ENDDO
221  ! ccccc      print*, 'pas de nuage dans le rayonnement'
222  ! ccccc      DO k = 1, klev
223  ! ccccc      DO i = 1, klon
224  ! ccccc         pclc(i,k) = 0.0
225  ! ccccc         pcltau(i,k) = 0.0
226  ! ccccc         pclemi(i,k) = 0.0
227  ! ccccc      ENDDO
228  ! ccccc      ENDDO
229
230  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
231
232  DO i = 1, klon
233    pct(i) = 1.0
234    pch(i) = 1.0
235    pcm(i) = 1.0
236    pcl(i) = 1.0
237    pctlwp(i) = 0.0
238  END DO
239
240  DO k = klev, 1, -1
241    DO i = 1, klon
242      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
243      pct(i) = pct(i)*(1.0-pclc(i,k))
244      IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))
245      IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &
246        pcm(i) = pcm(i)*(1.0-pclc(i,k))
247      IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))
248    END DO
249  END DO
250
251  DO i = 1, klon
252    pct(i) = 1. - pct(i)
253    pch(i) = 1. - pch(i)
254    pcm(i) = 1. - pcm(i)
255    pcl(i) = 1. - pcl(i)
256  END DO
257
258  RETURN
259END SUBROUTINE nuage
260SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
261  USE dimphy
262  IMPLICIT NONE
263
264  ! Laurent Li (LMD/CNRS), le 12 octobre 1998
265  ! (adaptation du code ECMWF)
266
267  ! Dans certains cas, le schema pronostique des nuages n'est
268  ! pas suffisament performant. On a donc besoin de diagnostiquer
269  ! ces nuages. Je dois avouer que c'est une frustration.
270
271  include "YOMCST.h"
272
273  ! Arguments d'entree:
274  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
275  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
276  REAL t(klon, klev) ! temperature (K)
277  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
278  REAL rain(klon) ! pluie convective (kg/m2/s)
279  REAL snow(klon) ! neige convective (kg/m2/s)
280  INTEGER ktop(klon) ! sommet de la convection
281  INTEGER kbot(klon) ! bas de la convection
282
283  ! Arguments de sortie:
284  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
285  REAL dialiq(klon, klev) ! eau liquide nuageuse
286
287  ! Constantes ajustables:
288  REAL canva, canvb, canvh
289  PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
290  REAL cca, ccb, ccc
291  PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
292  REAL ccfct, ccscal
293  PARAMETER (ccfct=0.400)
294  PARAMETER (ccscal=1.0E+11)
295  REAL cetahb, cetamb
296  PARAMETER (cetahb=0.45, cetamb=0.80)
297  REAL cclwmr
298  PARAMETER (cclwmr=1.E-04)
299  REAL zepscr
300  PARAMETER (zepscr=1.0E-10)
301
302  ! Variables locales:
303  INTEGER i, k
304  REAL zcc(klon)
305
306  ! Initialisation:
307
308  DO k = 1, klev
309    DO i = 1, klon
310      diafra(i, k) = 0.0
311      dialiq(i, k) = 0.0
312    END DO
313  END DO
314
315  DO i = 1, klon ! Calculer la fraction nuageuse
316    zcc(i) = 0.0
317    IF ((rain(i)+snow(i))>0.) THEN
318      zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
319      zcc(i) = min(ccc, max(0.0,zcc(i)))
320    END IF
321  END DO
322
323  DO i = 1, klon ! pour traiter les enclumes
324    diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
325    IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
326      1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
327      i)*ccfct,canva*(zcc(i)-canvb)))
328    dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
329  END DO
330
331  DO k = 1, klev ! nuages convectifs (sauf enclumes)
332    DO i = 1, klon
333      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
334        diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
335        dialiq(i, k) = cclwmr*diafra(i, k)
336      END IF
337    END DO
338  END DO
339
340  RETURN
341END SUBROUTINE diagcld1
342SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
343  USE dimphy
344  IMPLICIT NONE
345
346  include "YOMCST.h"
347
348  ! Arguments d'entree:
349  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
350  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
351  REAL t(klon, klev) ! temperature (K)
352  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
353
354  ! Arguments de sortie:
355  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
356  REAL dialiq(klon, klev) ! eau liquide nuageuse
357
358  REAL cetamb
359  PARAMETER (cetamb=0.80)
360  REAL cloia, cloib, cloic, cloid
361  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
362  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
363  REAL rgammas
364  PARAMETER (rgammas=0.05)
365  REAL crhl
366  PARAMETER (crhl=0.15)
367  ! cc      PARAMETER (CRHL=0.70)
368  REAL t_coup
369  PARAMETER (t_coup=234.0)
370
371  ! Variables locales:
372  INTEGER i, k, kb, invb(klon)
373  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
374  REAL zdelta, zcor
375
376  ! Fonctions thermodynamiques:
377  include "YOETHF.h"
378  include "FCTTRE.h"
379
380  ! Initialisation:
381
382  DO k = 1, klev
383    DO i = 1, klon
384      diafra(i, k) = 0.0
385      dialiq(i, k) = 0.0
386    END DO
387  END DO
388
389  DO i = 1, klon
390    invb(i) = klev
391    zdthmin(i) = 0.0
392  END DO
393
394  DO k = 2, klev/2 - 1
395    DO i = 1, klon
396      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
397        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
398      zdthdp = zdthdp*cloia
399      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
400        zdthmin(i) = zdthdp
401        invb(i) = k
402      END IF
403    END DO
404  END DO
405
406  DO i = 1, klon
407    kb = invb(i)
408    IF (thermcep) THEN
409      zdelta = max(0., sign(1.,rtt-t(i,kb)))
410      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
411      zqs = min(0.5, zqs)
412      zcor = 1./(1.-retv*zqs)
413      zqs = zqs*zcor
414    ELSE
415      IF (t(i,kb)<t_coup) THEN
416        zqs = qsats(t(i,kb))/pplay(i, kb)
417      ELSE
418        zqs = qsatl(t(i,kb))/pplay(i, kb)
419      END IF
420    END IF
421    zcll = cloib*zdthmin(i) + cloic
422    zcll = min(1.0, max(0.0,zcll))
423    zrhb = q(i, kb)/zqs
424    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
425    zcll = min(1.0, max(0.0,zcll))
426    diafra(i, kb) = max(diafra(i,kb), zcll)
427    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
428  END DO
429
430  RETURN
431END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.