source: LMDZ6/trunk/libf/phylmd/nuage.f90 @ 5300

Last change on this file since 5300 was 5285, checked in by abarral, 5 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules
  • 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 5285 2024-10-28 13:33:29Z 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 clesphys_mod_h
7  USE yomcst_mod_h
8  USE dimphy
9  USE lmdz_lscp_tools, only: icefrac_lscp
10  USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
11  USE lmdz_lscp_ini, only : iflag_t_glace
12  USE phys_local_var_mod, ONLY: ptconv
13  IMPLICIT NONE
14  ! ======================================================================
15  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
16  ! Objet: Calculer epaisseur optique et emmissivite des nuages
17  ! ======================================================================
18  ! Arguments:
19  ! t-------input-R-temperature
20  ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
21  ! picefra--inout-R-fraction de glace dans les nuages (-)
22  ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
23  ! ok_aie--input-L-apply aerosol indirect effect or not
24  ! mass_solu_aero-----input-R-total mass concentration for all soluble
25  ! aerosols[ug/m^3]
26  ! mass_solu_aero_pi--input-R-dito, pre-industrial value
27  ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
28  ! bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
29
30  ! cldtaupi-output-R-pre-industrial value of cloud optical thickness,
31  ! needed for the diagnostics of the aerosol indirect
32  ! radiative forcing (see radlwsw)
33  ! re------output-R-Cloud droplet effective radius multiplied by fl [um]
34  ! fl------output-R-Denominator to re, introduced to avoid problems in
35  ! the averaging of the output. fl is the fraction of liquid
36  ! water clouds within a grid cell
37
38  ! pcltau--output-R-epaisseur optique des nuages
39  ! pclemi--output-R-emissivite des nuages (0 a 1)
40  ! ======================================================================
41
42  include "nuage.h" ! JBM 3/14
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  USE yomcst_mod_h
263IMPLICIT NONE
264
265  ! Laurent Li (LMD/CNRS), le 12 octobre 1998
266  ! (adaptation du code ECMWF)
267
268  ! Dans certains cas, le schema pronostique des nuages n'est
269  ! pas suffisament performant. On a donc besoin de diagnostiquer
270  ! ces nuages. Je dois avouer que c'est une frustration.
271
272
273
274  ! Arguments d'entree:
275  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
276  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
277  REAL t(klon, klev) ! temperature (K)
278  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
279  REAL rain(klon) ! pluie convective (kg/m2/s)
280  REAL snow(klon) ! neige convective (kg/m2/s)
281  INTEGER ktop(klon) ! sommet de la convection
282  INTEGER kbot(klon) ! bas de la convection
283
284  ! Arguments de sortie:
285  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
286  REAL dialiq(klon, klev) ! eau liquide nuageuse
287
288  ! Constantes ajustables:
289  REAL canva, canvb, canvh
290  PARAMETER (canva=2.0, canvb=0.3, canvh=0.4)
291  REAL cca, ccb, ccc
292  PARAMETER (cca=0.125, ccb=1.5, ccc=0.8)
293  REAL ccfct, ccscal
294  PARAMETER (ccfct=0.400)
295  PARAMETER (ccscal=1.0E+11)
296  REAL cetahb, cetamb
297  PARAMETER (cetahb=0.45, cetamb=0.80)
298  REAL cclwmr
299  PARAMETER (cclwmr=1.E-04)
300  REAL zepscr
301  PARAMETER (zepscr=1.0E-10)
302
303  ! Variables locales:
304  INTEGER i, k
305  REAL zcc(klon)
306
307  ! Initialisation:
308
309  DO k = 1, klev
310    DO i = 1, klon
311      diafra(i, k) = 0.0
312      dialiq(i, k) = 0.0
313    END DO
314  END DO
315
316  DO i = 1, klon ! Calculer la fraction nuageuse
317    zcc(i) = 0.0
318    IF ((rain(i)+snow(i))>0.) THEN
319      zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb
320      zcc(i) = min(ccc, max(0.0,zcc(i)))
321    END IF
322  END DO
323
324  DO i = 1, klon ! pour traiter les enclumes
325    diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct)
326    IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, &
327      1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( &
328      i)*ccfct,canva*(zcc(i)-canvb)))
329    dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i))
330  END DO
331
332  DO k = 1, klev ! nuages convectifs (sauf enclumes)
333    DO i = 1, klon
334      IF (k<ktop(i) .AND. k>=kbot(i)) THEN
335        diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct)
336        dialiq(i, k) = cclwmr*diafra(i, k)
337      END IF
338    END DO
339  END DO
340
341  RETURN
342END SUBROUTINE diagcld1
343SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
344  USE dimphy
345  USE yomcst_mod_h
346  USE yoethf_mod_h
347IMPLICIT NONE
348
349
350
351  ! Arguments d'entree:
352  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
353  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
354  REAL t(klon, klev) ! temperature (K)
355  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
356
357  ! Arguments de sortie:
358  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
359  REAL dialiq(klon, klev) ! eau liquide nuageuse
360
361  REAL cetamb
362  PARAMETER (cetamb=0.80)
363  REAL cloia, cloib, cloic, cloid
364  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
365  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
366  REAL rgammas
367  PARAMETER (rgammas=0.05)
368  REAL crhl
369  PARAMETER (crhl=0.15)
370  ! cc      PARAMETER (CRHL=0.70)
371  REAL t_coup
372  PARAMETER (t_coup=234.0)
373
374  ! Variables locales:
375  INTEGER i, k, kb, invb(klon)
376  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
377  REAL zdelta, zcor
378
379  ! Fonctions thermodynamiques:
380  include "FCTTRE.h"
381
382  ! Initialisation:
383
384  DO k = 1, klev
385    DO i = 1, klon
386      diafra(i, k) = 0.0
387      dialiq(i, k) = 0.0
388    END DO
389  END DO
390
391  DO i = 1, klon
392    invb(i) = klev
393    zdthmin(i) = 0.0
394  END DO
395
396  DO k = 2, klev/2 - 1
397    DO i = 1, klon
398      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
399        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
400      zdthdp = zdthdp*cloia
401      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
402        zdthmin(i) = zdthdp
403        invb(i) = k
404      END IF
405    END DO
406  END DO
407
408  DO i = 1, klon
409    kb = invb(i)
410    IF (thermcep) THEN
411      zdelta = max(0., sign(1.,rtt-t(i,kb)))
412      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
413      zqs = min(0.5, zqs)
414      zcor = 1./(1.-retv*zqs)
415      zqs = zqs*zcor
416    ELSE
417      IF (t(i,kb)<t_coup) THEN
418        zqs = qsats(t(i,kb))/pplay(i, kb)
419      ELSE
420        zqs = qsatl(t(i,kb))/pplay(i, kb)
421      END IF
422    END IF
423    zcll = cloib*zdthmin(i) + cloic
424    zcll = min(1.0, max(0.0,zcll))
425    zrhb = q(i, kb)/zqs
426    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
427    zcll = min(1.0, max(0.0,zcll))
428    diafra(i, kb) = max(diafra(i,kb), zcll)
429    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
430  END DO
431
432  RETURN
433END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.