Context Navigation

newmicro.F @ 1235

Last change on this file since 1235 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

Rev	Line
[1179]	1	! $Id: newmicro.F 1220 2009-08-05 14:38:34Z musat $
	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,
[1183]	8	e mass_solu_aero, mass_solu_aero_pi,
[524]	9	e bl95_b0, bl95_b1,
	10	s cldtaupi, re, fl)
[766]	11	USE dimphy
[524]	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
[1183]	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
[524]	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
[766]	41	cym#include "dimensions.h"
	42	cym#include "dimphy.h"
[524]	43	#include "nuage.h"
[685]	44	cIM cf. CR: include pour NOVLP et ZEPSEC
	45	#include "radepsi.h"
	46	#include "radopt.h"
[524]	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
[1183]	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)
[524]	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
[685]	108	cIM cf. CR:parametres supplementaires
	109	REAL zclear(klon)
	110	REAL zcloud(klon)
[1092]	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
[524]	121	c
	122	c Calculer l'epaisseur optique et l'emmissivite des nuages
	123	c
[1092]	124	c IM inversion des DO
	125	xflwp = 0.d0
	126	xfiwp = 0.d0
	127	xflwc = 0.d0
	128	xfiwc = 0.d0
	129
[524]	130	DO k = 1, klev
[1092]	131	DO i = 1, klon
	132	diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
	133	ENDDO
	134	ENDDO
[524]	135
[1092]	136	IF (ok_newmicro) THEN
[524]	137
	138
[1092]	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
[524]	151
[1092]	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
[1183]	158	& log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
[1092]	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
[1183]	166	& log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
[1092]	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
[524]	194
[1092]	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
[524]	227
[1092]	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)
[524]	295
[1092]	296	ENDIF
	297
	298	ENDDO
	299	ENDDO
[524]	300
[1092]	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
[524]	307
[1092]	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
[524]	315
[1092]	316	pclc(i,k) = 0.0
	317	pcltau(i,k) = 0.0
	318	pclemi(i,k) = 0.0
	319	cldtaupi(i,k) = 0.0
[524]	320
[1092]	321	ELSE
[524]	322
[1092]	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
[524]	332
[1092]	333	pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
	334	pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
[524]	335
[1092]	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
[524]	349
[1092]	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
[685]	380	DO i=1,klon
	381	zclear(i)=1.
	382	zcloud(i)=0.
[524]	383	pch(i)=1.0
	384	pcm(i) = 1.0
	385	pcl(i) = 1.0
	386	pctlwp(i) = 0.0
	387	ENDDO
	388	C
[685]	389	cIM cf CR DO k=1,klev
[524]	390	DO k = klev, 1, -1
[1092]	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
[685]	400	zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1220]	401	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[685]	402	pct(i)=1.-zclear(i)
[1092]	403	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
[685]	404	pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1220]	405	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[1092]	406	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
	407	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
[685]	408	pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1220]	409	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[1092]	410	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
[685]	411	pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1220]	412	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[685]	413	endif
	414	zcloud(i)=pclc(i,k)
[1092]	415	ENDDO
	416	ENDDO
	417	ELSE IF (NOVLP.EQ.2) THEN
	418	DO k = klev, 1, -1
	419	DO i = 1, klon
[685]	420	zcloud(i)=MAX(pclc(i,k),zcloud(i))
	421	pct(i)=zcloud(i)
[1092]	422	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
[685]	423	pch(i) = MIN(pclc(i,k),pch(i))
[1092]	424	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
	425	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
[685]	426	pcm(i) = MIN(pclc(i,k),pcm(i))
[1092]	427	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
[685]	428	pcl(i) = MIN(pclc(i,k),pcl(i))
	429	endif
[1092]	430	ENDDO
	431	ENDDO
	432	ELSE IF (NOVLP.EQ.3) THEN
	433	DO k = klev, 1, -1
	434	DO i = 1, klon
[685]	435	zclear(i)=zclear(i)*(1.-pclc(i,k))
	436	pct(i)=1-zclear(i)
[1092]	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))
[685]	444	endif
[1092]	445	ENDDO
[685]	446	ENDDO
[1092]	447	ENDIF
	448
	449	C
[524]	450	DO i = 1, klon
[1092]	451	c IM cf. CR pct(i)=1.-pct(i)
[524]	452	pch(i)=1.-pch(i)
	453	pcm(i)=1.-pcm(i)
	454	pcl(i)=1.-pcl(i)
	455	ENDDO
[1092]	456
[524]	457	C
	458	RETURN
	459	END

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

Download in other formats:

Original Format