source: LMDZ6/trunk/libf/phylmd/nuage.F90 @ 3592

Last change on this file since 3592 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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