Context Navigation

improvedclouds.F @ 631

Last change on this file since 631 was 626, checked in by tnavarro, 13 years ago
New scheme for the clouds, no more sub-timestep. Clouds sedimentation is done with the dust one in callsedim.F like it was before. Added latent heat for sublimating ground ice. Bugs corrected. THIS VERSION OF THE WATER CYCLE SHOULD NOT BE USED WITH THERMALS DUE TO NEGATIVE TRACERS ISSUES.
File size: 31.6 KB

Rev	Line
[358]	1	subroutine improvedclouds(ngrid,nlay,ptimestep,
[626]	2	& pplev,pplay,pt,pdt,
[520]	3	& pq,pdq,pdqcloud,pdtcloud,
[358]	4	& nq,tauscaling,rdust,rice,nuice,
	5	& rsedcloud,rhocloud)
[520]	6	! to use 'getin'
	7	USE ioipsl_getincom
[358]	8	implicit none
	9	c------------------------------------------------------------------
	10	c This routine is used to form clouds when a parcel of the GCM is
	11	c saturated. It includes the ability to have supersaturation, a
	12	c computation of the nucleation rates, growthrates and the
	13	c scavenging of dust particles by clouds.
	14	c It is worth noting that the amount of dust is computed using the
	15	c dust optical depth computed in aeropacity.F. That's why
	16	c the variable called "tauscaling" is used to convert
	17	c pq(dust_mass) and pq(dust_number), which are relative
	18	c quantities, to absolute and realistic quantities stored in zq.
	19	c This has to be done to convert the inputs into absolute
	20	c values, but also to convert the outputs back into relative
	21	c values which are then used by the sedimentation and advection
	22	c schemes.
[626]	23	c A word about the radius growth ...
	24	c A word about nucleation and ice growth strategies ...
[358]	25
	26	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
	27	c (October 2011)
[626]	28	c T. Navarro, debug,correction, new scheme (October-April 2011)
[530]	29	c A. Spiga, optimization (February 2012)
[358]	30	c------------------------------------------------------------------
	31	#include "dimensions.h"
	32	#include "dimphys.h"
	33	#include "comcstfi.h"
	34	#include "callkeys.h"
	35	#include "tracer.h"
	36	#include "comgeomfi.h"
	37	#include "dimradmars.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 pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
	46	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
	47
[358]	48	REAL pt(ngrid,nlay) ! temperature at the middle of the
	49	! layers (K)
	50	REAL pdt(ngrid,nlay) ! tendance temperature des autres
	51	! param.
	52	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
	53	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
	54	! (kg/kg.s-1)
	55	REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust
	56	REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
	57
	58	c Outputs:
	59	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
	60	! (r_c in montmessin_2004)
	61	REAL nuice(ngrid,nlay) ! Estimated effective variance
	62	! of the size distribution
	63	REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
	64	REAL rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3)
	65	REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
	66	! H2O(kg/kg.s-1)
	67	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
	68	! a la chaleur latente
	69
	70	c------------------------------------------------------------------
	71	c Local variables:
	72
	73	LOGICAL firstcall
	74	DATA firstcall/.true./
	75	SAVE firstcall
	76
	77	REAL*8 derf ! Error function
	78	!external derf
	79
	80	REAL CBRT
	81	EXTERNAL CBRT
	82
	83	INTEGER ig,l,i
[520]	84
[358]	85
	86	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
	87	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
	88	REAL zt(ngridmx,nlayermx) ! local value of temperature
	89	REAL zqsat(ngridmx,nlayermx) ! saturation
	90	REAL lw !Latent heat of sublimation (J.kg-1)
	91	REAL Cste
	92	REAL dMice ! mass of condensated ice
	93	REAL sumcheck
	94	REAL*8 ph2o ! Water vapor partial pressure (Pa)
	95	REAL*8 satu ! Water vapor saturation ratio over ice
	96	REAL*8 Mo,No
	97	REAL*8 dN,dM,newvap
[626]	98	REAL*8 Rn, Rm, dev2, yeah, n_derf, m_derf
[358]	99	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
	100	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
	101	REAL*8 rate(nbin_cld) ! nucleation rate
