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

Last change on this file was 5274, checked in by abarral, 9 hours ago

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