Context Navigation

newmicro.F @ 5065

Last change on this file since 5065 was 1337, checked in by Laurent Fairhead, 15 years ago

Additions to aerosol outputs for CMIP5 exercise

Additions aux sorties aérosols pour l'exercice CMIP5

Michael, Anne

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

File size: 25.3 KB

Rev	Line
[1279]	1	! $Id: newmicro.F 1337 2010-04-02 11:31:05Z abarral $
	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
[1337]	13	USE phys_local_var_mod, only: scdnc,cldncl,reffclwtop,lcc,
	14	. reffclws,reffclwc,cldnvi,lcc3d,
	15	. lcc3dcon,lcc3dstra
	16	USE phys_state_var_mod, only: rnebcon,clwcon
[524]	17	IMPLICIT none
	18	c======================================================================
	19	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
	20	c Objet: Calculer epaisseur optique et emmissivite des nuages
	21	c======================================================================
	22	c Arguments:
	23	c t-------input-R-temperature
	24	c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
	25	c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
	26	c
	27	c ok_aie--input-L-apply aerosol indirect effect or not
[1279]	28	c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
	29	c mass_solu_aero_pi--input-R-dito, pre-industrial value
[524]	30	c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
	31	c bl95_b1-input-R-a parameter, may be varied for tests ( -"- )
	32	c
	33	c cldtaupi-output-R-pre-industrial value of cloud optical thickness,
	34	c needed for the diagnostics of the aerosol indirect
	35	c radiative forcing (see radlwsw)
	36	c re------output-R-Cloud droplet effective radius multiplied by fl [um]
	37	c fl------output-R-Denominator to re, introduced to avoid problems in
	38	c the averaging of the output. fl is the fraction of liquid
	39	c water clouds within a grid cell
	40	c pcltau--output-R-epaisseur optique des nuages
	41	c pclemi--output-R-emissivite des nuages (0 a 1)
	42	c======================================================================
	43	C
	44	#include "YOMCST.h"
	45	c
[766]	46	cym#include "dimensions.h"
	47	cym#include "dimphy.h"
[524]	48	#include "nuage.h"
[685]	49	cIM cf. CR: include pour NOVLP et ZEPSEC
	50	#include "radepsi.h"
	51	#include "radopt.h"
[1337]	52	c choix de l'hypothese de recouvrememnt nuageuse
	53	LOGICAL RANDOM,MAXIMUM_RANDOM,MAXIMUM,FIRST
	54	parameter (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.)
	55	c Hypoyhese de recouvrement : MAXIMUM_RANDOM
	56	INTEGER flag_max
	57	REAL phase3d(klon, klev),dh(klon, klev),pdel(klon, klev),
	58	. zrho(klon, klev)
	59	REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
	60	REAL thres_tau,thres_neb
	61	PARAMETER (thres_tau=0.3, thres_neb=0.001)
	62	REAL t_tmp
	63	REAL gravit
	64	PARAMETER (gravit=9.80616) !m/s2
	65	REAL pqlwpcon(klon, klev), pqlwpstra(klon, klev)
	66	c
[524]	67	REAL paprs(klon,klev+1), pplay(klon,klev)
	68	REAL t(klon,klev)
	69	c
	70	REAL pclc(klon,klev)
	71	REAL pqlwp(klon,klev)
	72	REAL pcltau(klon,klev), pclemi(klon,klev)
	73	c
	74	REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
	75	c
	76	LOGICAL lo
	77	c
	78	REAL cetahb, cetamb
	79	PARAMETER (cetahb = 0.45, cetamb = 0.80)
	80	C
	81	INTEGER i, k
	82	cIM: 091003 REAL zflwp, zradef, zfice, zmsac
	83	REAL zflwp(klon), zradef, zfice, zmsac
	84	cIM: 091003 rajout
	85	REAL xflwp(klon), xfiwp(klon)
	86	REAL xflwc(klon,klev), xfiwc(klon,klev)
	87	c
	88	REAL radius, rad_chaud
	89	cc PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
	90	ccc PARAMETER (rad_chaud=15.0, rad_froid=35.0)
	91	c sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0)
	92	REAL coef, coef_froi, coef_chau
	93	PARAMETER (coef_chau=0.13, coef_froi=0.09)