[626]	102	!REAL*8 up,dwn,Ctot,gr
	103	REAl*8 seq
[358]	104	REAL*8 sig ! Water-ice/air surface tension (N.m)
	105	EXTERNAL sig
	106
	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
	120	REAL sigma_ice ! Variance of the ice and CCN distributions
	121	SAVE sigma_ice
[520]	122
[626]	123	REAL tdicho, tmax, rmin, rmax, req, rdicho
	124	REAL coeff0, coeff1, coeff2
	125	REAL error_out(ngridmx,nlayermx)
	126	REAL error2d(ngridmx)
	127
	128	REAL var1,var2,var3
	129
	130	REAL rccn, epsilon
	131
[420]	132	c----------------------------------
[626]	133	c some outputs for 1D -- TESTS
[358]	134	REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
	135	REAL dN_out(ngridmx,nlayermx) ! mass variation for output
	136	REAL dM_out(ngridmx,nlayermx) ! number variation for output
	137	REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
[626]	138	REAL gr_out(ngridmx,nlayermx) ! for 1d output
[520]	139	REAL rice_out(ngridmx,nlayermx) ! ice radius change
[455]	140	REAL rate_out(ngridmx,nlayermx) ! nucleation rate
[626]	141	REAL satubf(ngridmx,nlayermx),satuaf(ngridmx,nlayermx)
	142	REAL ccnbf(ngridmx,nlayermx),ccnaf(ngridmx,nlayermx)
[420]	143	INTEGER count
	144
[626]	145	LOGICAL test_flag ! flag for test/debuging outputs
[520]	146
[626]	147	test_flag = .false.
[420]	148	c----------------------------------
	149	c----------------------------------
[358]	150
	151	c------------------------------------------------------------------
	152
	153	IF (firstcall) THEN
[626]	154	!=============================================================
	155	! 0. Definition of the size grid
	156	!=============================================================
[358]	157	c rad_cld is the primary radius grid used for microphysics computation.
	158	c The grid spacing is computed assuming a constant volume ratio
	159	c between two consecutive bins; i.e. vrat_cld.
	160	c vrat_cld is determined from the boundary values of the size grid:
	161	c rmin_cld and rmax_cld.
	162	c The rb_cld array contains the boundary values of each rad_cld bin.
	163	c dr_cld is the width of each rad_cld bin.
	164
	165	c Volume ratio between two adjacent bins
	166	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
	167	vrat_cld = dexp(vrat_cld)
	168	write(,) "vrat_cld", vrat_cld
	169
	170	rb_cld(1) = rbmin_cld
	171	rad_cld(1) = rmin_cld
[530]	172	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
[358]	173
	174	do i=1,nbin_cld-1
[531]	175	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
[358]	176	vol_cld(i+1) = vol_cld(i) * vrat_cld
	177	enddo
	178
	179	do i=1,nbin_cld
[531]	180	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
[358]	181	& rad_cld(i)
	182	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
	183	enddo
	184	rb_cld(nbin_cld+1) = rbmax_cld
	185	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
	186
	187	print*, ' '
	188	print*,'Microphysics: size bin information:'
	189	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
	190	print*,'-----------------------------------'
	191	do i=1,nbin_cld
	192	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
	193	& dr_cld(i)
	194	enddo
	195	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
	196	print*,'-----------------------------------'
	197
[541]	198	do i=1,nbin_cld+1
	199	rb_cld(i) = dlog(rb_cld(i)) !! we save that so that it is not computed
	200	!! at each timestep and gridpoint
	201	enddo
	202
[358]	203	c Contact parameter of water ice on dust ( m=cos(theta) )
	204	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	205	! mteta = 0.95
	206	write(,) 'water_param contact parameter:', mteta
	207
	208	c Volume of a water molecule (m3)
	209	vo1 = mh2o / dble(rho_ice)
	210	c Variance of the ice and CCN distributions
	211	sigma_ice = sqrt(log(1.+nuice_sed))
	212
	213	write(,) 'Variance of ice & CCN distribs :', sigma_ice
[455]	214	write(,) 'nuice for sedimentation:', nuice_sed
[358]	215	write(,) 'Volume of a water molecule:', vo1
	216
	217	firstcall=.false.
	218	END IF
	219
