Context Navigation

improvedclouds_mod.F @ 2156

Last change on this file since 2156 was 2151, checked in by aslmd, 5 years ago
changed old functions dexp dlog in generic exp and log
File size: 21.5 KB

Rev	Line
[1963]	1	MODULE improvedclouds_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[1909]	7	subroutine improvedclouds(ngrid,nlay,microtimestep,
	8	& pplay,pteff,sum_subpdt,
	9	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
[633]	10	& nq,tauscaling)
[520]	11	! to use 'getin'
	12	USE ioipsl_getincom
[740]	13	USE updaterad
[1996]	14	USE watersat_mod, ONLY: watersat
[1036]	15	use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
	16	& igcm_h2o_ice, igcm_dust_mass,
	17	& igcm_dust_number, igcm_ccn_mass,
	18	& igcm_ccn_number
[1047]	19	use conc_mod, only: mmean
[1226]	20	USE comcstfi_h
[358]	21	implicit none
[633]	22
	23
[358]	24	c------------------------------------------------------------------
	25	c This routine is used to form clouds when a parcel of the GCM is
	26	c saturated. It includes the ability to have supersaturation, a
	27	c computation of the nucleation rates, growthrates and the
	28	c scavenging of dust particles by clouds.
	29	c It is worth noting that the amount of dust is computed using the
	30	c dust optical depth computed in aeropacity.F. That's why
	31	c the variable called "tauscaling" is used to convert
	32	c pq(dust_mass) and pq(dust_number), which are relative
	33	c quantities, to absolute and realistic quantities stored in zq.
	34	c This has to be done to convert the inputs into absolute
	35	c values, but also to convert the outputs back into relative
	36	c values which are then used by the sedimentation and advection
	37	c schemes.
	38
	39	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
	40	c (October 2011)
[626]	41	c T. Navarro, debug,correction, new scheme (October-April 2011)
[530]	42	c A. Spiga, optimization (February 2012)
[358]	43	c------------------------------------------------------------------
	44	#include "callkeys.h"
	45	#include "microphys.h"
	46	c------------------------------------------------------------------
[1976]	47	c Inputs/outputs:
[358]	48
[1976]	49	INTEGER, INTENT(IN) :: ngrid,nlay
	50	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
	51	REAL, INTENT(IN) :: microtimestep ! pas de temps physique (s)
	52	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
	53	REAL, INTENT(IN) :: pteff(ngrid,nlay) ! temperature at the middle of the
	54	! layers (K)
	55	REAL, INTENT(IN) :: sum_subpdt(ngrid,nlay)! tendance temperature des autres
	56	! param.
	57	REAL, INTENT(IN) :: pqeff(ngrid,nlay,nq) ! traceur (kg/kg)
	58	REAL, INTENT(IN) :: sum_subpdq(ngrid,nlay,nq) ! tendance avant condensation
	59	! (kg/kg.s-1)
	60	REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
	61
	62	REAL, INTENT(OUT) :: subpdqcloud(ngrid,nlay,nq) ! tendance de la condensation
	63	! H2O(kg/kg.s-1)
	64	REAL, INTENT(OUT) :: subpdtcloud(ngrid,nlay) ! tendance temperature due
	65	! a la chaleur latente
[358]	66
	67	c------------------------------------------------------------------
	68	c Local variables:
	69
	70	LOGICAL firstcall
	71	DATA firstcall/.true./
	72	SAVE firstcall
	73
	74	REAL*8 derf ! Error function
	75	!external derf
[740]	76
[358]	77	INTEGER ig,l,i
[520]	78
[1047]	79	REAL zq(ngrid,nlay,nq) ! local value of tracers
	80	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
	81	REAL zt(ngrid,nlay) ! local value of temperature
	82	REAL zqsat(ngrid,nlay) ! saturation
[358]	83	REAL lw !Latent heat of sublimation (J.kg-1)
[633]	84	REAL cste
	85	REAL dMice ! mass of condensed ice
	86	! REAL sumcheck
[358]	87	REAL*8 ph2o ! Water vapor partial pressure (Pa)
	88	REAL*8 satu ! Water vapor saturation ratio over ice
	89	REAL*8 Mo,No
[633]	90	REAL*8 Rn, Rm, dev2, n_derf, m_derf
[358]	91	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
	92	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
[633]	93
[358]	94	REAL*8 sig ! Water-ice/air surface tension (N.m)
	95	EXTERNAL sig
	96
[633]	97	REAL dN,dM
	98	REAL rate(nbin_cld) ! nucleation rate
	99	REAL seq
	100
	101	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
	102	! (r_c in montmessin_2004)
[1047]	103	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
	104	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
[633]	105
	106	REAL res ! Resistance growth
