Context Navigation

improvedclouds_mod.F @ 2417

Last change on this file since 2417 was 2407, checked in by mvals, 5 years ago

Mars GCM:
Implementation of HDO in the microphysics of water ice clouds:

improvedclouds_mod.F: calculation of the HDO flux
growthrate.F, microphys.h: addings to take into account the effect of diffusion kinetics on fractionation
callsedim_mod.F: sedimentation of HDO as an isotope of water in the microphysics case

File size: 24.2 KB

Rev	Line
[1963]	1	MODULE improvedclouds_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[1909]	7	subroutine improvedclouds(ngrid,nlay,microtimestep,
	8	& pplay,pteff,sum_subpdt,
	9	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
[633]	10	& nq,tauscaling)
[2304]	11	USE updaterad, ONLY: updaterice_micro, updaterccn
[1996]	12	USE watersat_mod, ONLY: watersat
[1036]	13	use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
	14	& igcm_h2o_ice, igcm_dust_mass,
	15	& igcm_dust_number, igcm_ccn_mass,
[2407]	16	& igcm_ccn_number,
	17	& igcm_hdo_vap,igcm_hdo_ice,
	18	& qperemin
[1047]	19	use conc_mod, only: mmean
[2304]	20	use comcstfi_h, only: pi, cpp
[358]	21	implicit none
[633]	22
	23
[358]	24	c------------------------------------------------------------------
	25	c This routine is used to form clouds when a parcel of the GCM is
	26	c saturated. It includes the ability to have supersaturation, a
	27	c computation of the nucleation rates, growthrates and the
	28	c scavenging of dust particles by clouds.
	29	c It is worth noting that the amount of dust is computed using the
	30	c dust optical depth computed in aeropacity.F. That's why
	31	c the variable called "tauscaling" is used to convert
	32	c pq(dust_mass) and pq(dust_number), which are relative
	33	c quantities, to absolute and realistic quantities stored in zq.
	34	c This has to be done to convert the inputs into absolute
	35	c values, but also to convert the outputs back into relative
	36	c values which are then used by the sedimentation and advection
	37	c schemes.
	38
	39	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
	40	c (October 2011)
[626]	41	c T. Navarro, debug,correction, new scheme (October-April 2011)
[530]	42	c A. Spiga, optimization (February 2012)
[358]	43	c------------------------------------------------------------------
	44	#include "callkeys.h"
	45	#include "microphys.h"
	46	c------------------------------------------------------------------
[1976]	47	c Inputs/outputs:
[358]	48
[1976]	49	INTEGER, INTENT(IN) :: ngrid,nlay
	50	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
	51	REAL, INTENT(IN) :: microtimestep ! pas de temps physique (s)
	52	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
	53	REAL, INTENT(IN) :: pteff(ngrid,nlay) ! temperature at the middle of the
	54	! layers (K)
	55	REAL, INTENT(IN) :: sum_subpdt(ngrid,nlay)! tendance temperature des autres
	56	! param.
	57	REAL, INTENT(IN) :: pqeff(ngrid,nlay,nq) ! traceur (kg/kg)
	58	REAL, INTENT(IN) :: sum_subpdq(ngrid,nlay,nq) ! tendance avant condensation
	59	! (kg/kg.s-1)
	60	REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
	61
	62	REAL, INTENT(OUT) :: subpdqcloud(ngrid,nlay,nq) ! tendance de la condensation
	63	! H2O(kg/kg.s-1)
	64	REAL, INTENT(OUT) :: subpdtcloud(ngrid,nlay) ! tendance temperature due
	65	! a la chaleur latente
[358]	66
	67	c------------------------------------------------------------------
	68	c Local variables:
	69
	70	LOGICAL firstcall
	71	DATA firstcall/.true./
	72	SAVE firstcall
	73
	74	REAL*8 derf ! Error function
	75	!external derf
[740]	76
[358]	77	INTEGER ig,l,i
[520]	78
[1047]	79	REAL zq(ngrid,nlay,nq) ! local value of tracers
	80	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
	81	REAL zt(ngrid,nlay) ! local value of temperature
	82	REAL zqsat(ngrid,nlay) ! saturation
[358]	83	REAL lw !Latent heat of sublimation (J.kg-1)
[633]	84	REAL cste
	85	REAL dMice ! mass of condensed ice
[2407]	86	REAL dMice_hdo ! mass of condensed HDO ice
	87	REAL alpha(ngrid,nlay) ! HDO equilibrium fractionation coefficient (Saturation=1)
	88	REAL alpha_c(ngrid,nlay) ! HDO real fractionation coefficient
