Context Navigation

improvedclouds_mod.F @ 2393

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

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

Download in other formats:

Original Format