Context Navigation

improvedclouds.F @ 758

Last change on this file since 758 was 740, checked in by tnavarro, 13 years ago
module for ice and radius radius computation
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	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
	245
	246	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
	247	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
	248
[626]	249	!=============================================================
	250	! 2. Compute saturation
	251	!=============================================================
[358]	252
[633]	253
[541]	254	dev2 = 1. / ( sqrt(2.) * sigma_ice )
[358]	255
	256	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	257
[633]	258	countcells = 0
[358]	259
	260	c Main loop over the GCM's grid
	261	DO l=1,nlay
	262	DO ig=1,ngrid
	263
	264	c Get the partial pressure of water vapor and its saturation ratio
[663]	265	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
[689]	266	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
[358]	267
[626]	268	!=============================================================
	269	! 3. Nucleation
	270	!=============================================================
	271
[633]	272	IF ( satu .ge. 1. ) THEN ! if there is condensation
	273
[740]	274	call updaterccn(zq(ig,l,igcm_dust_mass),
	275	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
[633]	276
	277
[358]	278	c Expand the dust moments into a binned distribution
[633]	279	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
[626]	280	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
[358]	281	Rn = rdust(ig,l)
[530]	282	Rn = -dlog(Rn)
	283	Rm = Rn - 3. * sigma_ice*sigma_ice
[626]	284	n_derf = derf( (rb_cld(1)+Rn) *dev2)
	285	m_derf = derf( (rb_cld(1)+Rm) *dev2)
[358]	286	do i = 1, nbin_cld
[626]	287	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
	288	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
	289	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
	290	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
	291	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
	292	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
[358]	293	enddo
[530]	294
[358]	295	! sumcheck = 0
	296	! do i = 1, nbin_cld
	297	! sumcheck = sumcheck + n_aer(i)
	298	! enddo
	299	! sumcheck = abs(sumcheck/No - 1)
	300	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	301	! print*, "WARNING, No sumcheck PROBLEM"
	302	! print*, "sumcheck, No",sumcheck, No
	303	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	304	! print*, "Dust binned distribution", n_aer
	305	! endif
	306	!
	307	! sumcheck = 0
	308	! do i = 1, nbin_cld
[411]	309	! sumcheck = sumcheck + m_aer(i)
[358]	310	! enddo
	311	! sumcheck = abs(sumcheck/Mo - 1)
	312	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
	313	! print*, "WARNING, Mo sumcheck PROBLEM"
[411]	314	! print*, "sumcheck, Mo",sumcheck, Mo
[358]	315	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
	316	! print*, "Dust binned distribution", m_aer
	317	! endif
[633]	318
	319
[358]	320	c Get the rates of nucleation
	321	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
[411]	322
[358]	323	dN = 0.
	324	dM = 0.
	325	do i = 1, nbin_cld
[633]	326	n_aer(i) = n_aer(i)/( 1. + rate(i)*ptimestep)
	327	m_aer(i) = m_aer(i)/( 1. + rate(i)*ptimestep)
[358]	328	dN = dN + n_aer(i) * rate(i) * ptimestep
	329	dM = dM + m_aer(i) * rate(i) * ptimestep
	330	enddo
	331
[633]	332
	333	c Update Dust particles
[626]	334	zq(ig,l,igcm_dust_mass) =
[740]	335	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	336	zq(ig,l,igcm_dust_number) =
[740]	337	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[633]	338	c Update CCNs
[626]	339	zq(ig,l,igcm_ccn_mass) =
[740]	340	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	341	zq(ig,l,igcm_ccn_number) =
[740]	342	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	343
[633]	344	ENDIF ! of is satu >1
[626]	345
	346	!=============================================================
	347	! 4. Ice growth: scheme for radius evolution
	348	!=============================================================
	349
[633]	350	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
	351	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
	352	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
	353
	354	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
	355
[740]	356
	357	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
	358	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
	359	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
[633]	360
	361	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
[626]	362
	363	c saturation at equilibrium
[740]	364	c rice should not be too small, otherwise seq value is not valid
	365	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
	366	& max(rice(ig,l),1.e-7)))
	367
[633]	368	c get resistance growth
	369	call growthrate(zt(ig,l),pplay(ig,l),
	370	& real(ph2o/satu),rice(ig,l),res)
[358]	371
[740]	372	res_out(ig,l) = res
[626]	373
[633]	374	ccccccc implicit scheme of mass growth
[626]	375
[633]	376	dMice =
	377	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
	378	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
