Context Navigation

improvedclouds.F @ 693

Last change on this file since 693 was 689, checked in by flefevre, 13 years ago
Minor change. Better handling of very small water mixing ratio at saturation (for very low temperatures).
File size: 20.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
[358]	7	implicit none
[633]	8
	9
[358]	10	c------------------------------------------------------------------
	11	c This routine is used to form clouds when a parcel of the GCM is
	12	c saturated. It includes the ability to have supersaturation, a
	13	c computation of the nucleation rates, growthrates and the
	14	c scavenging of dust particles by clouds.
	15	c It is worth noting that the amount of dust is computed using the
	16	c dust optical depth computed in aeropacity.F. That's why
	17	c the variable called "tauscaling" is used to convert
	18	c pq(dust_mass) and pq(dust_number), which are relative
	19	c quantities, to absolute and realistic quantities stored in zq.
	20	c This has to be done to convert the inputs into absolute
	21	c values, but also to convert the outputs back into relative
	22	c values which are then used by the sedimentation and advection
	23	c schemes.
	24
	25	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
	26	c (October 2011)
[626]	27	c T. Navarro, debug,correction, new scheme (October-April 2011)
[530]	28	c A. Spiga, optimization (February 2012)
[358]	29	c------------------------------------------------------------------
	30	#include "dimensions.h"
	31	#include "dimphys.h"
	32	#include "comcstfi.h"
	33	#include "callkeys.h"
	34	#include "tracer.h"
	35	#include "comgeomfi.h"
	36	#include "dimradmars.h"
	37	#include "microphys.h"
[663]	38	#include "conc.h"
[358]	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)
	54	REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust
	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
	71
	72	REAL CBRT
	73	EXTERNAL CBRT
	74
	75	INTEGER ig,l,i
[520]	76
[358]	77
	78	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
	79	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
	80	REAL zt(ngridmx,nlayermx) ! local value of temperature
	81	REAL zqsat(ngridmx,nlayermx) ! saturation
	82	REAL lw !Latent heat of sublimation (J.kg-1)
[633]	83	REAL cste
	84	REAL dMice ! mass of condensed ice
	85	! REAL sumcheck
[358]	86	REAL*8 ph2o ! Water vapor partial pressure (Pa)
	87	REAL*8 satu ! Water vapor saturation ratio over ice
	88	REAL*8 Mo,No
[633]	89	REAL*8 Rn, Rm, dev2, n_derf, m_derf
[358]	90	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
	91	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
[633]	92
[358]	93	REAL*8 sig ! Water-ice/air surface tension (N.m)
	94	EXTERNAL sig
	95
[633]	96	REAL dN,dM
	97	REAL rate(nbin_cld) ! nucleation rate
	98	REAL seq
	99
	100	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
	101	! (r_c in montmessin_2004)
	102	REAL rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3)
	103	REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
	104
	105	REAL res ! Resistance growth
	106
	107
[358]	108	c Parameters of the size discretization
	109	c used by the microphysical scheme
	110	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
	111	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
	112	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
	113	! Minimum boundary radius (m)
	114	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
	115	DOUBLE PRECISION vrat_cld ! Volume ratio
	116	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
	117	SAVE rb_cld
[520]	118	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
	119	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
[358]	120
[633]	121
[358]	122	REAL sigma_ice ! Variance of the ice and CCN distributions
	123	SAVE sigma_ice
[633]	124
	125
[520]	126
[420]	127	c----------------------------------
[633]	128	c TESTS
	129
	130	INTEGER countcells
[420]	131
[626]	132	LOGICAL test_flag ! flag for test/debuging outputs
[633]	133	SAVE test_flag
	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
[633]	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
[633]	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
[633]	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
[411]	219	c Initialize the tendencies
	220	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	221	pdtcloud(1:ngrid,1:nlay)=0
[633]	222
	223	c Initialize the tendencies
	224	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	225	pdtcloud(1:ngrid,1:nlay)=0
[411]	226
[633]	227	c Initialize the tendencies
	228	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	229	pdtcloud(1:ngrid,1:nlay)=0
	230
	231
	232	zt(1:ngrid,1:nlay) =
	233	& pt(1:ngrid,1:nlay) +
	234	& pdt(1:ngrid,1:nlay) * ptimestep
[358]	235
[633]	236	zq(1:ngrid,1:nlay,1:nq) =
	237	& pq(1:ngrid,1:nlay,1:nq) +
	238	& pdq(1:ngrid,1:nlay,1:nq) * ptimestep
	239
	240	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
	241
	242	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
	243	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
	244