[1286]	94	REAL seuil_neb
	95	PARAMETER (seuil_neb=0.001)
[524]	96	INTEGER nexpo ! exponentiel pour glace/eau
	97	PARAMETER (nexpo=6)
	98	ccc PARAMETER (nexpo=1)
	99
	100	c -- sb:
	101	logical ok_newmicro
	102	c parameter (ok_newmicro=.FALSE.)
	103	cIM: 091003 real rel, tc, rei, zfiwp
	104	real rel, tc, rei, zfiwp(klon)
	105	real k_liq, k_ice0, k_ice, DF
	106	parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g
	107	parameter (DF=1.66) ! diffusivity factor
	108	c sb --
	109	cjq for the aerosol indirect effect
	110	cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
	111	cjq
	112	LOGICAL ok_aie ! Apply AIE or not?
	113	LOGICAL ok_a1lwpdep ! a1 LWP dependent?
	114
[1279]	115	REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols [ug m-3]
	116	REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value)
[524]	117	REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3]
	118	REAL re(klon, klev) ! cloud droplet effective radius [um]
	119	REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value)
	120	REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value)
	121
	122	REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
	123
	124	REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
	125
	126	REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
	127	cjq-end
[685]	128	cIM cf. CR:parametres supplementaires
	129	REAL zclear(klon)
	130	REAL zcloud(klon)
[1146]	131
	132	c **************************
	133	c * *
	134	c * DEBUT PARTIE OPTIMISEE *
	135	c * *
	136	c **************************
	137
	138	REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev)
	139	REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var
	140
[1279]	141	! Abderrahmane oct 2009
	142	Real reliq(klon, klev), reice(klon, klev)
	143
[524]	144	c
	145	c Calculer l'epaisseur optique et l'emmissivite des nuages
	146	c
[1146]	147	c IM inversion des DO
	148	xflwp = 0.d0
	149	xfiwp = 0.d0
	150	xflwc = 0.d0
	151	xfiwc = 0.d0
	152
[1334]	153	! Initialisation
	154	reliq=0.
	155	reice=0.
	156
[524]	157	DO k = 1, klev
[1146]	158	DO i = 1, klon
	159	diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
	160	ENDDO
	161	ENDDO
[524]	162
[1146]	163	IF (ok_newmicro) THEN
[524]	164
	165
[1146]	166	DO k = 1, klev
	167	DO i = 1, klon
[1286]	168	c zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
	169	zfice2(i,k) = 1.0 - (t(i,k)-t_glace_min) /
	170	& (t_glace_max-t_glace_min)
[1146]	171	zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
	172	c IM Total Liquid/Ice water content
	173	xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
	174	xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
	175	c IM In-Cloud Liquid/Ice water content
	176	c xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
	177	c xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
	178	ENDDO
	179	ENDDO
[524]	180
[1146]	181	IF (ok_aie) THEN
	182	DO k = 1, klev
	183	DO i = 1, klon
	184	! Formula "D" of Boucher and Lohmann, Tellus, 1995
	185	!
	186	cdnc(i,k) = 10.*(bl95_b0+bl95_b1
[1279]	187	& log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
[1146]	188	! Cloud droplet number concentration (CDNC) is restricted
	189	! to be within [20, 1000 cm^3]
	190	!
	191	cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
	192	!
	193	!
	194	cdnc_pi(i,k) = 10.*(bl95_b0+bl95_b1
[1305]	195	& log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))
	196	& *1.e6 !-m-3
