Context Navigation

improvedclouds.F @ 545

Last change on this file since 545 was 541, checked in by aslmd, 13 years ago
LMDZ.MARS. Added small optimizations in watercycle. 10% improvement.
File size: 25.7 KB

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

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

Download in other formats:

Original Format