source: LMDZ5/trunk/libf/phylmd/nuage.F90 @ 1998

Last change on this file since 1998 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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