[1146]	197	cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
	198	ENDDO
	199	ENDDO
	200	DO k = 1, klev
	201	DO i = 1, klon
	202	! rad_chaud_tab(i,k) =
	203	! & MAX(1.1e6
	204	! & ((pqlwp(i,k)pplay(i,k)/(RD * T(i,k)))
	205	! & /(4./3RPI1000.cdnc(i,k)) )*(1./3.),5.)
	206	rad_chaud_tab(i,k) =
	207	& 1.1
	208	& ((pqlwp(i,k)pplay(i,k)/(RD * T(i,k)))
	209	& /(4./3RPI1000.cdnc(i,k)) )*(1./3.)
	210	rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.)
	211	ENDDO
	212	ENDDO
	213	ELSE
	214	DO k = 1, MIN(3,klev)
	215	DO i = 1, klon
	216	rad_chaud_tab(i,k) = rad_chau2
	217	ENDDO
	218	ENDDO
	219	DO k = MIN(3,klev)+1, klev
	220	DO i = 1, klon
	221	rad_chaud_tab(i,k) = rad_chau1
	222	ENDDO
	223	ENDDO
[524]	224
[1146]	225	ENDIF
	226
	227	DO k = 1, klev
	228	! IF(.not.ok_aie) THEN
	229	rad_chaud = rad_chau1
	230	IF (k.LE.3) rad_chaud = rad_chau2
	231	! ENDIF
	232	DO i = 1, klon
	233	IF (pclc(i,k) .LE. seuil_neb) THEN
	234
	235	c -- effective cloud droplet radius (microns):
	236
	237	c for liquid water clouds:
	238	! For output diagnostics
	239	!
	240	! Cloud droplet effective radius [um]
	241	!
	242	! we multiply here with f * xl (fraction of liquid water
	243	! clouds in the grid cell) to avoid problems in the
	244	! averaging of the output.
	245	! In the output of IOIPSL, derive the real cloud droplet
	246	! effective radius as re/fl
	247	!
	248
	249	fl(i,k) = seuil_neb*(1.-zfice2(i,k))
	250	re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
	251
[1279]	252	rel = 0.
	253	rei = 0.
[1146]	254	pclc(i,k) = 0.0
	255	pcltau(i,k) = 0.0
	256	pclemi(i,k) = 0.0
	257	cldtaupi(i,k) = 0.0
	258	ELSE
[524]	259
[1146]	260	c -- liquid/ice cloud water paths:
	261
	262	zflwp_var= 1000.(1.-zfice2(i,k))pqlwp(i,k)/pclc(i,k)
	263	& *diff_paprs(i,k)
	264	zfiwp_var= 1000.zfice2(i,k)pqlwp(i,k)/pclc(i,k)
	265	& *diff_paprs(i,k)
	266
	267	c -- effective cloud droplet radius (microns):
	268
	269	c for liquid water clouds:
	270
	271	IF (ok_aie) THEN
	272	radius =
	273	& 1.1
	274	& ((pqlwp(i,k)pplay(i,k)/(RD * T(i,k)))
	275	& /(4./3.RPI1000.cdnc_pi(i,k)))*(1./3.)
	276	radius = MAX(radius*1e6, 5.)
	277
	278	tc = t(i,k)-273.15
	279	rei = 0.71*tc + 61.29
	280	if (tc.le.-81.4) rei = 3.5
	281	if (zflwp_var.eq.0.) radius = 1.
	282	if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
	283	cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
	284	& + zfiwp_var * (3.448e-03 + 2.431/rei)
[1279]	285
[1146]	286	ENDIF ! ok_aie
	287	! For output diagnostics
	288	!
	289	! Cloud droplet effective radius [um]
	290	!
	291	! we multiply here with f * xl (fraction of liquid water
	292	! clouds in the grid cell) to avoid problems in the
	293	! averaging of the output.
	294	! In the output of IOIPSL, derive the real cloud droplet
	295	! effective radius as re/fl
	296	!
	297
	298	fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))
	299	re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
	300
	301	rel = rad_chaud_tab(i,k)
	302	c for ice clouds: as a function of the ambiant temperature
	303	c [formula used by Iacobellis and Somerville (2000), with an
	304	c asymptotical value of 3.5 microns at T<-81.4 C added to be
	305	c consistent with observations of Heymsfield et al. 1986]:
	306	tc = t(i,k)-273.15
	307	rei = 0.71*tc + 61.29
	308	if (tc.le.-81.4) rei = 3.5
	309	c -- cloud optical thickness :
	310
	311	c [for liquid clouds, traditional formula,
	312	c for ice clouds, Ebert & Curry (1992)]
	313