[740]	107
[633]	108
[358]	109	c Parameters of the size discretization
	110	c used by the microphysical scheme
	111	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
	112	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
	113	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
	114	! Minimum boundary radius (m)
	115	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
	116	DOUBLE PRECISION vrat_cld ! Volume ratio
	117	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
	118	SAVE rb_cld
[520]	119	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
	120	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
[358]	121
[633]	122
[358]	123	REAL sigma_ice ! Variance of the ice and CCN distributions
	124	SAVE sigma_ice
[633]	125
	126
[520]	127
[420]	128	c----------------------------------
[633]	129	c TESTS
	130
	131	INTEGER countcells
[420]	132
[626]	133	LOGICAL test_flag ! flag for test/debuging outputs
[740]	134	SAVE test_flag
	135
	136
	137	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
	138	REAL res_out(ngrid,nlay)
[633]	139
[358]	140
	141	c------------------------------------------------------------------
	142
[1779]	143	! AS: firstcall OK absolute
[358]	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
[1307]	157	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
	158	! vrat_cld = exp(vrat_cld)
[2151]	159	vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
	160	vrat_cld = exp(vrat_cld)
[358]	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
[1307]	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))
[2151]	194	rb_cld(i) = log(rb_cld(i)) !! we save that so that it is not computed
[541]	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	!=============================================================
[1909]	223	cste = 4pirho_ice*microtimestep
[626]	224
[740]	225	res_out(:,:) = 0
	226	rice(:,:) = 1.e-8
	227
[411]	228	c Initialize the tendencies
[1909]	229	subpdqcloud(1:ngrid,1:nlay,1:nq)=0
	230	subpdtcloud(1:ngrid,1:nlay)=0
[633]	231
[411]	232
[633]	233	zt(1:ngrid,1:nlay) =
[1909]	234	& pteff(1:ngrid,1:nlay) +
	235	& sum_subpdt(1:ngrid,1:nlay) * microtimestep
[358]	236
[633]	237	zq(1:ngrid,1:nlay,1:nq) =
[1909]	238	& pqeff(1:ngrid,1:nlay,1:nq) +
	239	& sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
[633]	240
	241
	242	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
	243	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
[768]	244
	245	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
[633]	246
[626]	247	!=============================================================
	248	! 2. Compute saturation
	249	!=============================================================
[358]	250
[633]	251
[541]	252	dev2 = 1. / ( sqrt(2.) * sigma_ice )
[358]	253
[1047]	254	call watersat(ngrid*nlay,zt,pplay,zqsat)
[358]	255
[633]	256	countcells = 0
[358]	257
	258	c Main loop over the GCM's grid
	259	DO l=1,nlay
	260	DO ig=1,ngrid
	261
	262	c Get the partial pressure of water vapor and its saturation ratio
[663]	263	ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
[689]	264	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
[358]	265
[626]	266	!=============================================================
	267	! 3. Nucleation
	268	!=============================================================
	269
[633]	270	IF ( satu .ge. 1. ) THEN ! if there is condensation
	271
[740]	272	call updaterccn(zq(ig,l,igcm_dust_mass),
	273	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
[633]	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)
[2151]	280	Rn = -log(Rn)
[530]	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
[1909]	324	n_aer(i) = n_aer(i)/( 1. + rate(i)*microtimestep)
	325	m_aer(i) = m_aer(i)/( 1. + rate(i)*microtimestep)
	326	dN = dN + n_aer(i) * rate(i) * microtimestep
	327	dM = dM + m_aer(i) * rate(i) * microtimestep
[358]	328	enddo
	329
[633]	330
	331	c Update Dust particles
[626]	332	zq(ig,l,igcm_dust_mass) =
[740]	333	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	334	zq(ig,l,igcm_dust_number) =
[740]	335	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[633]	336	c Update CCNs
[626]	337	zq(ig,l,igcm_ccn_mass) =
[740]	338	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	339	zq(ig,l,igcm_ccn_number) =
[740]	340	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
[626]	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
[740]	354
	355	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
	356	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
	357	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
[633]	358
	359	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
[626]	360
	361	c saturation at equilibrium
[740]	362	c rice should not be too small, otherwise seq value is not valid
	363	seq = exp(2.sig(zt(ig,l))mh2o / (rho_icergpzt(ig,l)*
	364	& max(rice(ig,l),1.e-7)))
	365
[633]	366	c get resistance growth
	367	call growthrate(zt(ig,l),pplay(ig,l),
	368	& real(ph2o/satu),rice(ig,l),res)
