Context Navigation

improvedclouds_mod.F @ 3001

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

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

Download in other formats:

Original Format