[633]	89	! REAL sumcheck
[358]	90	REAL*8 ph2o ! Water vapor partial pressure (Pa)
	91	REAL*8 satu ! Water vapor saturation ratio over ice
	92	REAL*8 Mo,No
[633]	93	REAL*8 Rn, Rm, dev2, n_derf, m_derf
[358]	94	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
	95	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
[633]	96
[358]	97	REAL*8 sig ! Water-ice/air surface tension (N.m)
	98	EXTERNAL sig
	99
[633]	100	REAL dN,dM
	101	REAL rate(nbin_cld) ! nucleation rate
	102	REAL seq
	103
	104	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
	105	! (r_c in montmessin_2004)
[1047]	106	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
	107	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
[633]	108
	109	REAL res ! Resistance growth
[2407]	110	REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
[740]	111
[633]	112
[358]	113	c Parameters of the size discretization
	114	c used by the microphysical scheme
	115	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
	116	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
	117	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
	118	! Minimum boundary radius (m)
	119	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
	120	DOUBLE PRECISION vrat_cld ! Volume ratio
	121	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
	122	SAVE rb_cld
[520]	123	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
	124	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
[358]	125
[633]	126
[358]	127	REAL sigma_ice ! Variance of the ice and CCN distributions
	128	SAVE sigma_ice
[633]	129
	130
[520]	131
[420]	132	c----------------------------------
[633]	133	c TESTS
	134
	135	INTEGER countcells
[420]	136
[626]	137	LOGICAL test_flag ! flag for test/debuging outputs
[740]	138	SAVE test_flag
	139
	140
	141	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
	142	REAL res_out(ngrid,nlay)
[633]	143
[358]	144
	145	c------------------------------------------------------------------
	146
[1779]	147	! AS: firstcall OK absolute
[358]	148	IF (firstcall) THEN
[626]	149	!=============================================================
	150	! 0. Definition of the size grid
	151	!=============================================================
[358]	152	c rad_cld is the primary radius grid used for microphysics computation.
	153	c The grid spacing is computed assuming a constant volume ratio
	154	c between two consecutive bins; i.e. vrat_cld.
	155	c vrat_cld is determined from the boundary values of the size grid:
	156	c rmin_cld and rmax_cld.
	157	c The rb_cld array contains the boundary values of each rad_cld bin.
	158	c dr_cld is the width of each rad_cld bin.
	159
	160	c Volume ratio between two adjacent bins
[1307]	161	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
	162	! vrat_cld = exp(vrat_cld)
[2151]	163	vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
	164	vrat_cld = exp(vrat_cld)
[358]	165	write(,) "vrat_cld", vrat_cld
	166
	167	rb_cld(1) = rbmin_cld
	168	rad_cld(1) = rmin_cld
[530]	169	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
[1307]	170	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
[358]	171
	172	do i=1,nbin_cld-1
[531]	173	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
[358]	174	vol_cld(i+1) = vol_cld(i) * vrat_cld
	175	enddo
	176
	177	do i=1,nbin_cld
[531]	178	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
[358]	179	& rad_cld(i)
	180	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
	181	enddo
	182	rb_cld(nbin_cld+1) = rbmax_cld
	183	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
	184
	185	print*, ' '
	186	print*,'Microphysics: size bin information:'
	187	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
	188	print*,'-----------------------------------'
	189	do i=1,nbin_cld
	190	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
	191	& dr_cld(i)
	192	enddo
	193	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
	194	print*,'-----------------------------------'
	195
[541]	196	do i=1,nbin_cld+1
[740]	197	! rb_cld(i) = log(rb_cld(i))
[2151]	198	rb_cld(i) = log(rb_cld(i)) !! we save that so that it is not computed
[541]	199	!! at each timestep and gridpoint
	200	enddo
	201
[358]	202	c Contact parameter of water ice on dust ( m=cos(theta) )
	203	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	204	! mteta = 0.95
	205	write(,) 'water_param contact parameter:', mteta
	206
	207	c Volume of a water molecule (m3)
	208	vo1 = mh2o / dble(rho_ice)
	209	c Variance of the ice and CCN distributions
	210	sigma_ice = sqrt(log(1.+nuice_sed))
[633]	211
	212
[358]	213	write(,) 'Variance of ice & CCN distribs :', sigma_ice
[455]	214	write(,) 'nuice for sedimentation:', nuice_sed
[358]	215	write(,) 'Volume of a water molecule:', vo1
	216