[626]	220	!=============================================================
	221	! 1. Initialisation
	222	!=============================================================
	223
[411]	224	c Initialize the tendencies
	225	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	226	pdtcloud(1:ngrid,1:nlay)=0
	227
[358]	228	c Update the needed variables
	229	do l=1,nlay
	230	do ig=1,ngrid
	231	c Temperature
	232	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	233	c Dust mass mixing ratio
	234	zq(ig,l,igcm_dust_mass) =
	235	& pq(ig,l,igcm_dust_mass) +
	236	& pdq(ig,l,igcm_dust_mass) * ptimestep
	237	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
	238	c Dust particle number
	239	zq(ig,l,igcm_dust_number) =
[411]	240	& pq(ig,l,igcm_dust_number) +
	241	& pdq(ig,l,igcm_dust_number) * ptimestep
[358]	242	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
[411]	243	c Update rdust from last tendencies
	244	rdust(ig,l)=
	245	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
	246	& max(zq(ig,l,igcm_dust_number),0.01))
	247	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
[358]	248	c CCN mass mixing ratio
	249	zq(ig,l,igcm_ccn_mass)=
	250	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
	251	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
	252	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
	253	c CCN particle number
	254	zq(ig,l,igcm_ccn_number)=
	255	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
	256	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
	257	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
	258	c Water vapor
	259	zq(ig,l,igcm_h2o_vap)=
	260	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
	261	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
	262	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
	263	c Water ice
	264	zq(ig,l,igcm_h2o_ice)=
	265	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
[626]	266	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),1e-30) ! FF 12/2004
[358]	267	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
	268	enddo
	269	enddo
	270
[626]	271	!=============================================================
	272	! 2. Compute saturation
	273	!=============================================================
[358]	274
[541]	275	dev2 = 1. / ( sqrt(2.) * sigma_ice )
[626]	276	error_out(:,:) = 0.
[358]	277
[626]	278
[358]	279	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	280
[411]	281	count = 0
[358]	282
	283	c Main loop over the GCM's grid
	284	DO l=1,nlay
	285	DO ig=1,ngrid
	286
	287	c Get the partial pressure of water vapor and its saturation ratio
	288	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
	289	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
[626]	290
	291	IF (( satu .ge. 1. ) .or. ! if there is condensation
	292	& ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation
[358]	293
[411]	294
[626]	295	!=============================================================
	296	! 3. Nucleation
	297	!=============================================================
	298
[358]	299	c Expand the dust moments into a binned distribution
[626]	300	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) !+ 1.e-30
	301	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
[358]	302	Rn = rdust(ig,l)
[530]	303	Rn = -dlog(Rn)
	304	Rm = Rn - 3. * sigma_ice*sigma_ice
[626]	305	n_derf = derf( (rb_cld(1)+Rn) *dev2)
	306	m_derf = derf( (rb_cld(1)+Rm) *dev2)
[358]	307	do i = 1, nbin_cld
[626]	308	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
	309	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
	310	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
	311	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
	312	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
	313	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
[358]	314	enddo
[530]	315	!!! MORE EFFICIENT COMPUTATIONALLY THAN
	316	! Rn = rdust(ig,l)
	317	! Rm = Rn * exp( 3. * sigma_ice**2. )
	318	! Rn = 1. / Rn
	319	! Rm = 1. / Rm
	320	! do i = 1, nbin_cld
	321	! n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
	322	! & -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
	323	! m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
	324	! & -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
	325	! enddo
	326	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	327
[358]	328	! sumcheck = 0
	329	! do i = 1, nbin_cld
	330	! sumcheck = sumcheck + n_aer(i)
	331	! enddo
	332	! sumcheck = abs(sumcheck/No - 1)
	333	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	334	! print*, "WARNING, No sumcheck PROBLEM"
	335	! print*, "sumcheck, No",sumcheck, No
	336	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	337	! print*, "Dust binned distribution", n_aer
	338	! endif
	339	!
	340	! sumcheck = 0
	341	! do i = 1, nbin_cld