[1279]	314	if (zflwp_var.eq.0.) rel = 1.
	315	if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1.
	316	pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
[1146]	317	& + zfiwp_var * (3.448e-03 + 2.431/rei)
	318	c -- cloud infrared emissivity:
	319
	320	c [the broadband infrared absorption coefficient is parameterized
	321	c as a function of the effective cld droplet radius]
	322
	323	c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
	324	k_ice = k_ice0 + 1.0/rei
	325
	326	pclemi(i,k) = 1.0
	327	& - EXP( -coef_chauzflwp_var - DFk_ice*zfiwp_var)
[524]	328
[1146]	329	ENDIF
[1279]	330	reliq(i,k)=rel
	331	reice(i,k)=rei
	332	! if (i.eq.1) then
	333	! print*,'Dans newmicro rel, rei :',rel, rei
	334	! print*,'Dans newmicro reliq, reice :',
	335	! $ reliq(i,k),reice(i,k)
	336	! endif
	337
[1146]	338	ENDDO
	339	ENDDO
[524]	340
[1146]	341	DO k = 1, klev
	342	DO i = 1, klon
	343	xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
	344	xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
	345	ENDDO
	346	ENDDO
[524]	347
[1146]	348	ELSE
	349	DO k = 1, klev
	350	rad_chaud = rad_chau1
	351	IF (k.LE.3) rad_chaud = rad_chau2
	352	DO i = 1, klon
	353
	354	IF (pclc(i,k) .LE. seuil_neb) THEN
[524]	355
[1146]	356	pclc(i,k) = 0.0
	357	pcltau(i,k) = 0.0
	358	pclemi(i,k) = 0.0
	359	cldtaupi(i,k) = 0.0
[524]	360
[1146]	361	ELSE
[524]	362
[1146]	363	zflwp_var = 1000.pqlwp(i,k)diff_paprs(i,k)
	364	& /pclc(i,k)
	365
	366	zfice1 = MIN(
[1286]	367	& MAX( 1.0 - (t(i,k)-t_glace_min) /
	368	& (t_glace_max-t_glace_min),0.0),1.0)**nexpo
[1146]	369
	370	radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
	371	coef = coef_chau * (1.-zfice1) + coef_froi * zfice1
[524]	372
[1146]	373	pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
	374	pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
[524]	375
[1146]	376	ENDIF
	377
	378	ENDDO
	379	ENDDO
	380	ENDIF
	381
	382	IF (.NOT.ok_aie) THEN
	383	DO k = 1, klev
	384	DO i = 1, klon
	385	cldtaupi(i,k)=pcltau(i,k)
	386	ENDDO
	387	ENDDO
	388	ENDIF
