Context Navigation

improvedclouds.F @ 517

Last change on this file since 517 was 455, checked in by tnavarro, 14 years ago
tiny change : nuice_sed initialisation is now done in inifis. Also changed initracer and improvedclouds.
File size: 23.5 KB

Rev	Line
[358]	1	subroutine improvedclouds(ngrid,nlay,ptimestep,
[411]	2	& pplev,pplay,pt,pdt,
[358]	3	& pq,pdq,pdqcloud,pdqscloud,pdtcloud,
	4	& nq,tauscaling,rdust,rice,nuice,
	5	& rsedcloud,rhocloud)
	6	implicit none
	7	c------------------------------------------------------------------
	8	c This routine is used to form clouds when a parcel of the GCM is
	9	c saturated. It includes the ability to have supersaturation, a
	10	c computation of the nucleation rates, growthrates and the
	11	c scavenging of dust particles by clouds.
	12	c It is worth noting that the amount of dust is computed using the
	13	c dust optical depth computed in aeropacity.F. That's why
	14	c the variable called "tauscaling" is used to convert
	15	c pq(dust_mass) and pq(dust_number), which are relative
	16	c quantities, to absolute and realistic quantities stored in zq.
	17	c This has to be done to convert the inputs into absolute
	18	c values, but also to convert the outputs back into relative
	19	c values which are then used by the sedimentation and advection
	20	c schemes.
	21
	22	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
	23	c (October 2011)
	24	c------------------------------------------------------------------
	25	#include "dimensions.h"
	26	#include "dimphys.h"
	27	#include "comcstfi.h"
	28	#include "callkeys.h"
	29	#include "tracer.h"
	30	#include "comgeomfi.h"
	31	#include "dimradmars.h"
	32	#include "microphys.h"
	33	c------------------------------------------------------------------
	34	c Inputs:
	35
	36	INTEGER ngrid,nlay
	37	integer nq ! nombre de traceurs
	38	REAL ptimestep ! pas de temps physique (s)
[411]	39	REAL pplev(ngrid,nlay+1),pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
[358]	40	REAL pt(ngrid,nlay) ! temperature at the middle of the
	41	! layers (K)
	42	REAL pdt(ngrid,nlay) ! tendance temperature des autres
	43	! param.
	44	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
	45	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
	46	! (kg/kg.s-1)
	47	REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust
	48	REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
	49
	50	c Outputs:
	51	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
	52	! (r_c in montmessin_2004)
	53	REAL nuice(ngrid,nlay) ! Estimated effective variance
	54	! of the size distribution
	55	REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
	56	REAL rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3)
	57	REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
	58	! H2O(kg/kg.s-1)
	59	REAL pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
	60	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
	61	! a la chaleur latente
	62
	63	c------------------------------------------------------------------
	64	c Local variables:
	65
	66	LOGICAL firstcall
	67	DATA firstcall/.true./
	68	SAVE firstcall
	69
	70	REAL*8 derf ! Error function
	71	!external derf
	72
	73	REAL CBRT
	74	EXTERNAL CBRT
	75
	76	INTEGER ig,l,i
	77
	78	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
	79	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
	80	REAL zt(ngridmx,nlayermx) ! local value of temperature
	81	REAL zqsat(ngridmx,nlayermx) ! saturation
	82	REAL lw !Latent heat of sublimation (J.kg-1)
	83	REAL Cste
	84	REAL dMice ! mass of condensated ice
	85	REAL sumcheck
	86	REAL*8 ph2o ! Water vapor partial pressure (Pa)
	87	REAL*8 satu ! Water vapor saturation ratio over ice
	88	REAL*8 Mo,No
	89	REAL*8 dN,dM,newvap
	90	REAL*8 Rn, Rm, dev2
	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
	93	REAL*8 rate(nbin_cld) ! nucleation rate
	94	REAL*8 up,dwn,Ctot,gr,seq
	95	REAL*8 sig ! Water-ice/air surface tension (N.m)
	96	EXTERNAL sig
	97
	98	c Parameters of the size discretization
	99	c used by the microphysical scheme
	100	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
	101	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
	102	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
	103	! Minimum boundary radius (m)
	104	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
	105	DOUBLE PRECISION vrat_cld ! Volume ratio
	106	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
	107	SAVE rb_cld
	108	DOUBLE PRECISION dr_cld(nbin_cld)! width of each rad_cld bin (m)
	109	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
	110
	111	REAL sigma_ice ! Variance of the ice and CCN distributions
	112	SAVE sigma_ice