[633]	217
	218	test_flag = .false.
	219
[358]	220	firstcall=.false.
	221	END IF
	222
[633]	223
[626]	224	!=============================================================
	225	! 1. Initialisation
	226	!=============================================================
[1909]	227	cste = 4pirho_ice*microtimestep
[626]	228
[740]	229	res_out(:,:) = 0
	230	rice(:,:) = 1.e-8
	231
[411]	232	c Initialize the tendencies
[1909]	233	subpdqcloud(1:ngrid,1:nlay,1:nq)=0
	234	subpdtcloud(1:ngrid,1:nlay)=0
[633]	235
[411]	236
[633]	237	zt(1:ngrid,1:nlay) =
[1909]	238	& pteff(1:ngrid,1:nlay) +
	239	& sum_subpdt(1:ngrid,1:nlay) * microtimestep
[358]	240
[633]	241	zq(1:ngrid,1:nlay,1:nq) =
[1909]	242	& pqeff(1:ngrid,1:nlay,1:nq) +
	243	& sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
[633]	244
	245
	246	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
	247	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
[768]	248
	249	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
[633]	250
[626]	251	!=============================================================
	252	! 2. Compute saturation
	253	!=============================================================
[358]	254
[633]	255
[541]	256	dev2 = 1. / ( sqrt(2.) * sigma_ice )
[358]	257
[1047]	258	call watersat(ngrid*nlay,zt,pplay,zqsat)
[358]	259
[633]	260	countcells = 0
[358]	261
	262	c Main loop over the GCM's grid
	263	DO l=1,nlay
	264	DO ig=1,ngrid
	265
	266	c Get the partial pressure of water vapor and its saturation ratio
[663]	267	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
[689]	268	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
[358]	269
[626]	270	!=============================================================
	271	! 3. Nucleation
	272	!=============================================================
	273
[633]	274	IF ( satu .ge. 1. ) THEN ! if there is condensation
	275
[740]	276	call updaterccn(zq(ig,l,igcm_dust_mass),
	277	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
[633]	278
	279
[358]	280	c Expand the dust moments into a binned distribution
[633]	281	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
[626]	282	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
[358]	283	Rn = rdust(ig,l)
[2151]	284	Rn = -log(Rn)
[530]	285	Rm = Rn - 3. * sigma_ice*sigma_ice
[626]	286	n_derf = derf( (rb_cld(1)+Rn) *dev2)
	287	m_derf = derf( (rb_cld(1)+Rm) *dev2)
[358]	288	do i = 1, nbin_cld
[626]	289	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
	290	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
	291	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
	292	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
	293	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
	294	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
[358]	295	enddo
[530]	296
[358]	297	! sumcheck = 0
	298	! do i = 1, nbin_cld
	299	! sumcheck = sumcheck + n_aer(i)
	300	! enddo
	301	! sumcheck = abs(sumcheck/No - 1)
	302	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	303	! print*, "WARNING, No sumcheck PROBLEM"
	304	! print*, "sumcheck, No",sumcheck, No
	305	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	306	! print*, "Dust binned distribution", n_aer
	307	! endif
	308	!
	309	! sumcheck = 0
	310	! do i = 1, nbin_cld
[411]	311	! sumcheck = sumcheck + m_aer(i)
[358]	312	! enddo
	313	! sumcheck = abs(sumcheck/Mo - 1)
	314	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
	315	! print*, "WARNING, Mo sumcheck PROBLEM"
[411]	316	! print*, "sumcheck, Mo",sumcheck, Mo
[358]	317	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
	318	! print*, "Dust binned distribution", m_aer
	319	! endif
[633]	320
	321
[358]	322	c Get the rates of nucleation
	323	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
[411]	324
[358]	325	dN = 0.
	326	dM = 0.
	327	do i = 1, nbin_cld
[1909]	328	n_aer(i) = n_aer(i)/( 1. + rate(i)*microtimestep)
	329	m_aer(i) = m_aer(i)/( 1. + rate(i)*microtimestep)
	330	dN = dN + n_aer(i) * rate(i) * microtimestep
	331	dM = dM + m_aer(i) * rate(i) * microtimestep
[358]	332	enddo
	333
[633]	334
	335	c Update Dust particles
[626]	336	zq(ig,l,igcm_dust_mass) =
[740]	337	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	338	zq(ig,l,igcm_dust_number) =
[740]	339	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[633]	340	c Update CCNs
[626]	341	zq(ig,l,igcm_ccn_mass) =
[740]	342	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	343	zq(ig,l,igcm_ccn_number) =
[740]	344	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	345
[633]	346	ENDIF ! of is satu >1
[626]	347
	348	!=============================================================
	349	! 4. Ice growth: scheme for radius evolution
	350	!=============================================================
	351
[633]	352	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
	353	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
	354	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
	355
	356	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
	357
[740]	358
	359	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
	360	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
	361	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
[633]	362
	363	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
[626]	364
	365	c saturation at equilibrium
[740]	366	c rice should not be too small, otherwise seq value is not valid
	367	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
	368	& max(rice(ig,l),1.e-7)))
	369