[626]	245	!=============================================================
	246	! 2. Compute saturation
	247	!=============================================================
[358]	248
[633]	249
[541]	250	dev2 = 1. / ( sqrt(2.) * sigma_ice )
[358]	251
	252	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	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
	270	c Update rdust from last tendencies
	271	rdust(ig,l)=
	272	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
	273	& (zq(ig,l,igcm_dust_number)+1.e-30))
	274	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
	275
	276
[358]	277	c Expand the dust moments into a binned distribution
[633]	278	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) + 1.e-30
[626]	279	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
[358]	280	Rn = rdust(ig,l)
[530]	281	Rn = -dlog(Rn)
	282	Rm = Rn - 3. * sigma_ice*sigma_ice
[626]	283	n_derf = derf( (rb_cld(1)+Rn) *dev2)
	284	m_derf = derf( (rb_cld(1)+Rm) *dev2)
[358]	285	do i = 1, nbin_cld
[626]	286	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
	287	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
	288	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
	289	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
	290	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
	291	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
[358]	292	enddo
[530]	293
[358]	294	! sumcheck = 0
	295	! do i = 1, nbin_cld
	296	! sumcheck = sumcheck + n_aer(i)
	297	! enddo
	298	! sumcheck = abs(sumcheck/No - 1)
	299	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	300	! print*, "WARNING, No sumcheck PROBLEM"
	301	! print*, "sumcheck, No",sumcheck, No
	302	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	303	! print*, "Dust binned distribution", n_aer
	304	! endif
	305	!
	306	! sumcheck = 0
	307	! do i = 1, nbin_cld
[411]	308	! sumcheck = sumcheck + m_aer(i)
[358]	309	! enddo
	310	! sumcheck = abs(sumcheck/Mo - 1)
	311	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
	312	! print*, "WARNING, Mo sumcheck PROBLEM"
[411]	313	! print*, "sumcheck, Mo",sumcheck, Mo
[358]	314	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
	315	! print*, "Dust binned distribution", m_aer
	316	! endif
[633]	317
	318
[358]	319	c Get the rates of nucleation
	320	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
[411]	321
[358]	322	dN = 0.
	323	dM = 0.
	324	do i = 1, nbin_cld
[633]	325	n_aer(i) = n_aer(i)/( 1. + rate(i)*ptimestep)
	326	m_aer(i) = m_aer(i)/( 1. + rate(i)*ptimestep)
[358]	327	dN = dN + n_aer(i) * rate(i) * ptimestep
	328	dM = dM + m_aer(i) * rate(i) * ptimestep
	329	enddo
	330
[633]	331
	332	c Update Dust particles
[626]	333	zq(ig,l,igcm_dust_mass) =
	334	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig)
	335	zq(ig,l,igcm_dust_number) =
	336	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
[633]	337	c Update CCNs
[626]	338	zq(ig,l,igcm_ccn_mass) =
	339	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig)
	340	zq(ig,l,igcm_ccn_number) =
	341	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
	342
[633]	343	ENDIF ! of is satu >1
[626]	344
	345	!=============================================================
	346	! 4. Ice growth: scheme for radius evolution
	347	!=============================================================
	348
[633]	349	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
	350	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
	351	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
	352
	353	IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
	354
	355
[626]	356	Mo = zq(ig,l,igcm_h2o_ice) +
[633]	357	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
	358	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
[626]	359
[633]	360
[626]	361	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
	362	& + zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
	363	& / Mo * rho_dust
	364	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
	365
	366
	367	c ice crystal radius
	368	rice (ig,l) =
[633]	369	& CBRT( real(Mo/No) * 0.75 / pi / rhocloud(ig,l) )
[626]	370
	371
	372	c saturation at equilibrium
	373	seq = exp( 2.sig(zt(ig,l))mh2o /
[358]	374	& (rho_icergpzt(ig,l)*rice(ig,l)) )
	375
[633]	376	c get resistance growth
	377	call growthrate(zt(ig,l),pplay(ig,l),
	378	& real(ph2o/satu),rice(ig,l),res)
[358]	379
[626]	380
[633]	381	ccccccc implicit scheme of mass growth
[626]	382
[633]	383	dMice =
	384	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
	385	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