[420]	113
	114	c----------------------------------
	115	c some outputs for 1D -- TEMPORARY
[358]	116	REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
	117	REAL dN_out(ngridmx,nlayermx) ! mass variation for output
	118	REAL dM_out(ngridmx,nlayermx) ! number variation for output
[411]	119	REAL dM_col(ngridmx) ! total mass condensed in column
	120	REAL dN_col(ngridmx) ! total mass condensed in column
[358]	121	REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
	122	REAL gr_out(ngridmx,nlayermx) ! for 1d output
	123	REAL newvap_out(ngridmx,nlayermx) ! for 1d output
[411]	124	REAL Mdust_col(ngridmx) ! total column dust mass
	125	REAL Ndust_col(ngridmx) ! total column dust mass
	126	REAL Mccn_col(ngridmx) ! total column ccn mass
	127	REAL Nccn_col(ngridmx) ! total column ccn mass
[455]	128	REAL rate_out(ngridmx,nlayermx) ! nucleation rate
[420]	129	INTEGER count
	130
	131	LOGICAL output_sca ! scavenging outputs flag for tests
	132	output_sca = .false.
	133	c----------------------------------
	134	c----------------------------------
[358]	135
	136	c------------------------------------------------------------------
	137
	138	IF (firstcall) THEN
	139
	140	c Definition of the size grid
	141	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~
	142	c rad_cld is the primary radius grid used for microphysics computation.
	143	c The grid spacing is computed assuming a constant volume ratio
	144	c between two consecutive bins; i.e. vrat_cld.
	145	c vrat_cld is determined from the boundary values of the size grid:
	146	c rmin_cld and rmax_cld.
	147	c The rb_cld array contains the boundary values of each rad_cld bin.
	148	c dr_cld is the width of each rad_cld bin.
	149
	150	c Volume ratio between two adjacent bins
	151	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
	152	vrat_cld = dexp(vrat_cld)
	153	write(,) "vrat_cld", vrat_cld
	154
	155	rb_cld(1) = rbmin_cld
	156	rad_cld(1) = rmin_cld
	157	vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3.
	158
	159	do i=1,nbin_cld-1
	160	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
	161	vol_cld(i+1) = vol_cld(i) * vrat_cld
	162	enddo
	163
	164	do i=1,nbin_cld
	165	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
	166	& rad_cld(i)
	167	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
	168	enddo
	169	rb_cld(nbin_cld+1) = rbmax_cld
	170	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
	171
	172	print*, ' '
	173	print*,'Microphysics: size bin information:'
	174	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
	175	print*,'-----------------------------------'
	176	do i=1,nbin_cld
	177	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
	178	& dr_cld(i)
	179	enddo
	180	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
	181	print*,'-----------------------------------'
	182
	183	c Contact parameter of water ice on dust ( m=cos(theta) )
	184	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	185	! mteta = 0.95
	186	write(,) 'water_param contact parameter:', mteta
	187
	188	c Volume of a water molecule (m3)
	189	vo1 = mh2o / dble(rho_ice)
	190	c Variance of the ice and CCN distributions
	191	sigma_ice = sqrt(log(1.+nuice_sed))
	192
	193	write(,) 'Variance of ice & CCN distribs :', sigma_ice
[455]	194	write(,) 'nuice for sedimentation:', nuice_sed
[358]	195	write(,) 'Volume of a water molecule:', vo1
	196
	197	firstcall=.false.
	198	END IF
	199
	200	c-----------------------------------------------------------------------
	201	c 1. Initialization
	202	c-----------------------------------------------------------------------
