Context Navigation

improvedclouds_mod.F @ 2829

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

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

Download in other formats:

Original Format