source: LMDZ.3.3/branches/rel-LF/libf/phylmd/nuage.F @ 158

Last change on this file since 158 was 79, checked in by (none), 25 years ago

This commit was manufactured by cvs2svn to create branch 'rel-LF'.

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.3 KB
Line 
1      SUBROUTINE nuage (paprs, pplay,
2     .                  t, pqlwp, pclc, pcltau, pclemi,
3     .                  pch, pcl, pcm, pct, pctlwp)
4      IMPLICIT none
5c======================================================================
6c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
7c Objet: Calculer epaisseur optique et emmissivite des nuages
8c======================================================================
9c Arguments:
10c t-------input-R-temperature
11c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
12c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
13c
14c pcltau--output-R-epaisseur optique des nuages
15c pclemi--output-R-emissivite des nuages (0 a 1)
16c======================================================================
17C
18#include "YOMCST.h"
19c
20#include "dimensions.h"
21#include "dimphy.h"
22      REAL paprs(klon,klev+1), pplay(klon,klev)
23      REAL t(klon,klev)
24c
25      REAL pclc(klon,klev)
26      REAL pqlwp(klon,klev)
27      REAL pcltau(klon,klev), pclemi(klon,klev)
28c
29      REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
30c
31      LOGICAL lo
32c
33      REAL cetahb, cetamb
34      PARAMETER (cetahb = 0.45, cetamb = 0.80)
35C
36      INTEGER i, k
37      REAL zflwp, zradef, zfice, zmsac
38c
39      REAL radius, rad_froid, rad_chaud, rad_chau1, rad_chau2
40      PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
41ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
42c sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
43      REAL coef, coef_froi, coef_chau
44      PARAMETER (coef_chau=0.13, coef_froi=0.09)
45      REAL seuil_neb, t_glace
46      PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
47      INTEGER nexpo ! exponentiel pour glace/eau
48      PARAMETER (nexpo=6)
49ccc      PARAMETER (nexpo=1)
50c
51c Calculer l'epaisseur optique et l'emmissivite des nuages
52c
53      DO k = 1, klev
54      DO i = 1, klon
55         rad_chaud = rad_chau1
56         IF (k.LE.3) rad_chaud = rad_chau2
57         pclc(i,k) = MAX(pclc(i,k), seuil_neb)
58         zflwp = 1000.*pqlwp(i,k)/RG/pclc(i,k)
59     .          *(paprs(i,k)-paprs(i,k+1))
60         zfice = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
61         zfice = MIN(MAX(zfice,0.0),1.0)
62         zfice = zfice**nexpo
63         radius = rad_chaud * (1.-zfice) + rad_froid * zfice
64         coef = coef_chau * (1.-zfice) + coef_froi * zfice
65         pcltau(i,k) = 3.0/2.0 * zflwp / radius
66         pclemi(i,k) = 1.0 - EXP( - coef * zflwp)
67         lo = (pclc(i,k) .LE. seuil_neb)
68         IF (lo) pclc(i,k) = 0.0
69         IF (lo) pcltau(i,k) = 0.0
70         IF (lo) pclemi(i,k) = 0.0
71      ENDDO
72      ENDDO
73ccc      DO k = 1, klev
74ccc      DO i = 1, klon
75ccc         t(i,k) = t(i,k)
76ccc         pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
77ccc         lo = pclc(i,k) .GT. (2.*1.e-5)
78ccc         zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
79ccc     .          /(rg*pclc(i,k))
80ccc         zradef = 10.0 + (1.-sigs(k))*45.0
81ccc         pcltau(i,k) = 1.5 * zflwp / zradef
82ccc         zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
83ccc         zmsac = 0.13*(1.0-zfice) + 0.08*zfice
84ccc         pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
85ccc         if (.NOT.lo) pclc(i,k) = 0.0
86ccc         if (.NOT.lo) pcltau(i,k) = 0.0
87ccc         if (.NOT.lo) pclemi(i,k) = 0.0
88ccc      ENDDO
89ccc      ENDDO
90cccccc      print*, 'pas de nuage dans le rayonnement'
91cccccc      DO k = 1, klev
92cccccc      DO i = 1, klon
93cccccc         pclc(i,k) = 0.0
94cccccc         pcltau(i,k) = 0.0
95cccccc         pclemi(i,k) = 0.0
96cccccc      ENDDO
97cccccc      ENDDO
98C
99C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
100C
101      DO i = 1, klon
102         pct(i)=1.0
103         pch(i)=1.0
104         pcm(i) = 1.0
105         pcl(i) = 1.0
106         pctlwp(i) = 0.0
107      ENDDO
108C
109      DO k = klev, 1, -1
110      DO i = 1, klon
111         pctlwp(i) = pctlwp(i)
112     .             + pqlwp(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
113         pct(i) = pct(i)*(1.0-pclc(i,k))
114         if (pplay(i,k).LE.cetahb*paprs(i,1))
115     .      pch(i) = pch(i)*(1.0-pclc(i,k))
116         if (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
117     .       pplay(i,k).LE.cetamb*paprs(i,1))
118     .      pcm(i) = pcm(i)*(1.0-pclc(i,k))
119         if (pplay(i,k).GT.cetamb*paprs(i,1))
120     .      pcl(i) = pcl(i)*(1.0-pclc(i,k))
121      ENDDO
122      ENDDO
123C
124      DO i = 1, klon
125         pct(i)=1.-pct(i)
126         pch(i)=1.-pch(i)
127         pcm(i)=1.-pcm(i)
128         pcl(i)=1.-pcl(i)
129      ENDDO
130C
131      RETURN
132      END
133      SUBROUTINE diagcld1(paprs,pplay,rain,snow,kbot,ktop,
134     .                   diafra,dialiq)
135      IMPLICIT none
136c
137c Laurent Li (LMD/CNRS), le 12 octobre 1998
138c                        (adaptation du code ECMWF)
139c
140c Dans certains cas, le schema pronostique des nuages n'est
141c pas suffisament performant. On a donc besoin de diagnostiquer
142c ces nuages. Je dois avouer que c'est une frustration.
143c
144#include "dimensions.h"
145#include "dimphy.h"
146#include "YOMCST.h"
147c
148c Arguments d'entree:
149      REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
150      REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
151      REAL t(klon,klev) ! temperature (K)
152      REAL q(klon,klev) ! humidite specifique (Kg/Kg)
153      REAL rain(klon) ! pluie convective (kg/m2/s)
154      REAL snow(klon) ! neige convective (kg/m2/s)
155      INTEGER ktop(klon) ! sommet de la convection
156      INTEGER kbot(klon) ! bas de la convection
157c
158c Arguments de sortie:
159      REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
160      REAL dialiq(klon,klev) ! eau liquide nuageuse
161c
162c Constantes ajustables:
163      REAL CANVA, CANVB, CANVH
164      PARAMETER (CANVA=2.0, CANVB=0.3, CANVH=0.4)
165      REAL CCA, CCB, CCC
166      PARAMETER (CCA=0.125, CCB=1.5, CCC=0.8)
167      REAL CCFCT, CCSCAL
168      PARAMETER (CCFCT=0.400)
169      PARAMETER (CCSCAL=1.0E+11)
170      REAL CETAHB, CETAMB
171      PARAMETER (CETAHB=0.45, CETAMB=0.80)
172      REAL CCLWMR
173      PARAMETER (CCLWMR=1.E-04)
174      REAL ZEPSCR
175      PARAMETER (ZEPSCR=1.0E-10)
176c
177c Variables locales:
178      INTEGER i, k
179      REAL zcc(klon)
180c
181c Initialisation:
182c
183      DO k = 1, klev
184      DO i = 1, klon
185         diafra(i,k) = 0.0
186         dialiq(i,k) = 0.0
187      ENDDO
188      ENDDO
189c
190      DO i = 1, klon ! Calculer la fraction nuageuse
191      zcc(i) = 0.0
192      IF((rain(i)+snow(i)).GT.0.) THEN
193         zcc(i)= CCA * LOG(MAX(ZEPSCR,(rain(i)+snow(i))*CCSCAL))-CCB
194         zcc(i)= MIN(CCC,MAX(0.0,zcc(i)))
195      ENDIF
196      ENDDO
197c
198      DO i = 1, klon ! pour traiter les enclumes
199      diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),zcc(i)*CCFCT)
200      IF ((zcc(i).GE.CANVH) .AND.
201     .    (pplay(i,ktop(i)).LE.CETAHB*paprs(i,1)))
202     . diafra(i,ktop(i)) = MAX(diafra(i,ktop(i)),
203     .                         MAX(zcc(i)*CCFCT,CANVA*(zcc(i)-CANVB)))
204      dialiq(i,ktop(i))=CCLWMR*diafra(i,ktop(i))
205      ENDDO
206c
207      DO k = 1, klev ! nuages convectifs (sauf enclumes)
208      DO i = 1, klon
209      IF (k.LT.ktop(i) .AND. k.GE.kbot(i)) THEN
210         diafra(i,k)=MAX(diafra(i,k),zcc(i)*CCFCT)
211         dialiq(i,k)=CCLWMR*diafra(i,k)
212      ENDIF
213      ENDDO
214      ENDDO
215c
216      RETURN
217      END
218      SUBROUTINE diagcld2(paprs,pplay,t,q, diafra,dialiq)
219      IMPLICIT none
220c
221#include "dimensions.h"
222#include "dimphy.h"
223#include "YOMCST.h"
224c
225c Arguments d'entree:
226      REAL paprs(klon,klev+1) ! pression (Pa) a inter-couche
227      REAL pplay(klon,klev) ! pression (Pa) au milieu de couche
228      REAL t(klon,klev) ! temperature (K)
229      REAL q(klon,klev) ! humidite specifique (Kg/Kg)
230c
231c Arguments de sortie:
232      REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
233      REAL dialiq(klon,klev) ! eau liquide nuageuse
234c
235      REAL CETAMB
236      PARAMETER (CETAMB=0.80)
237      REAL CLOIA, CLOIB, CLOIC, CLOID
238      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.6, CLOID=5.0)
239ccc      PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0)
240      REAL RGAMMAS
241      PARAMETER (RGAMMAS=0.05)
242      REAL CRHL
243      PARAMETER (CRHL=0.15)
244ccc      PARAMETER (CRHL=0.70)
245      REAL t_coup
246      PARAMETER (t_coup=234.0)
247c
248c Variables locales:
249      INTEGER i, k, kb, invb(klon)
250      REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
251      REAL zdelta, zcor
252c
253c Fonctions thermodynamiques:
254#include "YOETHF.h"
255#include "FCTTRE.h"
256c
257c Initialisation:
258c
259      DO k = 1, klev
260      DO i = 1, klon
261         diafra(i,k) = 0.0
262         dialiq(i,k) = 0.0
263      ENDDO
264      ENDDO
265c
266      DO i = 1, klon
267         invb(i) = klev
268         zdthmin(i)=0.0
269      ENDDO
270
271      DO k = 2, klev/2-1
272      DO i = 1, klon
273         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
274     .          - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/paprs(i,k+1)
275         zdthdp = zdthdp * CLOIA
276         IF (pplay(i,k).GT.CETAMB*paprs(i,1) .AND.
277     .       zdthdp.LT.zdthmin(i) ) THEN
278            zdthmin(i) = zdthdp
279            invb(i) = k
280         ENDIF
281      ENDDO
282      ENDDO
283
284      DO i = 1, klon
285         kb=invb(i)
286         IF (thermcep) THEN
287           zdelta=MAX(0.,SIGN(1.,RTT-t(i,kb)))
288           zqs= R2ES*FOEEW(t(i,kb),zdelta)/pplay(i,kb)
289           zqs=MIN(0.5,zqs)
290           zcor=1./(1.-RETV*zqs)
291           zqs=zqs*zcor
292         ELSE
293           IF (t(i,kb) .LT. t_coup) THEN
294              zqs = qsats(t(i,kb)) / pplay(i,kb)
295           ELSE
296              zqs = qsatl(t(i,kb)) / pplay(i,kb)
297           ENDIF
298         ENDIF
299         zcll = CLOIB * zdthmin(i) + CLOIC
300         zcll = MIN(1.0,MAX(0.0,zcll))
301         zrhb= q(i,kb)/zqs
302         IF (zcll.GT.0.0.AND.zrhb.LT.CRHL)
303     .   zcll=zcll*(1.-(CRHL-zrhb)*CLOID)
304         zcll=MIN(1.0,MAX(0.0,zcll))
305         diafra(i,kb) = MAX(diafra(i,kb),zcll)
306         dialiq(i,kb)= diafra(i,kb) * RGAMMAS*zqs
307      ENDDO
308c
309      RETURN
310      END
Note: See TracBrowser for help on using the repository browser.