[411]	203	c Initialize the tendencies
	204	pdqscloud(1:ngrid,1:nq)=0
	205	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	206	pdtcloud(1:ngrid,1:nlay)=0
	207
[358]	208	c Update the needed variables
	209	do l=1,nlay
	210	do ig=1,ngrid
	211	c Temperature
	212	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	213	c Dust mass mixing ratio
	214	zq(ig,l,igcm_dust_mass) =
	215	& pq(ig,l,igcm_dust_mass) +
	216	& pdq(ig,l,igcm_dust_mass) * ptimestep
	217	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
	218	c Dust particle number
	219	zq(ig,l,igcm_dust_number) =
[411]	220	& pq(ig,l,igcm_dust_number) +
	221	& pdq(ig,l,igcm_dust_number) * ptimestep
[358]	222	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
[411]	223	c Update rdust from last tendencies
	224	rdust(ig,l)=
	225	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
	226	& max(zq(ig,l,igcm_dust_number),0.01))
	227	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
[358]	228	c CCN mass mixing ratio
	229	zq(ig,l,igcm_ccn_mass)=
	230	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
	231	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
	232	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
	233	c CCN particle number
	234	zq(ig,l,igcm_ccn_number)=
	235	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
	236	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
	237	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
	238	c Water vapor
	239	zq(ig,l,igcm_h2o_vap)=
	240	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
	241	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
	242	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
	243	c Water ice
	244	zq(ig,l,igcm_h2o_ice)=
	245	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
	246	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
	247	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
	248	enddo
	249	enddo
	250
[411]	251
[358]	252	c------------------------------------------------------------------
	253	c Cloud microphysical scheme
	254	c------------------------------------------------------------------
	255
	256	Cste = ptimestep * 4. * pi * rho_ice
	257
	258	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	259
[411]	260	count = 0
[358]	261
	262	c Main loop over the GCM's grid
	263	DO l=1,nlay
	264	DO ig=1,ngrid
	265
	266	c Get the partial pressure of water vapor and its saturation ratio
	267	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
	268	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
[411]	269
	270	IF ((satu .ge. 1)! ) THEN ! if we have condensation
	271	& .or. ( zq(ig,l,igcm_ccn_number)
	272	& .ge. 10) ) THEN ! or sublimation
[358]	273
[411]	274
[358]	275	c Expand the dust moments into a binned distribution
[411]	276	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
	277	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
[358]	278	Rn = rdust(ig,l)
	279	Rm = Rn * exp( 3. * sigma_ice**2. )
	280	Rn = 1. / Rn
	281	Rm = 1. / Rm
	282	dev2 = 1. / ( sqrt(2.) * sigma_ice )
	283	do i = 1, nbin_cld
	284	n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
	285	& -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
	286	m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
	287	& -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
	288	enddo
	289
	290	! sumcheck = 0
	291	! do i = 1, nbin_cld
	292	! sumcheck = sumcheck + n_aer(i)
	293	! enddo
	294	! sumcheck = abs(sumcheck/No - 1)
	295	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	296	! print*, "WARNING, No sumcheck PROBLEM"
	297	! print*, "sumcheck, No",sumcheck, No
	298	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	299	! print*, "Dust binned distribution", n_aer
	300	! endif
	301	!
	302	! sumcheck = 0
	303	! do i = 1, nbin_cld
[411]	304	! sumcheck = sumcheck + m_aer(i)
[358]	305	! enddo
	306	! sumcheck = abs(sumcheck/Mo - 1)
	307	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
	308	! print*, "WARNING, Mo sumcheck PROBLEM"
[411]	309	! print*, "sumcheck, Mo",sumcheck, Mo
[358]	310	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
	311	! print*, "Dust binned distribution", m_aer
	312	! endif
	313
	314	c Get the rates of nucleation
	315	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
[411]	316
[358]	317
	318	dN = 0.
	319	dM = 0.
[455]	320	rate_out(ig,l) = 0
[358]	321	do i = 1, nbin_cld
	322	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
	323	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
	324	dN = dN + n_aer(i) * rate(i) * ptimestep
	325	dM = dM + m_aer(i) * rate(i) * ptimestep