[358]	379
[626]	380
[633]	381	! With the above scheme, dMice cannot be bigger than vapor,
	382	! but can be bigger than all available ice.
	383	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
	384	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
	385
	386	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
	387	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
	388
	389
	390	countcells = countcells + 1
	391
	392	! latent heat release
	393	lw=(2834.3-0.28*(zt(ig,l)-To)-
	394	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
	395	pdtcloud(ig,l)= dMice*lw/cpp/ptimestep
[358]	396
[626]	397
[358]	398
[626]	399	!=============================================================
	400	! 5. Dust cores released, tendancies, latent heat, etc ...
	401	!=============================================================
	402
[358]	403	c If all the ice particles sublimate, all the condensation
[626]	404	c nuclei are released:
[633]	405	if (zq(ig,l,igcm_h2o_ice).le.1.e-8) then
	406
[626]	407	c Water
	408	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
	409	& + zq(ig,l,igcm_h2o_ice)
[358]	410	zq(ig,l,igcm_h2o_ice) = 0.
	411	c Dust particles
[626]	412	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
	413	& + zq(ig,l,igcm_ccn_mass)
	414	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
	415	& + zq(ig,l,igcm_ccn_number)
[358]	416	c CCNs
	417	zq(ig,l,igcm_ccn_mass) = 0.
	418	zq(ig,l,igcm_ccn_number) = 0.
[633]	419
[358]	420	endif
[411]	421
[633]	422	ENDIF !of if Nccn>1
	423
[626]	424	ENDDO ! of ig loop
	425	ENDDO ! of nlayer loop
[520]	426
[633]	427
	428	! Get cloud tendencies
[740]	429	pdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
	430	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
	431	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/ptimestep
	432	pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
	433	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
	434	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
	435	pdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
	436	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
	437	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
	438	pdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
	439	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
	440	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
[411]	441
[740]	442	if (scavenging) then
[633]	443
[740]	444	pdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
	445	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
	446	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
	447	pdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
	448	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
	449	& zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
	450
	451	endif
[633]	452
[740]	453
	454
[626]	455	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	456	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
[740]	457	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
[626]	458	IF (test_flag) then
	459
[633]	460	! error2d(:) = 0.
[740]	461	DO l=1,nlay
	462	DO ig=1,ngrid
[633]	463	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
[740]	464	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
	465	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
	466	ENDDO
	467	ENDDO
[420]	468
[633]	469	print*, 'count is ',countcells, ' i.e. ',
	470	& countcells100/(nlayngrid), '% for microphys computation'
[358]	471
[633]	472	! IF (ngrid.ne.1) THEN ! 3D
[626]	473	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
	474	! & satu_out)
	475	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
	476	! & dM_out)
	477	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
	478	! & dN_out)
[633]	479	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
	480	! & error2d)
[626]	481	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
	482	! & zqsat)
[633]	483	! ENDIF
[358]	484
[633]	485	! IF (ngrid.eq.1) THEN ! 1D
	486	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
	487	! & error_out)
[740]	488	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
	489	& res_out)
	490	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
	491	& satubf)
	492	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
	493	& satuaf)
	494	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
	495	& zq0(1,:,igcm_h2o_vap))
	496	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
	497	& zq(1,:,igcm_h2o_vap))
	498	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
	499	& zq0(1,:,igcm_h2o_ice))
	500	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
	501	& zq(1,:,igcm_h2o_ice))
	502	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
	503	& zq0(1,:,igcm_ccn_number))
	504	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
	505	& zq(1,:,igcm_ccn_number))
[626]	506	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
	507	c & gr_out)
	508	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
	509	c & rate_out)
	510	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
	511	c & dM_out)
	512	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
	513	c & dN_out)
[740]	514	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
	515	& zqsat)
[633]	516	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
	517	! & satu_out)
[740]	518	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
	519	& rice)
[633]	520	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
	521	! & rdust)
	522	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
	523	! & rsedcloud)
	524	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
	525	! & rhocloud)
	526	! ENDIF
[626]	527
	528	ENDIF ! endif test_flag
	529	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	530	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	531	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	532
[358]	533	return
	534	end
[626]	535
	536
	537
	538	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	539	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	540	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
	541	c It is an analytical solution to the ice radius growth equation,
	542	c with the approximation of a constant 'reduced' cunningham correction factor
	543	c (lambda in growthrate.F) taken at radius req instead of rice
	544	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	545	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	546
[633]	547	c subroutine phi(rice,req,coeff1,coeff2,time)
	548	c
	549	c implicit none
	550	c
	551	c ! inputs
	552	c real rice ! ice radius
	553	c real req ! ice radius at equilibirum
	554	c real coeff1 ! coeff for the log
	555	c real coeff2 ! coeff for the arctan
	556	c
	557	c ! output
	558	c real time
	559	c
	560	c !local
	561	c real var
	562	c
	563	c ! 1.73205 is sqrt(3)
	564	c
	565	c var = max(
	566	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
	567	c
	568	c time =
	569	c & coeff1 *
	570	c & log( var )
	571	c & + coeff2 * 1.73205 *
	572	c & atan( (2rice+req) / (1.73205req) )
	573	c
	574	c return
	575	c end
[626]	576
	577
	578

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

Download in other formats:

Original Format