Context Navigation

newmicro.F @ 1226

Last change on this file since 1226 was 1220, checked in by lguez, 15 years ago

-- Replaced "integer*4" declarations by "integer", "real*8" by

"real(kind=8)" and "real*4" by "real". Note that these are the only
modifications in the files "radiation_AR4.F" and "sw_aeroAR4.F90".

-- Corrected the kind of arguments to "max" and "min".

-- Replaced "nH" edit descriptors, which is a deleted feature in

Fortran 95, by character strings.

-- "regr_lat_time_climoz" now allows the pressure coordinate in the

input file to be in descending order.

-- Replaced call to not standard function "float" by call to intrinsic

function "real".

-- Included file "radepsi.h" in "physiq" was not used. Removed it.

The following set of modifications is related to the management of time.

-- In "gcm", "leapfrog" and "sortvarc0", "day_ini" was defined as 1

plus number of days between the reference date "(annee_ref,
day_ref)" and the first day of the current simulation. Changed
definition: "(annee_ref, day_ini)" is the first day of the current
simulation. There is an accompanying modification for "day_end".

-- Corrected bug in call to "ioconf_startdate" in "gcm".

-- Added call to "ioconf_calendar" in "create_etat0_limit".

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 16.8 KB

Line
1	! $Id: newmicro.F 1220 2009-08-05 14:38:34Z jghattas $
2	!
3	SUBROUTINE newmicro (paprs, pplay,ok_newmicro,
4	. t, pqlwp, pclc, pcltau, pclemi,
5	. pch, pcl, pcm, pct, pctlwp,
6	s xflwp, xfiwp, xflwc, xfiwc,
7	e ok_aie,
8	e mass_solu_aero, mass_solu_aero_pi,
9	e bl95_b0, bl95_b1,
10	s cldtaupi, re, fl)
11	USE dimphy
12	IMPLICIT none
13	c======================================================================
14	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
15	c Objet: Calculer epaisseur optique et emmissivite des nuages
16	c======================================================================
17	c Arguments:
18	c t-------input-R-temperature
19	c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
20	c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
21	c
22	c ok_aie--input-L-apply aerosol indirect effect or not
23	c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
24	c mass_solu_aero_pi--input-R-dito, pre-industrial value
25	c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
26	c bl95_b1-input-R-a parameter, may be varied for tests ( -"- )
27	c
28	c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
29	c needed for the diagnostics of the aerosol indirect
30	c radiative forcing (see radlwsw)
31	c re------output-R-Cloud droplet effective radius multiplied by fl [um]
32	c fl------output-R-Denominator to re, introduced to avoid problems in
33	c the averaging of the output. fl is the fraction of liquid
34	c water clouds within a grid cell
35	c pcltau--output-R-epaisseur optique des nuages
36	c pclemi--output-R-emissivite des nuages (0 a 1)
37	c======================================================================
38	C
39	#include "YOMCST.h"
40	c
41	cym#include "dimensions.h"
42	cym#include "dimphy.h"
43	#include "nuage.h"
44	cIM cf. CR: include pour NOVLP et ZEPSEC
45	#include "radepsi.h"
46	#include "radopt.h"
47	REAL paprs(klon,klev+1), pplay(klon,klev)
48	REAL t(klon,klev)
49	c
50	REAL pclc(klon,klev)
51	REAL pqlwp(klon,klev)
52	REAL pcltau(klon,klev), pclemi(klon,klev)
53	c
54	REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
55	c
56	LOGICAL lo
57	c
58	REAL cetahb, cetamb
59	PARAMETER (cetahb = 0.45, cetamb = 0.80)
60	C
61	INTEGER i, k
62	cIM: 091003 REAL zflwp, zradef, zfice, zmsac
63	REAL zflwp(klon), zradef, zfice, zmsac
64	cIM: 091003 rajout
65	REAL xflwp(klon), xfiwp(klon)
66	REAL xflwc(klon,klev), xfiwc(klon,klev)
67	c
68	REAL radius, rad_chaud
69	cc PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
70	ccc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
71	c sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
72	REAL coef, coef_froi, coef_chau
73	PARAMETER (coef_chau=0.13, coef_froi=0.09)
74	REAL seuil_neb, t_glace
75	PARAMETER (seuil_neb=0.001, t_glace=273.0-15.0)
76	INTEGER nexpo ! exponentiel pour glace/eau
77	PARAMETER (nexpo=6)
78	ccc PARAMETER (nexpo=1)
79
80	c -- sb:
81	logical ok_newmicro
82	c parameter (ok_newmicro=.FALSE.)
83	cIM: 091003 real rel, tc, rei, zfiwp
84	real rel, tc, rei, zfiwp(klon)
85	real k_liq, k_ice0, k_ice, DF
86	parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g
87	parameter (DF=1.66) ! diffusivity factor
88	c sb --
89	cjq for the aerosol indirect effect
90	cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
91	cjq
92	LOGICAL ok_aie ! Apply AIE or not?
93	LOGICAL ok_a1lwpdep ! a1 LWP dependent?
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 within the grid cell)
103
104	REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
105
106	REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
107	cjq-end
108	cIM cf. CR:parametres supplementaires
109	REAL zclear(klon)
110	REAL zcloud(klon)
111
112	c **************************
113	c * *
114	c * DEBUT PARTIE OPTIMISEE *
115	c * *
116	c **************************
117
118	REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev)
119	REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var
120
121	c
122	c Calculer l'epaisseur optique et l'emmissivite des nuages
123	c
124	c IM inversion des DO
125	xflwp = 0.d0
126	xfiwp = 0.d0
127	xflwc = 0.d0
128	xfiwc = 0.d0
129
130	DO k = 1, klev
131	DO i = 1, klon
132	diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
133	ENDDO
134	ENDDO
135
136	IF (ok_newmicro) THEN
137
138
139	DO k = 1, klev
140	DO i = 1, klon
141	zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
142	zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
143	c IM Total Liquid/Ice water content
144	xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
145	xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
146	c IM In-Cloud Liquid/Ice water content
147	c xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
148	c xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
149	ENDDO
150	ENDDO
151
152	IF (ok_aie) THEN
153	DO k = 1, klev
154	DO i = 1, klon
155	! Formula "D" of Boucher and Lohmann, Tellus, 1995
156	!
157	cdnc(i,k) = 10.*(bl95_b0+bl95_b1
158	& log(MAX(mass_solu_aero(i,k),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	!
164	!
165	cdnc_pi(i,k) = 10.*(bl95_b0+bl95_b1
166	& log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
167	cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
168	ENDDO
169	ENDDO
170	DO k = 1, klev
171	DO i = 1, klon
172	! rad_chaud_tab(i,k) =
173	! & MAX(1.1e6
174	! & ((pqlwp(i,k)pplay(i,k)/(RD * T(i,k)))
175	! & /(4./3RPI1000.cdnc(i,k)) )*(1./3.),5.)
176	rad_chaud_tab(i,k) =
177	& 1.1
178	& ((pqlwp(i,k)pplay(i,k)/(RD * T(i,k)))
179	& /(4./3RPI1000.cdnc(i,k)) )*(1./3.)
180	rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.)
181	ENDDO
182	ENDDO
183	ELSE
184	DO k = 1, MIN(3,klev)
185	DO i = 1, klon
186	rad_chaud_tab(i,k) = rad_chau2
187	ENDDO
188	ENDDO
189	DO k = MIN(3,klev)+1, klev
190	DO i = 1, klon
191	rad_chaud_tab(i,k) = rad_chau1
192	ENDDO
193	ENDDO
194
195	ENDIF
196
197	DO k = 1, klev
198	! IF(.not.ok_aie) THEN
199	rad_chaud = rad_chau1
200	IF (k.LE.3) rad_chaud = rad_chau2
201	! ENDIF
202	DO i = 1, klon
203	IF (pclc(i,k) .LE. seuil_neb) THEN
204
205	c -- effective cloud droplet radius (microns):
206
207	c for liquid water clouds:
208	! For output diagnostics
209	!
210	! Cloud droplet effective radius [um]
211	!
212	! we multiply here with f * xl (fraction of liquid water
213	! clouds in the grid cell) to avoid problems in the
214	! averaging of the output.
215	! In the output of IOIPSL, derive the real cloud droplet
216	! effective radius as re/fl
217	!
218
219	fl(i,k) = seuil_neb*(1.-zfice2(i,k))
220	re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
221
222	pclc(i,k) = 0.0
223	pcltau(i,k) = 0.0
224	pclemi(i,k) = 0.0
225	cldtaupi(i,k) = 0.0
226	ELSE
227
228	c -- liquid/ice cloud water paths:
229
230	zflwp_var= 1000.(1.-zfice2(i,k))pqlwp(i,k)/pclc(i,k)
231	& *diff_paprs(i,k)
232	zfiwp_var= 1000.zfice2(i,k)pqlwp(i,k)/pclc(i,k)
233	& *diff_paprs(i,k)
234
235	c -- effective cloud droplet radius (microns):
236
237	c for liquid water clouds:
238
239	IF (ok_aie) THEN
240	radius =
241	& 1.1
242	& ((pqlwp(i,k)pplay(i,k)/(RD * T(i,k)))
243	& /(4./3.RPI1000.cdnc_pi(i,k)))*(1./3.)
244	radius = MAX(radius*1e6, 5.)
245
246	tc = t(i,k)-273.15
247	rei = 0.71*tc + 61.29
248	if (tc.le.-81.4) rei = 3.5
249	if (zflwp_var.eq.0.) radius = 1.
250	if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
251	cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
252	& + zfiwp_var * (3.448e-03 + 2.431/rei)
253	ENDIF ! ok_aie
254	! For output diagnostics
255	!
256	! Cloud droplet effective radius [um]
257	!
258	! we multiply here with f * xl (fraction of liquid water
259	! clouds in the grid cell) to avoid problems in the
260	! averaging of the output.
261	! In the output of IOIPSL, derive the real cloud droplet
262	! effective radius as re/fl
263	!
264
265	fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))
266	re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
267
268	rel = rad_chaud_tab(i,k)
269	c for ice clouds: as a function of the ambiant temperature
270	c [formula used by Iacobellis and Somerville (2000), with an
271	c asymptotical value of 3.5 microns at T<-81.4 C added to be
272	c consistent with observations of Heymsfield et al. 1986]:
273	tc = t(i,k)-273.15
274	rei = 0.71*tc + 61.29
275	if (tc.le.-81.4) rei = 3.5
276	c -- cloud optical thickness :
277
278	c [for liquid clouds, traditional formula,
279	c for ice clouds, Ebert & Curry (1992)]
280
281	if (zflwp_var.eq.0.) rel = 1.
282	if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
283	pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
284	& + zfiwp_var * (3.448e-03 + 2.431/rei)
285	c -- cloud infrared emissivity:
286
287	c [the broadband infrared absorption coefficient is parameterized
288	c as a function of the effective cld droplet radius]
289
290	c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
291	k_ice = k_ice0 + 1.0/rei
292
293	pclemi(i,k) = 1.0
294	& - EXP( -coef_chauzflwp_var - DFk_ice*zfiwp_var)
295
296	ENDIF
297
298	ENDDO
299	ENDDO
300
301	DO k = 1, klev
302	DO i = 1, klon
303	xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
304	xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
305	ENDDO
306	ENDDO
307
308	ELSE
309	DO k = 1, klev
310	rad_chaud = rad_chau1
311	IF (k.LE.3) rad_chaud = rad_chau2
312	DO i = 1, klon
313
314	IF (pclc(i,k) .LE. seuil_neb) THEN
315
316	pclc(i,k) = 0.0
317	pcltau(i,k) = 0.0
318	pclemi(i,k) = 0.0
319	cldtaupi(i,k) = 0.0
320
321	ELSE
322
323	zflwp_var = 1000.pqlwp(i,k)diff_paprs(i,k)
324	& /pclc(i,k)
325
326	zfice1 = MIN(
327	& MAX( 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
328	& ,0.0),1.0)**nexpo
329
330	radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
331	coef = coef_chau * (1.-zfice1) + coef_froi * zfice1
332
333	pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
334	pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
335
336	ENDIF
337
338	ENDDO
339	ENDDO
340	ENDIF
341
342	IF (.NOT.ok_aie) THEN
343	DO k = 1, klev
344	DO i = 1, klon
345	cldtaupi(i,k)=pcltau(i,k)
346	ENDDO
347	ENDDO
348	ENDIF
349
350	ccc DO k = 1, klev
351	ccc DO i = 1, klon
352	ccc t(i,k) = t(i,k)
353	ccc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
354	ccc lo = pclc(i,k) .GT. (2.*1.e-5)
355	ccc zflwp = pqlwp(i,k)1000.(paprs(i,k)-paprs(i,k+1))
356	ccc . /(rg*pclc(i,k))
357	ccc zradef = 10.0 + (1.-sigs(k))*45.0
358	ccc pcltau(i,k) = 1.5 * zflwp / zradef
359	ccc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
360	ccc zmsac = 0.13(1.0-zfice) + 0.08zfice
361	ccc pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
362	ccc if (.NOT.lo) pclc(i,k) = 0.0
363	ccc if (.NOT.lo) pcltau(i,k) = 0.0
364	ccc if (.NOT.lo) pclemi(i,k) = 0.0
365	ccc ENDDO
366	ccc ENDDO
367	ccccc print*, 'pas de nuage dans le rayonnement'
368	ccccc DO k = 1, klev
369	ccccc DO i = 1, klon
370	ccccc pclc(i,k) = 0.0
371	ccccc pcltau(i,k) = 0.0
372	ccccc pclemi(i,k) = 0.0
373	ccccc ENDDO
374	ccccc ENDDO
375	C
376	C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
377	C
378	c IM cf. CR:test: calcul prenant ou non en compte le recouvrement
379	c initialisations
380	DO i=1,klon
381	zclear(i)=1.
382	zcloud(i)=0.
383	pch(i)=1.0
384	pcm(i) = 1.0
385	pcl(i) = 1.0
386	pctlwp(i) = 0.0
387	ENDDO
388	C
389	cIM cf CR DO k=1,klev
390	DO k = klev, 1, -1
391	DO i = 1, klon
392	pctlwp(i) = pctlwp(i)
393	& + pqlwp(i,k)*diff_paprs(i,k)
394	ENDDO
395	ENDDO
396	c IM cf. CR
397	IF (NOVLP.EQ.1) THEN
398	DO k = klev, 1, -1
399	DO i = 1, klon
400	zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
401	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
402	pct(i)=1.-zclear(i)
403	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
404	pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
405	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
406	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
407	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
408	pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
409	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
410	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
411	pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
412	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
413	endif
414	zcloud(i)=pclc(i,k)
415	ENDDO
416	ENDDO
417	ELSE IF (NOVLP.EQ.2) THEN
418	DO k = klev, 1, -1
419	DO i = 1, klon
420	zcloud(i)=MAX(pclc(i,k),zcloud(i))
421	pct(i)=zcloud(i)
422	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
423	pch(i) = MIN(pclc(i,k),pch(i))
424	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
425	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
426	pcm(i) = MIN(pclc(i,k),pcm(i))
427	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
428	pcl(i) = MIN(pclc(i,k),pcl(i))
429	endif
430	ENDDO
431	ENDDO
432	ELSE IF (NOVLP.EQ.3) THEN
433	DO k = klev, 1, -1
434	DO i = 1, klon
435	zclear(i)=zclear(i)*(1.-pclc(i,k))
436	pct(i)=1-zclear(i)
437	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
438	pch(i) = pch(i)*(1.0-pclc(i,k))
439	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
440	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
441	pcm(i) = pcm(i)*(1.0-pclc(i,k))
442	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
443	pcl(i) = pcl(i)*(1.0-pclc(i,k))
444	endif
445	ENDDO
446	ENDDO
447	ENDIF
448
449	C
450	DO i = 1, klon
451	c IM cf. CR pct(i)=1.-pct(i)
452	pch(i)=1.-pch(i)
453	pcm(i)=1.-pcm(i)
454	pcl(i)=1.-pcl(i)
455	ENDDO
456
457	C
458	RETURN
459	END

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

Download in other formats:

Original Format