[524]	389
[1146]	390	ccc DO k = 1, klev
	391	ccc DO i = 1, klon
	392	ccc t(i,k) = t(i,k)
	393	ccc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
	394	ccc lo = pclc(i,k) .GT. (2.*1.e-5)
	395	ccc zflwp = pqlwp(i,k)1000.(paprs(i,k)-paprs(i,k+1))
	396	ccc . /(rg*pclc(i,k))
	397	ccc zradef = 10.0 + (1.-sigs(k))*45.0
	398	ccc pcltau(i,k) = 1.5 * zflwp / zradef
	399	ccc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
	400	ccc zmsac = 0.13(1.0-zfice) + 0.08zfice
	401	ccc pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
	402	ccc if (.NOT.lo) pclc(i,k) = 0.0
	403	ccc if (.NOT.lo) pcltau(i,k) = 0.0
	404	ccc if (.NOT.lo) pclemi(i,k) = 0.0
	405	ccc ENDDO
	406	ccc ENDDO
	407	ccccc print*, 'pas de nuage dans le rayonnement'
	408	ccccc DO k = 1, klev
	409	ccccc DO i = 1, klon
	410	ccccc pclc(i,k) = 0.0
	411	ccccc pcltau(i,k) = 0.0
	412	ccccc pclemi(i,k) = 0.0
	413	ccccc ENDDO
	414	ccccc ENDDO
	415	C
	416	C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
	417	C
	418	c IM cf. CR:test: calcul prenant ou non en compte le recouvrement
	419	c initialisations
[685]	420	DO i=1,klon
	421	zclear(i)=1.
	422	zcloud(i)=0.
[524]	423	pch(i)=1.0
	424	pcm(i) = 1.0
	425	pcl(i) = 1.0
	426	pctlwp(i) = 0.0
	427	ENDDO
	428	C
[685]	429	cIM cf CR DO k=1,klev
[524]	430	DO k = klev, 1, -1
[1146]	431	DO i = 1, klon
	432	pctlwp(i) = pctlwp(i)
	433	& + pqlwp(i,k)*diff_paprs(i,k)
	434	ENDDO
	435	ENDDO
	436	c IM cf. CR
	437	IF (NOVLP.EQ.1) THEN
	438	DO k = klev, 1, -1
	439	DO i = 1, klon
[685]	440	zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1279]	441	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[685]	442	pct(i)=1.-zclear(i)
[1146]	443	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
[685]	444	pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1279]	445	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[1146]	446	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
	447	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
[685]	448	pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1279]	449	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[1146]	450	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
[685]	451	pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloud(i)))
[1279]	452	& /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
[685]	453	endif
	454	zcloud(i)=pclc(i,k)
[1146]	455	ENDDO
	456	ENDDO
	457	ELSE IF (NOVLP.EQ.2) THEN
	458	DO k = klev, 1, -1
	459	DO i = 1, klon
[685]	460	zcloud(i)=MAX(pclc(i,k),zcloud(i))
	461	pct(i)=zcloud(i)
[1146]	462	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
[685]	463	pch(i) = MIN(pclc(i,k),pch(i))
[1146]	464	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
	465	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
[685]	466	pcm(i) = MIN(pclc(i,k),pcm(i))
[1146]	467	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
[685]	468	pcl(i) = MIN(pclc(i,k),pcl(i))
	469	endif
[1146]	470	ENDDO
	471	ENDDO
	472	ELSE IF (NOVLP.EQ.3) THEN
	473	DO k = klev, 1, -1
	474	DO i = 1, klon
[685]	475	zclear(i)=zclear(i)*(1.-pclc(i,k))
	476	pct(i)=1-zclear(i)
[1146]	477	IF (pplay(i,k).LE.cetahb*paprs(i,1)) THEN
	478	pch(i) = pch(i)*(1.0-pclc(i,k))
	479	ELSE IF (pplay(i,k).GT.cetahb*paprs(i,1) .AND.
	480	& pplay(i,k).LE.cetamb*paprs(i,1)) THEN
	481	pcm(i) = pcm(i)*(1.0-pclc(i,k))
	482	ELSE IF (pplay(i,k).GT.cetamb*paprs(i,1)) THEN
	483	pcl(i) = pcl(i)*(1.0-pclc(i,k))
[685]	484	endif
[1146]	485	ENDDO
[685]	486	ENDDO
[1146]	487	ENDIF
	488
	489	C
[524]	490	DO i = 1, klon
[1146]	491	c IM cf. CR pct(i)=1.-pct(i)
[524]	492	pch(i)=1.-pch(i)
	493	pcm(i)=1.-pcm(i)
	494	pcl(i)=1.-pcl(i)
	495	ENDDO
