Context Navigation

improvedclouds.F @ 1884

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

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

Download in other formats:

Original Format