Context Navigation

improvedclouds.F @ 1009

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

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

Download in other formats:

Original Format