Context Navigation

nuage.F @ 416

Last change on this file since 416 was 44, checked in by lmdz, 25 years ago
Ajustement de parametres, separation de la routine "nuages diagnostiques" en deux routines (une appelee uniquement dans le cas Tiedtke, l'autre dans le cas ok-stratus) L.Li 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
5	c======================================================================
6	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
7	c Objet: Calculer epaisseur optique et emmissivite des nuages
8	c======================================================================
9	c Arguments:
10	c t-------input-R-temperature
11	c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
12	c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
13	c
14	c pcltau--output-R-epaisseur optique des nuages
15	c pclemi--output-R-emissivite des nuages (0 a 1)
16	c======================================================================
17	C
18	#include "YOMCST.h"
19	c
20	#include "dimensions.h"
21	#include "dimphy.h"
22	REAL paprs(klon,klev+1), pplay(klon,klev)
23	REAL t(klon,klev)
24	c
25	REAL pclc(klon,klev)
26	REAL pqlwp(klon,klev)
27	REAL pcltau(klon,klev), pclemi(klon,klev)
28	c
29	REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
30	c
31	LOGICAL lo
32	c
33	REAL cetahb, cetamb
34	PARAMETER (cetahb = 0.45, cetamb = 0.80)
35	C
36	INTEGER i, k
37	REAL zflwp, zradef, zfice, zmsac
38	c
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)
41	ccc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
42	c 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)
49	ccc PARAMETER (nexpo=1)
50	c
51	c Calculer l'epaisseur optique et l'emmissivite des nuages
52	c
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
73	ccc DO k = 1, klev
74	ccc DO i = 1, klon
75	ccc t(i,k) = t(i,k)
76	ccc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
77	ccc lo = pclc(i,k) .GT. (2.*1.e-5)
78	ccc zflwp = pqlwp(i,k)1000.(paprs(i,k)-paprs(i,k+1))
79	ccc . /(rg*pclc(i,k))
80	ccc zradef = 10.0 + (1.-sigs(k))*45.0
81	ccc pcltau(i,k) = 1.5 * zflwp / zradef
82	ccc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
83	ccc zmsac = 0.13(1.0-zfice) + 0.08zfice
84	ccc pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
85	ccc if (.NOT.lo) pclc(i,k) = 0.0
86	ccc if (.NOT.lo) pcltau(i,k) = 0.0
87	ccc if (.NOT.lo) pclemi(i,k) = 0.0
88	ccc ENDDO
89	ccc ENDDO
90	cccccc print*, 'pas de nuage dans le rayonnement'
91	cccccc DO k = 1, klev
92	cccccc DO i = 1, klon
93	cccccc pclc(i,k) = 0.0
94	cccccc pcltau(i,k) = 0.0
95	cccccc pclemi(i,k) = 0.0
96	cccccc ENDDO
97	cccccc ENDDO
98	C
99	C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
100	C
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
108	C
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
123	C
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
130	C
131	RETURN
132	END
133	SUBROUTINE diagcld1(paprs,pplay,rain,snow,kbot,ktop,
134	. diafra,dialiq)
135	IMPLICIT none
136	c
137	c Laurent Li (LMD/CNRS), le 12 octobre 1998
138	c (adaptation du code ECMWF)
139	c
140	c Dans certains cas, le schema pronostique des nuages n'est
141	c pas suffisament performant. On a donc besoin de diagnostiquer
142	c ces nuages. Je dois avouer que c'est une frustration.
143	c
144	#include "dimensions.h"
145	#include "dimphy.h"
146	#include "YOMCST.h"
147	c
148	c 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
157	c
158	c Arguments de sortie:
159	REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
160	REAL dialiq(klon,klev) ! eau liquide nuageuse
161	c
162	c 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)
176	c
177	c Variables locales:
178	INTEGER i, k
179	REAL zcc(klon)
180	c
181	c Initialisation:
182	c
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
189	c
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
197	c
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
206	c
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
215	c
216	RETURN
217	END
218	SUBROUTINE diagcld2(paprs,pplay,t,q, diafra,dialiq)
219	IMPLICIT none
220	c
221	#include "dimensions.h"
222	#include "dimphy.h"
223	#include "YOMCST.h"
224	c
225	c 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)
230	c
231	c Arguments de sortie:
232	REAL diafra(klon,klev) ! fraction nuageuse diagnostiquee
233	REAL dialiq(klon,klev) ! eau liquide nuageuse
234	c
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)
239	ccc 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)
244	ccc PARAMETER (CRHL=0.70)
245	REAL t_coup
246	PARAMETER (t_coup=234.0)
247	c
248	c Variables locales:
249	INTEGER i, k, kb, invb(klon)
250	REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp
251	REAL zdelta, zcor
252	c
253	c Fonctions thermodynamiques:
254	#include "YOETHF.h"
255	#include "FCTTRE.h"
256	c
257	c Initialisation:
258	c
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
265	c
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
308	c
309	RETURN
310	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format