[358]	386
[626]	387
[633]	388	! With the above scheme, dMice cannot be bigger than vapor,
	389	! but can be bigger than all available ice.
	390	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
	391	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
	392
	393	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
	394	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
	395
	396
	397	countcells = countcells + 1
	398
	399	! latent heat release
	400	lw=(2834.3-0.28*(zt(ig,l)-To)-
	401	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
	402	pdtcloud(ig,l)= dMice*lw/cpp/ptimestep
[358]	403
[626]	404
[358]	405
[626]	406	!=============================================================
	407	! 5. Dust cores released, tendancies, latent heat, etc ...
	408	!=============================================================
	409
[358]	410	c If all the ice particles sublimate, all the condensation
[626]	411	c nuclei are released:
[633]	412	if (zq(ig,l,igcm_h2o_ice).le.1.e-8) then
	413
[626]	414	c Water
	415	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
	416	& + zq(ig,l,igcm_h2o_ice)
[358]	417	zq(ig,l,igcm_h2o_ice) = 0.
	418	c Dust particles
[626]	419	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
	420	& + zq(ig,l,igcm_ccn_mass)
	421	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
	422	& + zq(ig,l,igcm_ccn_number)
[358]	423	c CCNs
	424	zq(ig,l,igcm_ccn_mass) = 0.
	425	zq(ig,l,igcm_ccn_number) = 0.
[633]	426
[358]	427	endif
[411]	428
[633]	429	ENDIF !of if Nccn>1
	430
[626]	431	ENDDO ! of ig loop
	432	ENDDO ! of nlayer loop
[520]	433
[633]	434
	435	! Get cloud tendencies
	436	pdqcloud(1:ngrid,1:nlay,1:nq) =
	437	& (zq(1:ngrid,1:nlay,1:nq)-zq0(1:ngrid,1:nlay,1:nq))/ptimestep
[411]	438
[633]	439
	440
[626]	441	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	442	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	443	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	444	IF (test_flag) then
	445
[633]	446	! error2d(:) = 0.
	447	! DO l=1,nlay
	448	! DO ig=1,ngrid
	449	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
	450	! ENDDO
	451	! ENDDO
[420]	452
[633]	453	print*, 'count is ',countcells, ' i.e. ',
	454	& countcells100/(nlayngrid), '% for microphys computation'
[358]	455
[633]	456	! IF (ngrid.ne.1) THEN ! 3D
[626]	457	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
	458	! & satu_out)
	459	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
	460	! & dM_out)
	461	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
	462	! & dN_out)
[633]	463	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
	464	! & error2d)
[626]	465	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
	466	! & zqsat)
[633]	467	! ENDIF
[358]	468
[633]	469	! IF (ngrid.eq.1) THEN ! 1D
	470	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
	471	! & error_out)
	472	! call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
	473	! & satubf)
	474	! call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
	475	! & satuaf)
	476	! call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
	477	! & zq0(1,:,igcm_h2o_vap))
	478	! call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
	479	! & zq(1,:,igcm_h2o_vap))
	480	! call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
	481	! & zq0(1,:,igcm_h2o_ice))
	482	! call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
	483	! & zq(1,:,igcm_h2o_ice))
	484	! call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
	485	! & ccnbf)
	486	! call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
	487	! & ccnaf)
[626]	488	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
	489	c & gr_out)
	490	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
	491	c & rate_out)
	492	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
	493	c & dM_out)
	494	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
	495	c & dN_out)
[633]	496	! call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
	497	! & zqsat)
	498	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
	499	! & satu_out)
	500	! call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
	501	! & rice)
	502	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
	503	! & rdust)
	504	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
	505	! & rsedcloud)
	506	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
	507	! & rhocloud)
	508	! ENDIF
[626]	509
	510	ENDIF ! endif test_flag
	511	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	512	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	513	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	514
[358]	515	return
	516	end
[626]	517
	518
	519
	520	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	521	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	522	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
	523	c It is an analytical solution to the ice radius growth equation,
	524	c with the approximation of a constant 'reduced' cunningham correction factor
	525	c (lambda in growthrate.F) taken at radius req instead of rice
	526	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	527	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	528
[633]	529	c subroutine phi(rice,req,coeff1,coeff2,time)
	530	c
	531	c implicit none
	532	c
	533	c ! inputs
	534	c real rice ! ice radius
	535	c real req ! ice radius at equilibirum
	536	c real coeff1 ! coeff for the log
	537	c real coeff2 ! coeff for the arctan
	538	c
	539	c ! output
	540	c real time
	541	c
	542	c !local
	543	c real var
	544	c
	545	c ! 1.73205 is sqrt(3)
	546	c
	547	c var = max(
	548	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
	549	c
	550	c time =
	551	c & coeff1 *
	552	c & log( var )
	553	c & + coeff2 * 1.73205 *
	554	c & atan( (2rice+req) / (1.73205req) )
	555	c
	556	c return
	557	c end
[626]	558
	559
	560

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

Download in other formats:

Original Format