[1337]	496
	497	c ========================================================
	498	! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
	499	c ========================================================
	500	!! change by Nicolas Yan (LSCE)
	501	!! Cloud Droplet Number Concentration (CDNC) : 3D variable
	502	!! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
	503	!! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
	504	!! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
	505	!! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
	506	IF (ok_newmicro) THEN
	507	IF (ok_aie) THEN
	508	DO k = 1, klev
	509	DO i = 1, klon
	510	phase3d(i,k)=1-zfice2(i,k)
	511	IF (pclc(i,k) .LE. seuil_neb) THEN
	512	lcc3d(i,k)=seuil_neb*phase3d(i,k)
	513	ELSE
	514	lcc3d(i,k)=pclc(i,k)*phase3d(i,k)
	515	ENDIF
	516	scdnc(i,k)=lcc3d(i,k)*cdnc(i,k) ! m-3
	517	ENDDO
	518	ENDDO
	519
	520	DO i=1,klon
	521	lcc(i)=0.
	522	reffclwtop(i)=0.
	523	cldncl(i)=0.
	524	IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i) = 1.
	525	IF(MAXIMUM) tcc(i) = 0.
	526	ENDDO
	527
	528	FIRST=.TRUE.
	529
	530	DO i=1,klon
	531	DO k=klev-1,1,-1 !From TOA down
	532
	533
	534	! Test, if the cloud optical depth exceeds the necessary
	535	! threshold:
	536
	537	IF (pcltau(i,k).GT.thres_tau .AND. pclc(i,k).GT.thres_neb)
	538	. THEN
	539	! To calculate the right Temperature at cloud top,
	540	! interpolate it between layers:
	541	t_tmp = t(i,k) +
	542	. (paprs(i,k+1)-pplay(i,k))/(pplay(i,k+1)-pplay(i,k))
	543	. * ( t(i,k+1) - t(i,k) )
	544
	545	IF(MAXIMUM) THEN
	546	IF(FIRST) THEN
	547	write(,)'Hypothese de recouvrement: MAXIMUM'
	548	FIRST=.FALSE.
	549	ENDIF
	550	flag_max= -1.
	551	ftmp(i) = MAX(tcc(i),pclc(i,k))
	552	ENDIF
	553
	554	IF(RANDOM) THEN
	555	IF(FIRST) THEN
	556	write(,)'Hypothese de recouvrement: RANDOM'
	557	FIRST=.FALSE.
	558	ENDIF
	559	flag_max= 1.
	560	ftmp(i) = tcc(i) * (1-pclc(i,k))
	561	ENDIF
	562
	563	IF(MAXIMUM_RANDOM) THEN
	564	IF(FIRST) THEN
	565	write(,)'Hypothese de recouvrement: MAXIMUM_
	566	. RANDOM'
	567	FIRST=.FALSE.
	568	ENDIF
	569	flag_max= 1.
	570	ftmp(i) = tcc(i) *
	571	. (1. - MAX(pclc(i,k),pclc(i,k+1))) /
	572	. (1. - MIN(pclc(i,k+1),1.-thres_neb))
	573	ENDIF
	574	c Effective radius of cloud droplet at top of cloud (m)
	575	reffclwtop(i) = reffclwtop(i) + rad_chaud_tab(i,k) *
	576	. 1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max
	577	c CDNC at top of cloud (m-3)
	578	cldncl(i) = cldncl(i) + cdnc(i,k) * phase3d(i,k) *
	579	. (tcc(i) - ftmp(i))*flag_max
	580	c Liquid Cloud Content at top of cloud
	581	lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))*
	582	. flag_max
	583	c Total Cloud Content at top of cloud
	584	tcc(i)=ftmp(i)
	585
	586	ENDIF ! is there a visible, not-too-small cloud?
	587	ENDDO ! loop over k
	588
	589	IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
	590	ENDDO ! loop over i
	591
	592	!! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC REFFCLWS)
	593	DO i = 1, klon
	594	DO k = 1, klev
	595	pqlwpcon(i,k)=rnebcon(i,k)*clwcon(i,k) ! fraction eau liquide convective
	596	pqlwpstra(i,k)=pclc(i,k)*phase3d(i,k)-pqlwpcon(i,k) ! fraction eau liquide stratiforme
	597	IF (pqlwpstra(i,k) .LE. 0.0) pqlwpstra(i,k)=0.0
	598	! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
	599	reffclwc(i,k)=1.1
	600	& ((pqlwpcon(i,k)pplay(i,k)/(RD * T(i,k)))
	601	& /(4./3RPI1000.cdnc(i,k)) )*(1./3.)
	602	reffclwc(i,k) = MAX(reffclwc(i,k) * 1e6, 5.)
	603
	604	! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
	605	IF ((pclc(i,k)-rnebcon(i,k)) .LE. seuil_neb) THEN ! tout sous la forme convective
	606	reffclws(i,k)=0.0
	607	lcc3dstra(i,k)= 0.0
	608	ELSE
	609	reffclws(i,k) = (pclc(i,k)phase3d(i,k)
	610	& rad_chaud_tab(i,k)-
	611	& pqlwpcon(i,k)*reffclwc(i,k))
	612	IF(reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
	613	lcc3dstra(i,k)=pqlwpstra(i,k)
	614	ENDIF
	615	!Convertion from um to m
	616	IF(rnebcon(i,k). LE. seuil_neb) THEN
	617	reffclwc(i,k) = reffclwc(i,k)seuil_nebclwcon(i,k)
	618	& *1.0E-06
	619	lcc3dcon(i,k)= seuil_neb*clwcon(i,k)
	620	ELSE
	621	reffclwc(i,k) = reffclwc(i,k)*pqlwpcon(i,k)
	622	& *1.0E-06
	623	lcc3dcon(i,k) = pqlwpcon(i,k)
	624	ENDIF
	625
	626	reffclws(i,k) = reffclws(i,k)*1.0E-06
	627
	628	ENDDO !klev
	629	ENDDO !klon
	630
	631	!! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
	632	DO k = 1, klev
	633	DO i = 1, klon
	634	pdel(i,k) = paprs(i,k)-paprs(i,k+1)
	635	zrho(i,k)=pplay(i,k)/t(i,k)/RD ! kg/m3
	636	dh(i,k)=pdel(i,k)/(gravit*zrho(i,k)) ! hauteur de chaque boite (m)
	637	ENDDO
	638	ENDDO
	639	c
	640	DO i = 1, klon
	641	cldnvi(i)=0.
	642	lcc_integrat(i)=0.
	643	height(i)=0.
	644	DO k = 1, klev
	645	cldnvi(i)=cldnvi(i)+cdnc(i,k)lcc3d(i,k)dh(i,k)
	646	lcc_integrat(i)=lcc_integrat(i)+lcc3d(i,k)*dh(i,k)
	647	height(i)=height(i)+dh(i,k)
	648	ENDDO ! klev
	649	lcc_integrat(i)=lcc_integrat(i)/height(i)
	650	IF (lcc_integrat(i) .LE. 1.0E-03) THEN
	651	cldnvi(i)=cldnvi(i)*lcc(i)/seuil_neb
	652	ELSE
	653	cldnvi(i)=cldnvi(i)*lcc(i)/lcc_integrat(i)
	654	ENDIF
	655	ENDDO ! klon
	656
	657	DO i = 1, klon
	658	DO k = 1, klev
	659	IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0
	660	IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
	661	IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0
	662	IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0
	663	IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0
	664	IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
	665	ENDDO
	666	IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0
	667	IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0
	668	IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0
	669	IF (lcc(i) .LE. 0.0) lcc(i)=0.0
	670	ENDDO
	671
	672	ENDIF !ok_aie
	673	ENDIF !ok newmicro
	674	c
[524]	675	C
	676	RETURN
	677	END

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

Download in other formats:

Original Format