Context Navigation

improvedclouds.F @ 524

Last change on this file since 524 was 520, checked in by tnavarro, 13 years ago

10/02/12 == TN

Major update on watercycle: a smaller integration timestep is now used

in watercloud.F, sedimentation of clouds is done in watercloud instead of
callsedim.F

Temperature-dependant contact parameter in nuclea.F
No dust lifting if CO2 ice in vdif.c
Ice integrated column opacity is written in diagfi from physiq.F, instead

of aeropacity.F. Mandatory if iradia is not 1.

New definition of permanent ice in surfini.F and possibility to have an ice

cap in it in 1d.

Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging

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

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

Download in other formats:

Original Format