[358]	369
[740]	370	res_out(ig,l) = res
[626]	371
[633]	372	ccccccc implicit scheme of mass growth
[626]	373
[633]	374	dMice =
	375	& (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
	376	& /(reszqsat(ig,l)/(csteNo*rice(ig,l)) + 1.)
[358]	377
[626]	378
[633]	379	! With the above scheme, dMice cannot be bigger than vapor,
	380	! but can be bigger than all available ice.
	381	dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
	382	dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
	383
	384	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
	385	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
	386
	387
	388	countcells = countcells + 1
	389
	390	! latent heat release
	391	lw=(2834.3-0.28*(zt(ig,l)-To)-
	392	& 0.004(zt(ig,l)-To)(zt(ig,l)-To))*1.e+3
[1909]	393	subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
[358]	394
[626]	395
[358]	396
[626]	397	!=============================================================
	398	! 5. Dust cores released, tendancies, latent heat, etc ...
	399	!=============================================================
	400
[358]	401	c If all the ice particles sublimate, all the condensation
[626]	402	c nuclei are released:
[1307]	403	if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
[633]	404
[626]	405	c Water
	406	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
	407	& + zq(ig,l,igcm_h2o_ice)
[358]	408	zq(ig,l,igcm_h2o_ice) = 0.
	409	c Dust particles
[626]	410	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
	411	& + zq(ig,l,igcm_ccn_mass)
	412	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
	413	& + zq(ig,l,igcm_ccn_number)
[358]	414	c CCNs
	415	zq(ig,l,igcm_ccn_mass) = 0.
	416	zq(ig,l,igcm_ccn_number) = 0.
[633]	417
[358]	418	endif
[411]	419
[633]	420	ENDIF !of if Nccn>1
	421
[626]	422	ENDDO ! of ig loop
	423	ENDDO ! of nlayer loop
[520]	424
[633]	425
	426	! Get cloud tendencies
[1909]	427	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
[740]	428	& (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
[1909]	429	& zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
	430	subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
[740]	431	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
[1909]	432	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
	433	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
[740]	434	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
[1909]	435	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
	436	subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
[740]	437	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
[1909]	438	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
[411]	439
[740]	440	if (scavenging) then
[633]	441
[1909]	442	subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
[740]	443	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
[1909]	444	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
	445	subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
[740]	446	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
[1909]	447	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
[740]	448
	449	endif
[633]	450
[740]	451
	452
[626]	453	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	454	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
[740]	455	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
[626]	456	IF (test_flag) then
	457
[633]	458	! error2d(:) = 0.
[740]	459	DO l=1,nlay
	460	DO ig=1,ngrid
[633]	461	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
[740]	462	satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
	463	satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
	464	ENDDO
	465	ENDDO
[420]	466
[633]	467	print*, 'count is ',countcells, ' i.e. ',
	468	& countcells100/(nlayngrid), '% for microphys computation'
[358]	469
[1212]	470	#ifndef MESOSCALE
[633]	471	! IF (ngrid.ne.1) THEN ! 3D
[626]	472	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
	473	! & satu_out)
	474	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
	475	! & dM_out)
	476	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
	477	! & dN_out)
[633]	478	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
	479	! & error2d)
[626]	480	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
	481	! & zqsat)
[633]	482	! ENDIF
[358]	483
[633]	484	! IF (ngrid.eq.1) THEN ! 1D
	485	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
	486	! & error_out)
[740]	487	call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
	488	& res_out)
	489	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
	490	& satubf)
	491	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
	492	& satuaf)
	493	call WRITEdiagfi(ngrid,"vapbf","h2ovap before","kg/kg",1,
[1020]	494	& zq0(1,1,igcm_h2o_vap))
[740]	495	call WRITEdiagfi(ngrid,"vapaf","h2ovap after","kg/kg",1,
[1020]	496	& zq(1,1,igcm_h2o_vap))
[740]	497	call WRITEdiagfi(ngrid,"icebf","h2oice before","kg/kg",1,
[1020]	498	& zq0(1,1,igcm_h2o_ice))
[740]	499	call WRITEdiagfi(ngrid,"iceaf","h2oice after","kg/kg",1,
[1020]	500	& zq(1,1,igcm_h2o_ice))
[740]	501	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
[1020]	502	& zq0(1,1,igcm_ccn_number))
[740]	503	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
[1020]	504	& zq(1,1,igcm_ccn_number))
[626]	505	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
	506	c & gr_out)
	507	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
	508	c & rate_out)
	509	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
	510	c & dM_out)
	511	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
	512	c & dN_out)
[740]	513	call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
	514	& zqsat)
[633]	515	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
	516	! & satu_out)
[740]	517	call WRITEdiagfi(ngrid,"rice","ice radius","m",1,
	518	& rice)
[633]	519	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
	520	! & rdust)
	521	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
	522	! & rsedcloud)
	523	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
	524	! & rhocloud)
	525	! ENDIF
[1212]	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
[1963]	534
[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
[1963]	579	END SUBROUTINE improvedclouds
	580
	581	END MODULE improvedclouds_mod

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

Download in other formats:

Original Format