Context Navigation

newmicro.F @ 1292

Last change on this file since 1292 was 1286, checked in by Laurent Fairhead, 15 years ago
Les parametres definissant la transition eau glace dans les nuages sont maintenant lus dans physiq.def. Les valeurs par defaut donnent les memes resultats que precedemment. JLD
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.3 KB

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

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

Download in other formats:

Original Format