Context Navigation

improvedclouds_mod.F @ 3026

Last change on this file since 3026 was 3017, checked in by jnaar, 17 months ago
Change the way zdq is estimated to fasten (by a lot) the adaptative subtimestep of water ice clouds. JN
File size: 29.3 KB

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

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

Download in other formats:

Original Format