[455]	326	rate_out(ig,l)=rate_out(ig,l)+rate(i)
[358]	327	enddo
	328
	329	! dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ),
	330	! & zq(ig,l,igcm_dust_number) )
	331	!
	332	! dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ),
	333	! & zq(ig,l,igcm_dust_mass) )
	334
	335	Mo = zq0(ig,l,igcm_h2o_ice) +
[411]	336	& zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
	337	No = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
[358]	338	!write(,) "l,cloud particles,cloud mass",l, No, Mo
	339	rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
[411]	340	& +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig)
	341	& / Mo * rho_dust
[358]	342	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
	343	rice(ig,l) =
	344	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
[411]	345	c nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
	346	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
[358]	347	seq = exp( 2.sig(zt(ig,l))mh2o /
	348	& (rho_icergpzt(ig,l)*rice(ig,l)) )
	349
	350	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
	351	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
	352
	353	up = Cste * gr * rice(ig,l) * No * seq +
	354	& zq(ig,l,igcm_h2o_vap)
	355	dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
	356
	357	Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
	358
	359	newvap = min(up/dwn,Ctot)
	360
	361	gr = gr * ( newvap/zqsat(ig,l) - seq )
	362
	363
	364	dMice = min( max(Cste * No * rice(ig,l) * gr,
	365	& -zq(ig,l,igcm_h2o_ice) ),
	366	& zq(ig,l,igcm_h2o_vap) )
	367
	368	c----------- TESTS 1D output ---------
[411]	369	satu_out(ig,l) = satu
[358]	370	Mcon_out(ig,l) = dMice
	371	newvap_out(ig,l) = newvap
	372	gr_out(ig,l) = gr
	373	dN_out(ig,l) = dN
	374	dM_out(ig,l) = dM
	375	c-------------------------------------
	376
	377	c Water ice
	378	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) +
	379	& dMice
	380
	381	c Water vapor
	382	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) -
	383	& dMice
	384
	385	c If all the ice particles sublimate, all the condensation
	386	c nuclei are released:
	387	if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
[411]	388	c for coherence
	389	dM = 0
	390	dN = 0
	391	dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
	392	dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
[358]	393	c Water ice particles
	394	zq(ig,l,igcm_h2o_ice) = 0.
	395	c Dust particles
[411]	396	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) +
	397	& zq(ig,l,igcm_ccn_mass)
	398	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) +
	399	& zq(ig,l,igcm_ccn_number)
[358]	400	c CCNs
	401	zq(ig,l,igcm_ccn_mass) = 0.
	402	zq(ig,l,igcm_ccn_number) = 0.
	403	endif
	404
[411]	405	dN = dN/ tauscaling(ig)
	406	dM = dM/ tauscaling(ig)
[358]	407	c Dust particles
[411]	408	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
	409	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
[358]	410	c CCNs
	411	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
	412	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
[411]	413
	414	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
	415	& -zq0(ig,l,igcm_dust_mass))/ptimestep
	416	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
	417	& -zq0(ig,l,igcm_dust_number))/ptimestep
	418	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
	419	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
	420	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
	421	& -zq0(ig,l,igcm_ccn_number))/ptimestep
	422	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
	423	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
	424	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
	425	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
	426
	427	count = count +1
	428	ELSE ! for coherence (rdust, rice computations etc ..)
	429	zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass)
	430	zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
	431	zq(ig,l,igcm_ccn_mass) = zq0(ig,l,igcm_ccn_mass)
	432	zq(ig,l,igcm_ccn_number) = zq0(ig,l,igcm_ccn_number)
	433	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice)
	434	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap)
[358]	435
[411]	436	! pour les sorties de test
	437	satu_out(ig,l) = satu
	438	Mcon_out(ig,l) = 0
[420]	439	newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap)
	440	gr_out(ig,l) = 0
	441	dN_out(ig,l) = 0
	442	dM_out(ig,l) = 0