[411]	342	! sumcheck = sumcheck + m_aer(i)
[358]	343	! enddo
	344	! sumcheck = abs(sumcheck/Mo - 1)
	345	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
	346	! print*, "WARNING, Mo sumcheck PROBLEM"
[411]	347	! print*, "sumcheck, Mo",sumcheck, Mo
[358]	348	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
	349	! print*, "Dust binned distribution", m_aer
	350	! endif
	351
	352	c Get the rates of nucleation
	353	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
[411]	354
[358]	355
	356	dN = 0.
	357	dM = 0.
[455]	358	rate_out(ig,l) = 0
[358]	359	do i = 1, nbin_cld
	360	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
	361	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
	362	dN = dN + n_aer(i) * rate(i) * ptimestep
	363	dM = dM + m_aer(i) * rate(i) * ptimestep
[520]	364	!rate_out(ig,l)=rate_out(ig,l)+rate(i)
[358]	365	enddo
	366
[626]	367	c Update CCNs, can also be done after the radius growth
	368	c Dust particles
	369	zq(ig,l,igcm_dust_mass) =
	370	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig)
	371	zq(ig,l,igcm_dust_number) =
	372	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
	373	c CCNs
	374	zq(ig,l,igcm_ccn_mass) =
	375	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig)
	376	zq(ig,l,igcm_ccn_number) =
	377	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
	378
	379
	380	!=============================================================
	381	! 4. Ice growth: scheme for radius evolution
	382	!=============================================================
	383
	384	Mo = zq(ig,l,igcm_h2o_ice) +
	385	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
	386	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
	387
	388	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
	389	& + zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
	390	& / Mo * rho_dust
	391	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
	392
	393	c nuclei radius
	394	rccn = CBRT(zq(ig,l,igcm_ccn_mass)/
	395	& (pirho_dustzq(ig,l,igcm_ccn_number)*4./3.))
	396
	397	c ice crystal radius
	398	rice (ig,l) =
	399	& CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) )
	400
	401	c enforce physical value : crystal cannot be smaller than its nuclei !
	402	rice(ig,l) = max(rice(ig,l), rccn)
	403
	404	c saturation at equilibrium
	405	seq = exp( 2.sig(zt(ig,l))mh2o /
[358]	406	& (rho_icergpzt(ig,l)*rice(ig,l)) )
	407
	408
[626]	409	c If there is more than on nuclei, we peform ice growth
	410	var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
	411	IF (var1 .ge. -1) THEN
	412
	413
	414	if (test_flag) then
	415	print*, ' '
	416	print*, ptimestep
	417	print*, 'satu,seq', real(satu), real(seq), ig,l
	418	print*, 'dN,dM', real(dN),real(dM)
	419	print*,'rccn', rccn
	420	print, 'Nccn, Mccn', zq(ig,l,igcm_ccn_number)tauscaling(ig),
	421	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
	422	endif
	423
	424	c crystal radius to reach saturation at equilibrium (i.e. satu=seq)
	425	req = ( zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
	426	& + zq(ig,l,igcm_ccn_mass)tauscaling(ig)
	427	& rho_ice/rho_dust - seq * zqsat(ig,l))
	428	& / ( zq(ig,l,igcm_ccn_number)tauscaling(ig)
	429	& pirho_ice4./3. )
	430	req = CBRT(req)
[358]	431
[626]	432	c compute ice radius growth resistances (diffusive and latent heat resistancea)
	433	call growthrate(zt(ig,l),pplay(ig,l),
	434	& ph2o/satu,seq,req,coeff1,coeff2)
	435
	436	coeff0 = -zqsat(ig,l) / (4.pireq*rho_ice
	437	& * zq(ig,l,igcm_ccn_number)*tauscaling(ig))
	438
	439	c compute tmax, the time needed to reach req
	440	call phi(rice(ig,l),req,coeff1,coeff2,tmax)
	441
	442	if (test_flag) then
	443	print*, 'coeffs', coeff0,coeff1,coeff2
	444	print, 'req,tmax', req,tmaxcoeff0
	445	print*, 'i,rmin,rdicho,rmax,tdicho'
	446	endif