[633]	370	c get resistance growth
	371	call growthrate(zt(ig,l),pplay(ig,l),
[2407]	372	& real(ph2o/satu),rice(ig,l),res,Dv)
[358]	373
[740]	374	res_out(ig,l) = res
[626]	375
[633]	376	ccccccc implicit scheme of mass growth
[626]	377
[633]	378	dMice =
	379	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
	380	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
[358]	381
[626]	382
[633]	383	! With the above scheme, dMice cannot be bigger than vapor,
	384	! but can be bigger than all available ice.
	385	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
	386	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
	387
	388	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
	389	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
	390
	391
	392	countcells = countcells + 1
	393
	394	! latent heat release
	395	lw=(2834.3-0.28*(zt(ig,l)-To)-
	396	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
[1909]	397	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
[358]	398
[2407]	399	c Special case of the isotope of water HDO
	400	if (hdo) then
	401	!! condensation
	402	if (dMice.gt.0.0) then
	403	!! do we use fractionation?
	404	if (hdofrac) then
	405	!! Calculation of the HDO vapor coefficient
	406	Dv_hdo = 1./3. * sqrt( 8kbzzt(ig,l)/(pi*mhdo/nav) )
	407	& * kbz * zt(ig,l) /
	408	& ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo)
	409	& * sqrt(1.+mhdo/mco2) )
	410	!! Calculation of the fractionnation coefficient at equilibrium
	411	alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
	412	c alpha = exp(13525./zt(ig,l)**2.-5.59d-2) !Lamb
	413	!! Calculation of the 'real' fractionnation coefficient
	414	alpha_c(ig,l) = (alpha(ig,l)*satu)/
	415	& ( (alpha(ig,l)(Dv/Dv_hdo)(satu-1.)) + 1.)
	416	c alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
	417	else
	418	alpha_c(ig,l) = 1.d0
	419	endif
	420	if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
	421	dMice_hdo=
	422	& dMicealpha_c(ig,l)
	423	& ( zq0(ig,l,igcm_hdo_vap)
	424	& /zq0(ig,l,igcm_h2o_vap) )
	425	else
	426	dMice_hdo=0.
	427	endif
	428	!! sublimation
	429	else
	430	if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
	431	dMice_hdo=
	432	& dMice*
	433	& ( zq0(ig,l,igcm_hdo_ice)
	434	& /zq0(ig,l,igcm_h2o_ice) )
	435	else
	436	dMice_hdo=0.
	437	endif
	438	endif !if (dMice.gt.0.0)
	439
	440	dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
	441	dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
	442
	443	zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
	444	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
	445
	446	endif ! if (hdo)
[358]	447
[626]	448	!=============================================================
	449	! 5. Dust cores released, tendancies, latent heat, etc ...
	450	!=============================================================
	451
[358]	452	c If all the ice particles sublimate, all the condensation
[626]	453	c nuclei are released:
[1307]	454	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
[633]	455
[626]	456	c Water
	457	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
	458	& + zq(ig,l,igcm_h2o_ice)
[358]	459	zq(ig,l,igcm_h2o_ice) = 0.
[2407]	460	if (hdo) then
	461	zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)
	462	& + zq(ig,l,igcm_hdo_ice)
	463	zq(ig,l,igcm_hdo_ice) = 0.
	464	endif
[358]	465	c Dust particles
[626]	466	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
	467	& + zq(ig,l,igcm_ccn_mass)
	468	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
	469	& + zq(ig,l,igcm_ccn_number)
[358]	470	c CCNs
	471	zq(ig,l,igcm_ccn_mass) = 0.
	472	zq(ig,l,igcm_ccn_number) = 0.
[633]	473
[358]	474	endif
[411]	475
[633]	476	ENDIF !of if Nccn>1
	477
[626]	478	ENDDO ! of ig loop
	479	ENDDO ! of nlayer loop
[520]	480
[633]	481
	482	! Get cloud tendencies