[411]	443
	444	ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
	445
	446	c-----update temperature
	447	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
	448	pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
	449
	450	c----- update rice & rhocloud AFTER scavenging
	451	Mo = zq(ig,l,igcm_h2o_ice) +
	452	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
	453	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
	454	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
	455	& +zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
	456	& / Mo * rho_dust
	457	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
	458	rice(ig,l) =
	459	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
	460	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
	461
	462	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
[358]	463
[411]	464	c----- update rdust and sedimentation radius
[358]	465	rdust(ig,l)=
	466	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
	467	& max(zq(ig,l,igcm_dust_number),0.01))
	468	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
[411]	469
[358]	470	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)*3.,
	471	& rdust(ig,l) )
	472	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	473
[411]	474	ENDDO
	475	ENDDO
	476
	477
[358]	478	c------------------------------------------------------------------
	479	c------------------------------------------------------------------
	480	c------------------------------------------------------------------
	481	c------------------------------------------------------------------
[411]	482	c------------------------------------------------------------------
	483	c TESTS
[420]	484
	485	IF (output_sca) then
[411]	486
	487	print*, 'count is ',count, ' i.e. ',
	488	& count100/(nlayngrid), '% for microphys computation'
[358]	489
[420]	490	dM_col(:) = 0
	491	dN_col(:) = 0
	492	Mdust_col(:) = 0
	493	Ndust_col(:) = 0
	494	Mccn_col(:) = 0
	495	Nccn_col(:) = 0
	496	do l=1, nlay
	497	do ig=1,ngrid
	498	dM_col(ig) = dM_col(ig) +
	499	& dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
	500	dN_col(ig) = dN_col(ig) +
	501	& dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
	502	Mdust_col(ig) = Mdust_col(ig) +
	503	& zq(ig,l,igcm_dust_mass)*tauscaling(ig)
	504	& *(pplev(ig,l) - pplev(ig,l+1)) / g
	505	Ndust_col(ig) = Ndust_col(ig) +
	506	& zq(ig,l,igcm_dust_number)*tauscaling(ig)
	507	& *(pplev(ig,l) - pplev(ig,l+1)) / g
	508	Mccn_col(ig) = Mccn_col(ig) +
	509	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
	510	& *(pplev(ig,l) - pplev(ig,l+1)) / g
	511	Nccn_col(ig) = Nccn_col(ig) +
	512	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
	513	& *(pplev(ig,l) - pplev(ig,l+1)) / g
	514	enddo ! of do ig=1,ngrid
	515	enddo ! of do l=1,nlay
[358]	516
	517
[420]	518	IF (ngrid.ne.1) THEN ! 3D
[411]	519	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
	520	& satu_out)
	521	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2,
	522	& dM_col)
	523	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2,
	524	& dN_col)
	525	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2,
	526	& Ndust_col)
	527	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2,
	528	& Mdust_col)
	529	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2,
	530	& Nccn_col)
	531	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2,
	532	& Mccn_col)
	533	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
	534	& dM_out)
	535	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
	536	& dN_out)
	537	ENDIF
[358]	538
[411]	539	IF (ngrid.eq.1) THEN ! 1D
[358]	540
	541	call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
	542	& newvap_out)
	543	call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
	544	& gr_out)
[455]	545	call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
	546	& rate_out)
[358]	547	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
	548	& dM_out)
	549	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
	550	& dN_out)
	551	call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
	552	& Mcon_out)
	553	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
	554	& zqsat)
	555	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
	556	& satu_out)
[411]	557	call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
[358]	558	& rice)
[420]	559	call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
[358]	560	& rdust)
	561	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
	562	& rsedcloud)
	563	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
	564	& rhocloud)
[411]	565	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0,
	566	& dM_col)
	567	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0,
	568	& dN_col)
	569	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0,
	570	& Ndust_col)
	571	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0,
	572	& Mdust_col)
	573	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0,
	574	& Nccn_col)
	575	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0,
	576	& Mccn_col)
[358]	577	ENDIF
[420]	578	ENDIF ! endif output_sca
[358]	579	c------------------------------------------------------------------
	580	return
	581	end

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

Download in other formats:

Original Format