[358]	447
[626]	448	c rmin is rice if r increases (satu >1) or req if it decreases (satu<1)
	449	c if req is lower than rccn (ie not enough ice to reach saturation), rmin is forced to rccn
	450	if (satu .ge. seq) then
	451	! crystal size is increasing
	452	rmin = max(min(rice(ig,l),req),rccn)
	453	rmax = max(rice(ig,l),req)
	454	else
	455	! crystal size is decreasing
	456	rmax = max(min(rice(ig,l),req),rccn)
	457	rmin = max(rice(ig,l),req)
	458	endif
	459	!rmax = min(rmax,1.e-3) ! au max on a des rayons de 1 mm pour la dicho ...
	460	rdicho = 0.5*(rmin+rmax)
[358]	461
[626]	462	! for output
	463	var1 = rice(ig,l)
	464	var2 = rmin
	465	var3 = rmax
	466
	467	c Given the phi function is monotonous, we perform a simple dichotomy to find the raidus at t+1
	468	do i = 1,10 ! dichotomy loop
	469
	470	c compute tdicho, the time needed to reach rdicho
	471	call phi(rdicho,req,coeff1,coeff2,tdicho)
	472	!print*, tdicho,tmax
	473	tdicho = coeff0*(tdicho - tmax)
	474
	475	if (test_flag) print*, i,rmin,rdicho,rmax,tdicho
	476
	477	if (tdicho .ge. ptimestep) then
	478	rmax = rdicho
	479	else
	480	rmin = rdicho
	481	endif
	482
	483	rdicho = 0.5*(rmin+rmax)
	484
	485	enddo ! of dichotomy loop
	486
	487	if (test_flag) then
	488	if (abs(rdicho - rccn) .ge. 1e-15) then ! to avoid infinite values
	489	epsilon = (rmax - rmin)/(2**10)
	490	error_out(ig,l) =
	491	& 100(epsilon3 +3epsilon*2rdicho +3epsilonrdicho**2)
	492	& / (rdicho3-rccn3)
	493	endif
	494	print*, 'error masse glace %', error_out(ig,l)
	495	print*, 'rice,ice,vap bf',
	496	& rice(ig,l),zq0(ig,l,igcm_h2o_ice),zq0(ig,l,igcm_h2o_vap)
	497	endif
	498
	499	c If the initial condition is subsaturated and there is not enough ice available for sublimation
	500	c to reach equilibrium, req is neagtive. Therefore, enforce physical value.
	501	rice(ig,l) = max(rdicho,rccn)
	502
	503	!!!!! Water ice mass is computed with radius at t+1 and their current number
	504	!!!!! Nccn is at t or t+1, depending on what has been done before
	505	! zq(ig,l,igcm_h2o_ice) =
	506	! & (pirho_icezq(ig,l,igcm_ccn_number)*4./3.
	507	! & * rdichordichordicho -
	508	! & zq(ig,l,igcm_ccn_mass)*rho_ice/rho_dust)
	509	! & * tauscaling(ig)
	510
	511	!!!!! Water ice mass is computed with radius at t+1 and their number at t+1
	512	!!!!! that is dirty, but this enforces the use of Nccn at t+1, whatever is done before
	513	!!!!! TO BE CLEANED ONE DAY
	514	var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
	515	var2 = zq0(ig,l,igcm_ccn_mass)*tauscaling(ig) + dM
	516	zq(ig,l,igcm_h2o_ice) =
	517	& (pirho_icevar1*4./3.
	518	& * rdichordichordicho -
	519	& var2*rho_ice/rho_dust)
	520
