source: LMDZ6/branches/Amaury_dev/libf/phylmd/nuage.F90 @ 5137

Last change on this file since 5137 was 5137, checked in by abarral, 8 weeks ago

Put gradsdef.h, tracstoke.h, clesphys.h into 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.6 KB
Line 
1! $Id: nuage.F90 5137 2024-07-28 20:25:12Z 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  USE lmdz_clesphys
12
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 "YOMCST.h"
43  include "nuage.h" ! JBM 3/14
44
45  REAL paprs(klon, klev+1), pplay(klon, klev)
46  REAL t(klon, klev)
47
48  REAL pclc(klon, klev)
49  REAL pqlwp(klon, klev), picefra(klon,klev)
50  REAL pcltau(klon, klev), pclemi(klon, klev)
51
52  REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
53  REAL distcltop(klon,klev)
54  REAL temp_cltop(klon,klev)
55  LOGICAL lo
56
57  REAL cetahb, cetamb
58  PARAMETER (cetahb=0.45, cetamb=0.80)
59
60  INTEGER i, k
61  REAL zflwp, zradef, zfice(klon), zmsac
62
63  REAL radius, rad_chaud
64! JBM (3/14) parameters already defined in nuage.h:
65! REAL rad_froid, rad_chau1, rad_chau2
66! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
67  ! cc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
68  ! sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
69  REAL coef, coef_froi, coef_chau
70  PARAMETER (coef_chau=0.13, coef_froi=0.09)
71  REAL seuil_neb
72  PARAMETER (seuil_neb=0.001)
73! JBM (3/14) nexpo is replaced by exposant_glace
74! REAL nexpo ! exponentiel pour glace/eau
75! PARAMETER (nexpo=6.)
76  REAL, PARAMETER :: t_glace_min_old = 258.
77  INTEGER, PARAMETER :: exposant_glace_old = 6
78
79
80  ! jq for the aerosol indirect effect
81  ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
82  ! jq
83  LOGICAL ok_aie ! Apply AIE or not?
84
85  REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols[ug m-3]
86  REAL mass_solu_aero_pi(klon, klev) ! - " - pre-industrial value
87  REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
88  REAL re(klon, klev) ! cloud droplet effective radius [um]
89  REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
90  REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
91
92  REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds
93  ! within the grid cell)
94
95  REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
96
97  REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
98  REAl dzfice(klon)
99  ! jq-end
100
101  ! cc      PARAMETER (nexpo=1)
102
103  ! Calculer l'epaisseur optique et l'emmissivite des nuages
104
105  DO k = 1, klev
106     IF (iflag_t_glace==0) THEN
107       DO i = 1, klon
108        zfice(i) = 1.0 - (t(i,k)-t_glace_min_old)/(273.13-t_glace_min_old)
109        zfice(i) = min(max(zfice(i),0.0), 1.0)
110        zfice(i) = zfice(i)**exposant_glace_old
111       ENDDO
112     ELSE ! of IF (iflag_t_glace.EQ.0)
113! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod
114!       zfice(i) = icefrac_lsc(t(i,k), t_glace_min, &
115!                           t_glace_max, exposant_glace)
116        IF (ok_new_lscp) THEN
117            CALL icefrac_lscp(klon,t(:,k),iflag_ice_thermo,distcltop(:,k),temp_cltop(:,k),zfice(:),dzfice(:))
118        ELSE
119            CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:))
120
121        ENDIF
122
123    IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
124    ! EV: take the ice fraction directly from the lscp code
125    ! consistent only for non convective grid points
126    ! critical for mixed phase clouds
127        DO i=1,klon
128        IF (.NOT. ptconv(i,k)) THEN
129           zfice(i)=picefra(i,k)
130        ENDIF
131        ENDDO
132    ENDIF
133
134
135     ENDIF
136
137    DO i = 1, klon
138      rad_chaud = rad_chau1
139      IF (k<=3) rad_chaud = rad_chau2
140
141      pclc(i, k) = max(pclc(i,k), seuil_neb)
142      zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1))
143
144      IF (ok_aie) THEN
145          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
146
147        cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
148          1.E-4))/log(10.))*1.E6 !-m-3
149          ! Cloud droplet number concentration (CDNC) is restricted
150          ! to be within [20, 1000 cm^3]
151
152        cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
153        cdnc_pi(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero_pi(i,k), &
154          1.E-4))/log(10.))*1.E6 !-m-3
155        cdnc_pi(i, k) = min(1000.E6, max(20.E6,cdnc_pi(i,k)))
156
157
158          ! air density: pplay(i,k) / (RD * zT(i,k))
159          ! factor 1.1: derive effective radius from volume-mean radius
160          ! factor 1000 is the water density
161          ! _chaud means that this is the CDR for liquid water clouds
162
163        rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. &
164          *cdnc(i,k)))**(1./3.)
165
166          ! Convert to um. CDR shall be at least 3 um.
167
168        rad_chaud = max(rad_chaud*1.E6, 3.)
169
170          ! For output diagnostics
171
172          ! Cloud droplet effective radius [um]
173
174          ! we multiply here with f * xl (fraction of liquid water
175          ! clouds in the grid cell) to avoid problems in the
176          ! averaging of the output.
177          ! In the output of IOIPSL, derive the real cloud droplet
178          ! effective radius as re/fl
179
180        fl(i, k) = pclc(i, k)*(1.-zfice(i))
181        re(i, k) = rad_chaud*fl(i, k)
182
183          ! Pre-industrial cloud opt thickness
184
185          ! "radius" is calculated as rad_chaud above (plus the
186          ! ice cloud contribution) but using cdnc_pi instead of
187          ! cdnc.
188        radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* &
189          1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i)
190        cldtaupi(i, k) = 3.0/2.0*zflwp/radius
191      END IF ! ok_aie
192
193      radius = rad_chaud*(1.-zfice(i)) + rad_froid*zfice(i)
194      coef = coef_chau*(1.-zfice(i)) + coef_froi*zfice(i)
195      pcltau(i, k) = 3.0/2.0*zflwp/radius
196      pclemi(i, k) = 1.0 - exp(-coef*zflwp)
197      lo = (pclc(i,k)<=seuil_neb)
198      IF (lo) pclc(i, k) = 0.0
199      IF (lo) pcltau(i, k) = 0.0
200      IF (lo) pclemi(i, k) = 0.0
201
202      IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k)
203    END DO
204  END DO
205  ! cc      DO k = 1, klev
206  ! cc      DO i = 1, klon
207  ! cc         t(i,k) = t(i,k)
208  ! cc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
209  ! cc         lo = pclc(i,k) .GT. (2.*1.e-5)
210  ! cc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
211  ! cc     .          /(rg*pclc(i,k))
212  ! cc         zradef = 10.0 + (1.-sigs(k))*45.0
213  ! cc         pcltau(i,k) = 1.5 * zflwp / zradef
214  ! cc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
215  ! cc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
216  ! cc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
217  ! cc         if (.NOT.lo) pclc(i,k) = 0.0
218  ! cc         if (.NOT.lo) pcltau(i,k) = 0.0
219  ! cc         if (.NOT.lo) pclemi(i,k) = 0.0
220  ! cc      ENDDO
221  ! cc      ENDDO
222  ! ccccc      PRINT*, 'pas de nuage dans le rayonnement'
223  ! ccccc      DO k = 1, klev
224  ! ccccc      DO i = 1, klon
225  ! ccccc         pclc(i,k) = 0.0
226  ! ccccc         pcltau(i,k) = 0.0
227  ! ccccc         pclemi(i,k) = 0.0
228  ! ccccc      ENDDO
229  ! ccccc      ENDDO
230
231  ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
232
233  DO i = 1, klon
234    pct(i) = 1.0
235    pch(i) = 1.0
236    pcm(i) = 1.0
237    pcl(i) = 1.0
238    pctlwp(i) = 0.0
239  END DO
240
241  DO k = klev, 1, -1
242    DO i = 1, klon
243      pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg
244      pct(i) = pct(i)*(1.0-pclc(i,k))
245      IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k))
246      IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) &
247        pcm(i) = pcm(i)*(1.0-pclc(i,k))
248      IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k))
249    END DO
250  END DO
251
252  DO i = 1, klon
253    pct(i) = 1. - pct(i)
254    pch(i) = 1. - pch(i)
255    pcm(i) = 1. - pcm(i)
256    pcl(i) = 1. - pcl(i)
257  END DO
258
259
260END SUBROUTINE nuage
261SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq)
262  USE dimphy
263  IMPLICIT 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  include "YOMCST.h"
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
342END SUBROUTINE diagcld1
343SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq)
344  USE dimphy
345  IMPLICIT NONE
346
347  include "YOMCST.h"
348
349  ! Arguments d'entree:
350  REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche
351  REAL pplay(klon, klev) ! pression (Pa) au milieu de couche
352  REAL t(klon, klev) ! temperature (K)
353  REAL q(klon, klev) ! humidite specifique (Kg/Kg)
354
355  ! Arguments de sortie:
356  REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee
357  REAL dialiq(klon, klev) ! eau liquide nuageuse
358
359  REAL cetamb
360  PARAMETER (cetamb=0.80)
361  REAL cloia, cloib, cloic, cloid
362  PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0)
363  ! cc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
364  REAL rgammas
365  PARAMETER (rgammas=0.05)
366  REAL crhl
367  PARAMETER (crhl=0.15)
368  ! cc      PARAMETER (CRHL=0.70)
369  REAL t_coup
370  PARAMETER (t_coup=234.0)
371
372  ! Variables locales:
373  INTEGER i, k, kb, invb(klon)
374  REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
375  REAL zdelta, zcor
376
377  ! Fonctions thermodynamiques:
378  include "YOETHF.h"
379  include "FCTTRE.h"
380
381  ! Initialisation:
382
383  DO k = 1, klev
384    DO i = 1, klon
385      diafra(i, k) = 0.0
386      dialiq(i, k) = 0.0
387    END DO
388  END DO
389
390  DO i = 1, klon
391    invb(i) = klev
392    zdthmin(i) = 0.0
393  END DO
394
395  DO k = 2, klev/2 - 1
396    DO i = 1, klon
397      zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - &
398        rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1)
399      zdthdp = zdthdp*cloia
400      IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp<zdthmin(i)) THEN
401        zdthmin(i) = zdthdp
402        invb(i) = k
403      END IF
404    END DO
405  END DO
406
407  DO i = 1, klon
408    kb = invb(i)
409    IF (thermcep) THEN
410      zdelta = max(0., sign(1.,rtt-t(i,kb)))
411      zqs = r2es*foeew(t(i,kb), zdelta)/pplay(i, kb)
412      zqs = min(0.5, zqs)
413      zcor = 1./(1.-retv*zqs)
414      zqs = zqs*zcor
415    ELSE
416      IF (t(i,kb)<t_coup) THEN
417        zqs = qsats(t(i,kb))/pplay(i, kb)
418      ELSE
419        zqs = qsatl(t(i,kb))/pplay(i, kb)
420      END IF
421    END IF
422    zcll = cloib*zdthmin(i) + cloic
423    zcll = min(1.0, max(0.0,zcll))
424    zrhb = q(i, kb)/zqs
425    IF (zcll>0.0 .AND. zrhb<crhl) zcll = zcll*(1.-(crhl-zrhb)*cloid)
426    zcll = min(1.0, max(0.0,zcll))
427    diafra(i, kb) = max(diafra(i,kb), zcll)
428    dialiq(i, kb) = diafra(i, kb)*rgammas*zqs
429  END DO
430
431
432END SUBROUTINE diagcld2
Note: See TracBrowser for help on using the repository browser.