Context Navigation

improvedclouds.F @ 1242

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

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

Download in other formats:

Original Format