[358]	521
[626]	522	!!!! enforce realistic values (if the case of growth on Nccn(t) and condensation on Nccn(t+1))
	523	zq(ig,l,igcm_h2o_ice) =
	524	& min(zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap),
	525	& zq(ig,l,igcm_h2o_ice))
	526	zq(ig,l,igcm_h2o_ice) =
	527	& max(1e-30,zq(ig,l,igcm_h2o_ice))
	528
	529	zq(ig,l,igcm_h2o_vap) =
	530	& zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
	531	& - zq(ig,l,igcm_h2o_ice)
	532
	533
	534	if (test_flag) then
	535	print*, 'rice,ice,vap af',
	536	& rice(ig,l),zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_h2o_vap)
	537	print*, 'satu bf, af',
	538	& zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l),
	539	& zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
	540	endif
	541
	542
	543	!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
	544	if ((zq(ig,l,igcm_h2o_ice).le. -1e-8)
	545	& .or. (zq(ig,l,igcm_h2o_vap).le. -1e-8)) then
	546	print*, 'NEGATIVE WATER'
	547	print*, 'ig,l', ig,l
	548	print*, 'satu', satu
	549	print*, 'vap, ice bf',
	550	& zq0(ig,l,igcm_h2o_vap), zq0(ig,l,igcm_h2o_ice)
	551	print*, 'vap, ice af',
	552	& zq(ig,l,igcm_h2o_vap), zq(ig,l,igcm_h2o_ice)
	553	print*, 'ccn N,M bf',
	554	& zq0(ig,l,igcm_ccn_number), zq0(ig,l,igcm_ccn_mass)
	555	print*, 'ccn N,M af',
	556	& zq(ig,l,igcm_ccn_number), zq(ig,l,igcm_ccn_mass)
	557	print*, 'tauscaling',
	558	& tauscaling(ig)
	559	print*, 'req,rccn,rice bf,rice af',
	560	& req,rccn,var1,rice(ig,l)
	561	print*, 'rmin, rmax', var2,var3
	562	print*, 'error_out,tdicho,ptimestep',
	563	& error_out(ig,l),tdicho,ptimestep
	564	print*, ' '
	565	endif
	566	!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
[358]	567
[626]	568
	569	ENDIF !of if Nccn >1
	570
[358]	571
[626]	572	!=============================================================
	573	! 5. Dust cores released, tendancies, latent heat, etc ...
	574	!=============================================================
	575
[358]	576	c If all the ice particles sublimate, all the condensation
[626]	577	c nuclei are released:
	578	if (zq(ig,l,igcm_h2o_ice).le.1e-10) then
	579	! for coherence
	580	! dM = 0
	581	! dN = 0
	582	! dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
	583	! dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
	584	c Water
	585	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
	586	& + zq(ig,l,igcm_h2o_ice)
[358]	587	zq(ig,l,igcm_h2o_ice) = 0.
	588	c Dust particles
[626]	589	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
	590	& + zq(ig,l,igcm_ccn_mass)
	591	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
	592	& + zq(ig,l,igcm_ccn_number)
[358]	593	c CCNs
	594	zq(ig,l,igcm_ccn_mass) = 0.
	595	zq(ig,l,igcm_ccn_number) = 0.
	596	endif
[626]	597
[358]	598
[626]	599	! dN = dN/ tauscaling(ig)
	600	! dM = dM/ tauscaling(ig)
	601	!c Dust particles
	602	! zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
	603	! zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
	604	!c CCNs
	605	! zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
	606	! zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
[520]	607
	608
[411]	609	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
	610	& -zq0(ig,l,igcm_dust_mass))/ptimestep
	611	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
	612	& -zq0(ig,l,igcm_dust_number))/ptimestep
	613	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
	614	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
	615	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
	616	& -zq0(ig,l,igcm_ccn_number))/ptimestep
	617	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
	618	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
	619	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
	620	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
	621
[520]	622
[411]	623	count = count +1
[626]	624
	625
[411]	626	ELSE ! for coherence (rdust, rice computations etc ..)
	627	zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass)
	628	zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
	629	zq(ig,l,igcm_ccn_mass) = zq0(ig,l,igcm_ccn_mass)
	630	zq(ig,l,igcm_ccn_number) = zq0(ig,l,igcm_ccn_number)
	631	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice)
	632	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap)
[626]	633
[411]	634	! pour les sorties de test
[626]	635	! satu_out(ig,l) = satu
	636	! gr_out(ig,l) = 0
	637	! dN_out(ig,l) = 0
	638	! dM_out(ig,l) = 0
[411]	639
	640	ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
	641
[626]	642	! ccnbf(ig,l) = zq0(ig,l,igcm_ccn_number) * tauscaling(ig)
	643	! ccnaf(ig,l) = zq(ig,l,igcm_ccn_number) * tauscaling(ig)
	644	!
	645	! satubf(ig,l) = zq0(ig,l,igcm_h2o_vap) / zqsat(ig,l)
	646	! satuaf(ig,l) = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
	647
	648
	649	c-----update temperature (latent heat relase)
