Context Navigation

improvedclouds.F @ 1635

Last change on this file since 1635 was 1307, checked in by tnavarro, 10 years ago
Changed two microphysics thresholds for timesteps used in Mesoscale and LES
File size: 20.9 KB

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

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

Download in other formats:

Original Format