[1909]	483	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
[740]	484	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
[1909]	485	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
	486	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
[740]	487	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
[1909]	488	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
[2407]	489	if (hdo) then
	490	subpdqcloud(1:ngrid,1:nlay,igcm_hdo_vap) =
	491	& (zq(1:ngrid,1:nlay,igcm_hdo_vap) -
	492	& zq0(1:ngrid,1:nlay,igcm_hdo_vap))/microtimestep
	493	subpdqcloud(1:ngrid,1:nlay,igcm_hdo_ice) =
	494	& (zq(1:ngrid,1:nlay,igcm_hdo_ice) -
	495	& zq0(1:ngrid,1:nlay,igcm_hdo_ice))/microtimestep
	496	endif
[1909]	497	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
[740]	498	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
[1909]	499	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
	500	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
[740]	501	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
[1909]	502	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
[411]	503
[740]	504	if (scavenging) then
[633]	505
[1909]	506	subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
[740]	507	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
[1909]	508	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
	509	subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
[740]	510	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
[1909]	511	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
[740]	512
	513	endif
[633]	514
[740]	515
	516
[626]	517	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	518	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
[740]	519	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
[626]	520	IF (test_flag) then
	521
[633]	522	! error2d(:) = 0.
[740]	523	DO l=1,nlay
	524	DO ig=1,ngrid
[633]	525	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
[740]	526	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
	527	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
	528	ENDDO
	529	ENDDO
[420]	530
[633]	531	print*, 'count is ',countcells, ' i.e. ',
	532	& countcells100/(nlayngrid), '% for microphys computation'
[358]	533
[1212]	534	#ifndef MESOSCALE
[633]	535	! IF (ngrid.ne.1) THEN ! 3D
[626]	536	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
	537	! & satu_out)
	538	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
	539	! & dM_out)
	540	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
	541	! & dN_out)
[633]	542	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
	543	! & error2d)
[626]	544	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
	545	! & zqsat)
[633]	546	! ENDIF
[358]	547
[633]	548	! IF (ngrid.eq.1) THEN ! 1D
	549	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
	550	! & error_out)
[740]	551	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
	552	& res_out)
	553	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
	554	& satubf)
	555	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
	556	& satuaf)
	557	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
[1020]	558	& zq0(1,1,igcm_h2o_vap))
[740]	559	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
[1020]	560	& zq(1,1,igcm_h2o_vap))
[740]	561	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
[1020]	562	& zq0(1,1,igcm_h2o_ice))
[740]	563	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
[1020]	564	& zq(1,1,igcm_h2o_ice))
[740]	565	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
[1020]	566	& zq0(1,1,igcm_ccn_number))
[740]	567	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
[1020]	568	& zq(1,1,igcm_ccn_number))
[626]	569	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
	570	c & gr_out)
	571	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
	572	c & rate_out)
	573	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
	574	c & dM_out)
	575	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
	576	c & dN_out)
[740]	577	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
	578	& zqsat)
[633]	579	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
	580	! & satu_out)
[740]	581	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
	582	& rice)
[633]	583	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
	584	! & rdust)
	585	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
	586	! & rsedcloud)
	587	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
	588	! & rhocloud)
	589	! ENDIF
[1212]	590	#endif
[626]	591
	592	ENDIF ! endif test_flag
	593	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	594	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	595	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	596
[358]	597	return
[1963]	598
[626]	599
	600
	601
	602	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	603	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	604	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
	605	c It is an analytical solution to the ice radius growth equation,
	606	c with the approximation of a constant 'reduced' cunningham correction factor
	607	c (lambda in growthrate.F) taken at radius req instead of rice
	608	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	609	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	610
[633]	611	c subroutine phi(rice,req,coeff1,coeff2,time)
	612	c
	613	c implicit none
	614	c
	615	c ! inputs
	616	c real rice ! ice radius
	617	c real req ! ice radius at equilibirum
	618	c real coeff1 ! coeff for the log
	619	c real coeff2 ! coeff for the arctan
	620	c
	621	c ! output
	622	c real time
	623	c
	624	c !local
	625	c real var
	626	c
	627	c ! 1.73205 is sqrt(3)
	628	c
	629	c var = max(
	630	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
	631	c
	632	c time =
	633	c & coeff1 *
	634	c & log( var )
	635	c & + coeff2 * 1.73205 *
	636	c & atan( (2rice+req) / (1.73205req) )
	637	c
	638	c return
	639	c end
[626]	640
	641
	642
[1963]	643	END SUBROUTINE improvedclouds
	644
	645	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format