Context Navigation

improvedclouds.F @ 658

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

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

Download in other formats:

Original Format