[411]	650	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
[626]	651	pdtcloud(ig,l)= -pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
[411]	652
[520]	653	c----- update rice & rhocloud AFTER microphysic
[411]	654	Mo = zq(ig,l,igcm_h2o_ice) +
	655	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
	656	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
	657	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
	658	& +zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
	659	& / Mo * rho_dust
	660	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
[520]	661
[541]	662	if ((Mo.lt.1.e-20) .or. (No.le.1)) then
	663	rice(ig,l) = 1.e-8
	664	else
	665	rice(ig,l) =
	666	& CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
	667	endif
[411]	668
	669	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
[358]	670
[411]	671	c----- update rdust and sedimentation radius
[358]	672	rdust(ig,l)=
	673	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
	674	& max(zq(ig,l,igcm_dust_number),0.01))
	675	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
[411]	676
[530]	677	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)
	678	& (1.+nuice_sed)*(1.+nuice_sed),
[358]	679	& rdust(ig,l) )
	680	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	681
[626]	682	ENDDO ! of ig loop
	683	ENDDO ! of nlayer loop
[520]	684
[411]	685
[626]	686	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	687	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	688	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	689	IF (test_flag) then
	690
	691	error2d(:) = 0.
	692	DO l=1,nlay
	693	DO ig=1,ngrid
	694	error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
	695	ENDDO
	696	ENDDO
[420]	697
[411]	698	print*, 'count is ',count, ' i.e. ',
	699	& count100/(nlayngrid), '% for microphys computation'
[358]	700
[420]	701	IF (ngrid.ne.1) THEN ! 3D
[626]	702	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
	703	! & satu_out)
	704	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
	705	! & dM_out)
	706	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
	707	! & dN_out)
	708	call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
	709	& error2d)
	710	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
	711	! & zqsat)
[411]	712	ENDIF
[358]	713
[411]	714	IF (ngrid.eq.1) THEN ! 1D
[626]	715	call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
	716	& error_out)
	717	call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
	718	& satubf)
	719	call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
	720	& satuaf)
	721	call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
	722	& zq0(1,:,igcm_h2o_vap))
	723	call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
	724	& zq(1,:,igcm_h2o_vap))
	725	call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
	726	& zq0(1,:,igcm_h2o_ice))
	727	call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
	728	& zq(1,:,igcm_h2o_ice))
	729	call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
	730	& ccnbf)
	731	call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
	732	& ccnaf)
	733	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
	734	c & gr_out)
	735	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
	736	c & rate_out)
	737	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
	738	c & dM_out)
	739	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
	740	c & dN_out)
[358]	741	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
	742	& zqsat)
	743	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
	744	& satu_out)
[411]	745	call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
[358]	746	& rice)
[420]	747	call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
[358]	748	& rdust)
	749	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
	750	& rsedcloud)
	751	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
	752	& rhocloud)
	753	ENDIF
[626]	754
	755	ENDIF ! endif test_flag
	756	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	757	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	758	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
	759
[358]	760	return
	761	end
[626]	762
	763
	764
	765	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	766	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	767	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
	768	c It is an analytical solution to the ice radius growth equation,
	769	c with the approximation of a constant 'reduced' cunningham correction factor
	770	c (lambda in growthrate.F) taken at radius req instead of rice
	771	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	772	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	773
	774	subroutine phi(rice,req,coeff1,coeff2,time)
	775
	776	implicit none
	777
	778	! inputs
	779	real rice ! ice radius
	780	real req ! ice radius at equilibirum
	781	real coeff1 ! coeff for the log
	782	real coeff2 ! coeff for the arctan
	783
	784	! output
	785	real time
	786
	787	!local
	788	real var
	789
	790	! 1.73205 is sqrt(3)
	791
	792	var = max(
	793	& abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
	794
	795	time =
	796	& coeff1 *
	797	& log( var )
	798	& + coeff2 * 1.73205 *
	799	& atan( (2rice+req) / (1.73205req) )
	800
	801	return
	802	end
	803
	804
	805

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

Download in other formats:

Original Format