Context Navigation

improvedclouds.F @ 370

Last change on this file since 370 was 358, checked in by aslmd, 14 years ago

LMDZ.MARS:

A FIRST VERSION WITH SCAVENGING
TRANSPARENT TO CASUAL USER
FOR WATER CYCLE FOLKS, SEE BELOW

[TN and AS solved conflicts due to r330 >> r357]

07/11/11 == JBM

Changed watercloud.F to call two separate routines,

simpleclouds.F or improvedclouds.F, which are a simplified and
full microphysical scheme for cloud formation, respectively.
Removed the tag called "improved" in watercloud.F, and added
another tag called "microphys" which is defined in callphys.F
instead. Changed routines: callkeys, inifis, physiq, watercloud.

Reimplemented the use of typical profiles of dust particle sizes

and CCNs in simpleclouds.F. Corrected the previously used
analytical CCN profile. Moved ccn_factor to simpleclouds.F,
which won't be used in the new microphysical scheme. Changed
routines: aeropacity, physiq, simpleclouds, watercloud.

Computed rdust at the end of callsedim instead of updatereffrad,

except at firstcall. Renamed rfall into rsedcloud and moved it
in simpleclouds. Moved nuice_sed to initracer.F and added it to
tracer.h. Changed routines: callsedim, physiq, tracer.h,
watercloud, initracer, simpleclouds, updatereffrad.

Added two tracers for the CCNs. Added the different tests in

initracer.F to make sure that, for example, microphys is not
used without doubleq, etc. Corrected an inconsistency in
callsedim.F, and changed the way r0 is updated. Changes
routines: callsedim, inifis, initracer, physiq, testphys1d,
tracer.h.

Added the ability to have a spatially variable density in

newsedim (same method as that used for the radius of
sedimentation). Required the addition of one input to newsedim,
which is the size of the density variable "rho". Changed
routines: callsedim, newsedim.

Added an output to aeropacity called "tauscaling", which is a

factor to convert dust_mass and dust_number to true quantities,
based on tauvis. Removed ccn and qdust from aeropacity, which
are obsolete now.

Wrote improvedcloud.F which includes all the microphysical

scheme, and connected it to the sedimentation scheme. Added and
changed routines: callsedim, physiq, growthrate, nuclea,
improvedclouds, initracer, watercloud, microphys.h.

07/11/11 == TN

Added CCN & dust tracers updates in physiq.F

Corrected a bug that can give negative CCN numbers, due to the
use of single precision values (only 6 decimals) whereas up to 10E+10
CCN can be formed or disappears...

Corrected physical bug that causes h2o_vap or h2o_ice

to be negative in improvedclouds.F.

Corrected physical bug that causes CCN & dust tracers

to be negative in improvedclouds.F when all ice sublimates and
releases dust

Added parameter contact mteta in callphys.def

Default value is still 0.95, see inifis.F

Changed tendancy computation for dust_number in improvedclouds.F

that was not the right one. Indeed, the scavenged dust_number tracer
is computed from the dust_mass one, and its tendancy before scavenging
must be taken into account to compute its scavenging's tendancy.

File size: 22.5 KB

Rev	Line
[358]	1	subroutine improvedclouds(ngrid,nlay,ptimestep,
	2	& pplay,pt,pdt,
	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)
	39	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
	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
	113
	114	c some outputs for 1D
	115	REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
	116	REAL dN_out(ngridmx,nlayermx) ! mass variation for output
	117	REAL dM_out(ngridmx,nlayermx) ! number variation for output
	118	REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
	119	REAL gr_out(ngridmx,nlayermx) ! for 1d output
	120	REAL newvap_out(ngridmx,nlayermx) ! for 1d output
	121
	122	c------------------------------------------------------------------
	123
	124	IF (firstcall) THEN
	125
	126	c Definition of the size grid
	127	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~
	128	c rad_cld is the primary radius grid used for microphysics computation.
	129	c The grid spacing is computed assuming a constant volume ratio
	130	c between two consecutive bins; i.e. vrat_cld.
	131	c vrat_cld is determined from the boundary values of the size grid:
	132	c rmin_cld and rmax_cld.
	133	c The rb_cld array contains the boundary values of each rad_cld bin.
	134	c dr_cld is the width of each rad_cld bin.
	135
	136	c Volume ratio between two adjacent bins
	137	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
	138	vrat_cld = dexp(vrat_cld)
	139	write(,) "vrat_cld", vrat_cld
	140
	141	rb_cld(1) = rbmin_cld
	142	rad_cld(1) = rmin_cld
	143	vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3.
	144
	145	do i=1,nbin_cld-1
	146	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
	147	vol_cld(i+1) = vol_cld(i) * vrat_cld
	148	enddo
	149
	150	do i=1,nbin_cld
	151	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
	152	& rad_cld(i)
	153	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
	154	enddo
	155	rb_cld(nbin_cld+1) = rbmax_cld
	156	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
	157
	158	print*, ' '
	159	print*,'Microphysics: size bin information:'
	160	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
	161	print*,'-----------------------------------'
	162	do i=1,nbin_cld
	163	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
	164	& dr_cld(i)
	165	enddo
	166	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
	167	print*,'-----------------------------------'
	168
	169	c Contact parameter of water ice on dust ( m=cos(theta) )
	170	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	171	! mteta = 0.95
	172	write(,) 'water_param contact parameter:', mteta
	173
	174	c Volume of a water molecule (m3)
	175	vo1 = mh2o / dble(rho_ice)
	176	c Variance of the ice and CCN distributions
	177	sigma_ice = sqrt(log(1.+nuice_sed))
	178
	179	write(,) 'Variance of ice & CCN distribs :', sigma_ice
	180	write(,) 'Volume of a water molecule:', vo1
	181
	182	firstcall=.false.
	183	END IF
	184
	185	c-----------------------------------------------------------------------
	186	c 1. Initialization
	187	c-----------------------------------------------------------------------
	188	c Update the needed variables
	189
	190	c write(,) "tauscaling", tauscaling
	191
	192	!write(,) "pq ccn_mass", pq(ig,:,igcm_ccn_mass)
	193	!write(,) "pdq ccn_mass", pdq(ig,:,igcm_ccn_mass)
	194	!write(,) "pq ccn_number", pq(ig,:,igcm_ccn_number)
	195	!write(,) "pdq ccn_number", pdq(ig,:,igcm_ccn_number)
	196
	197	! print*, "improvedcloud debut pdq",
	198	! & pdq(:,:,igcm_ccn_number)*ptimestep
	199	! print*, "improvedcloud debut pq", pq(:,:,igcm_ccn_number)
	200
	201
	202	do l=1,nlay
	203	do ig=1,ngrid
	204	c Temperature
	205	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	206	c Dust mass mixing ratio
	207	c (converted to the true value using tauscaling)
	208	zq(ig,l,igcm_dust_mass) =
	209	& pq(ig,l,igcm_dust_mass) +
	210	& pdq(ig,l,igcm_dust_mass) * ptimestep
	211	zq(ig,l,igcm_dust_mass) =
	212	& zq(ig,l,igcm_dust_mass) * tauscaling(ig)
	213	zq(ig,l,igcm_dust_mass)=max(zq(ig,l,igcm_dust_mass),1.E-30)
	214	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
	215	c Dust particle number
	216	c (converted to the true value using rdust and tauscaling)
	217	! zq(ig,l,igcm_dust_number) =
	218	! & pq(ig,l,igcm_dust_number) +
	219	! & pdq(ig,l,igcm_dust_number) * ptimestep
	220	zq(ig,l,igcm_dust_number) =
	221	& (1.e0/rdust(ig,l))*3. r3n_q * zq(ig,l,igcm_dust_mass)
	222	zq(ig,l,igcm_dust_number)=
	223	& max(zq(ig,l,igcm_dust_number),1.E-30)
	224	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
	225	c CCN mass mixing ratio
	226	zq(ig,l,igcm_ccn_mass)=
	227	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
	228	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
	229	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
	230	c write(,) "pq,zq ccn_mass", pq(ig,l,igcm_ccn_mass),
	231	c & zq(ig,l,igcm_ccn_mass)
	232	c CCN particle number
	233	zq(ig,l,igcm_ccn_number)=
	234	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
	235	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
	236	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
	237	c write(,) "pq,zq ccn_number", pq(ig,l,igcm_ccn_number),
	238	c & zq(ig,l,igcm_ccn_number)
	239	c Water vapor
	240	zq(ig,l,igcm_h2o_vap)=
	241	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
	242	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
	243	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
	244	c Water ice
	245	zq(ig,l,igcm_h2o_ice)=
	246	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
	247	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
	248	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
	249	enddo
	250	enddo
	251
	252	!print*, "improvedcloud debut pq", pq(1,:,igcm_dust_number)
	253
	254	c------------------------------------------------------------------
	255	c Cloud microphysical scheme
	256	c------------------------------------------------------------------
	257
	258	Cste = ptimestep * 4. * pi * rho_ice
	259
	260	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	261
	262	c write(,) "ccn_number avant loop phy",
	263	c & zq(ig,:,igcm_ccn_number)
	264	c write(,) "ccn_mass avant loop phy",
	265	c & zq(ig,:,igcm_ccn_mass)
	266
	267	c Main loop over the GCM's grid
	268	DO l=1,nlay
	269	!ig = 1
	270	DO ig=1,ngrid
	271
	272	c Get the partial pressure of water vapor and its saturation ratio
	273	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
	274	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
	275	satu_out(ig,l) = satu
	276	!write(,) "l \| h2o_vap \| zqsat \| satu \| ph2o"
	277	!write(,) l,zq(ig,l,igcm_h2o_vap), zqsat(ig,l), satu, ph2o
	278
	279	c Expand the dust moments into a binned distribution
	280	Mo = zq(ig,l,igcm_dust_mass)
	281	No = zq(ig,l,igcm_dust_number)+ 1.e-30
	282	Rn = rdust(ig,l)
	283	Rm = Rn * exp( 3. * sigma_ice**2. )
	284	Rn = 1. / Rn
	285	Rm = 1. / Rm
	286	dev2 = 1. / ( sqrt(2.) * sigma_ice )
	287	do i = 1, nbin_cld
	288	n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
	289	& -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
	290	m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
	291	& -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
	292	enddo
	293
	294	! sumcheck = 0
	295	! do i = 1, nbin_cld
	296	! sumcheck = sumcheck + n_aer(i)
	297	! enddo
	298	! sumcheck = abs(sumcheck/No - 1)
	299	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	300	! print*, "WARNING, No sumcheck PROBLEM"
	301	! print*, "sumcheck, No",sumcheck, No
	302	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	303	! print*, "Dust binned distribution", n_aer
	304	! endif
	305	!
	306	! sumcheck = 0
	307	! do i = 1, nbin_cld
	308	! sumcheck = sumcheck + M_aer(i)
	309	! enddo
	310	! sumcheck = abs(sumcheck/Mo - 1)
	311	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
	312	! print*, "WARNING, Mo sumcheck PROBLEM"
	313	! print*, "sumcheck, No",sumcheck, No
	314	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
	315	! print*, "Dust binned distribution", m_aer
	316	! endif
	317
	318	c Get the rates of nucleation
	319	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
	320
	321	dN = 0.
	322	dM = 0.
	323	do i = 1, nbin_cld
	324	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
	325	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
	326	dN = dN + n_aer(i) * rate(i) * ptimestep
	327	dM = dM + m_aer(i) * rate(i) * ptimestep
	328	enddo
	329
	330	! dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ),
	331	! & zq(ig,l,igcm_dust_number) )
	332	!
	333	! dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ),
	334	! & zq(ig,l,igcm_dust_mass) )
	335
	336
	337	c IF (zq(ig,l,igcm_ccn_number).ge.1.e-20) THEN
	338
	339	Mo = zq0(ig,l,igcm_h2o_ice) +
	340	& zq0(ig,l,igcm_ccn_mass) + 1.e-30
	341	No = zq0(ig,l,igcm_ccn_number)
	342	!write(,) "l,cloud particles,cloud mass",l, No, Mo
	343	rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
	344	& +zq0(ig,l,igcm_ccn_mass) / Mo * rho_dust
	345	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
	346	rice(ig,l) =
	347	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
	348	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
	349	if (Mo.lt.1.e-20) rice(ig,l) = 1.e-8
	350	seq = exp( 2.sig(zt(ig,l))mh2o /
	351	& (rho_icergpzt(ig,l)*rice(ig,l)) )
	352
	353	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
	354	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
	355
	356	up = Cste * gr * rice(ig,l) * No * seq +
	357	& zq(ig,l,igcm_h2o_vap)
	358	dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
	359
	360	Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
	361
	362	newvap = min(up/dwn,Ctot)
	363
	364	gr = gr * ( newvap/zqsat(ig,l) - seq )
	365
	366
	367	dMice = min( max(Cste * No * rice(ig,l) * gr,
	368	& -zq(ig,l,igcm_h2o_ice) ),
	369	& zq(ig,l,igcm_h2o_vap) )
	370
	371	c----------- TESTS 1D output ---------
	372	if (ngrid.eq.1) then
	373	Mcon_out(ig,l) = dMice
	374	newvap_out(ig,l) = newvap
	375	gr_out(ig,l) = gr
	376	dN_out(ig,l) = dN
	377	dM_out(ig,l) = dM
	378	endif
	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
	388
	389
	390	c If all the ice particles sublimate, all the condensation
	391	c nuclei are released:
	392	if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
	393	c Water ice particles
	394	zq(ig,l,igcm_h2o_ice) = 0.
	395	c Dust particles
	396	zq(ig,l,igcm_dust_mass ) = zq0(ig,l,igcm_dust_mass) +
	397	& zq0(ig,l,igcm_ccn_mass)
	398	zq(ig,l,igcm_dust_number ) = zq0(ig,l,igcm_dust_number) +
	399	& zq0(ig,l,igcm_ccn_number)
	400	c CCNs
	401	zq(ig,l,igcm_ccn_mass) = 0.
	402	zq(ig,l,igcm_ccn_number) = 0.
	403	c for coherence
	404	dM = 0
	405	dN = 0
	406	endif
	407	c ELSE
	408	cc Initialize rhocloud and rice to avoid divisions by 0
	409	c rhocloud(ig,l) = 1.e-10
	410	c rice(ig,l) = 1.e-8
	411	c dM = 0
	412	c dN = 0
	413	c ENDIF ! of if (zq(ig,l,igcm_ccn_number).ge.1.e-20)
	414
	415	c Dust particles
	416	zq(ig,l,igcm_dust_mass ) = zq(ig,l,igcm_dust_mass ) - dM
	417	zq(ig,l,igcm_dust_number ) = zq(ig,l,igcm_dust_number ) - dN
	418	c CCNs
	419	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
	420	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
	421
	422	ENDDO
	423	ENDDO
	424
	425
	426	! print*, "improvedclouds zq0 abs.", zq0(:,:,igcm_ccn_number)
	427	! print*, "improvedclouds zq abs.", zq(:,:,igcm_ccn_number)
	428	c------------------------------------------------------------------
	429	c Convert the initial values back into relative values
	430	c (has to be done before updating rdust!)
	431	c------------------------------------------------------------------
	432
	433	do l=1, nlay
	434	do ig=1,ngrid
	435	zq0(ig,l,igcm_dust_mass) =
	436	& zq0(ig,l,igcm_dust_mass) / tauscaling(ig)
	437	zq0(ig,l,igcm_dust_number) =
	438	& (1.e0/rdust(ig,l))*3. r3n_q *
	439	& zq0(ig,l,igcm_dust_mass)
	440	enddo ! of do ig=1,ngrid
	441	enddo ! of do l=1,nlay
	442
	443	c------------------------------------------------------------------
	444	c Update the dust radius
	445	c------------------------------------------------------------------
	446
	447	DO l=1,nlay
	448	DO ig=1,ngrid
	449	rdust(ig,l)=
	450	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
	451	& max(zq(ig,l,igcm_dust_number),0.01))
	452	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
	453	ENDDO
	454	ENDDO
	455
	456	c------------------------------------------------------------------
	457	c Convert zq back into relative values (only applies to dust)
	458	c------------------------------------------------------------------
	459
	460	do l=1, nlay
	461	do ig=1,ngrid
	462	c Dust mass mixing ratio
	463	c (converted back into relative value using tauscaling)
	464	zq(ig,l,igcm_dust_mass) =
	465	& zq(ig,l,igcm_dust_mass) / tauscaling(ig)
	466	zq(ig,l,igcm_dust_mass)=max(zq(ig,l,igcm_dust_mass),1.E-30)
	467	c Dust particle number
	468	c (converted back into relative value)
	469	zq(ig,l,igcm_dust_number) =
	470	& (1.e0/rdust(ig,l))*3. r3n_q * zq(ig,l,igcm_dust_mass)
	471	zq(ig,l,igcm_dust_number)=
	472	& max(zq(ig,l,igcm_dust_number),1.E-30)
	473	enddo ! of do ig=1,ngrid
	474	enddo ! of do l=1,nlay
	475
	476	! print*, "improvedclouds zq0 rel.", zq0(1,:,igcm_ccn_number)
	477	! print*, "improvedclouds zq rel.", zq(1,:,igcm_ccn_number)
	478
	479	c------------------------------------------------------------------
	480	c Compute the sedimentation radius
	481	c------------------------------------------------------------------
	482
	483	do l=1, nlay
	484	do ig=1,ngrid
	485	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)*3.,
	486	& rdust(ig,l) )
	487	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	488	enddo ! of do ig=1,ngrid
	489	enddo ! of do l=1,nlay
	490
	491	c------------------------------------------------------------------
	492	c Force positive values
	493	c------------------------------------------------------------------
	494
	495	! do l=1, nlay
	496	! do ig=1,ngrid
	497	! zq(ig,l,igcm_ccn_mass)=
	498	! & max(zq(ig,l,igcm_ccn_mass),1e-30)
	499	! zq(ig,l,igcm_ccn_number)=
	500	! & max(zq(ig,l,igcm_ccn_number),1e-30)
	501	! zq(ig,l,igcm_h2o_vap)=
	502	! & max(zq(ig,l,igcm_h2o_vap),1e-30)
	503	! zq(ig,l,igcm_h2o_ice)=
	504	! & max(zq(ig,l,igcm_h2o_ice),1e-30)
	505	! enddo ! of do ig=1,ngrid
	506	! enddo ! of do l=1,nlay
	507
	508
	509	c------------------------------------------------------------------
	510	c Compute the final tendencies
	511	c------------------------------------------------------------------
	512
	513	c Initialize the tendencies
	514	pdqscloud(1:ngrid,1:nq)=0
	515	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	516	pdtcloud(1:ngrid,1:nlay)=0
	517
	518	c Update the tendencies
	519	do l=1, nlay
	520	do ig=1,ngrid
	521	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
	522	& -zq0(ig,l,igcm_dust_mass))/ptimestep
	523	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
	524	& -pq(ig,l,igcm_dust_number))/ptimestep
	525	& - pdq(ig,l,igcm_dust_number) !!! AJOUT TN
	526	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
	527	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
	528	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
	529	& -zq0(ig,l,igcm_ccn_number))/ptimestep
	530	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
	531	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
	532	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
	533	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
	534	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
	535	pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
	536	end do
	537	end do
	538	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
	539	! & satu_out)
	540
	541	! print, "improvedclouds pdqdt",
	542	! & pdqcloud(:,:,igcm_ccn_number)*ptimestep
	543
	544
	545	c------------------------------------------------------------------
	546	c TEST_JBM
	547	IF (ngrid.eq.1) THEN
	548	c call WRITEDIAGFI(ngrid,"tausca","tauscaling","",0,
	549	c & tauscaling)
	550	call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
	551	& newvap_out)
	552	call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
	553	& gr_out)
	554	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
	555	& dM_out)
	556	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
	557	& dN_out)
	558	call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
	559	& Mcon_out)
	560	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
	561	& zqsat)
	562	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
	563	& satu_out)
	564	call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
	565	& rice)
	566	call WRITEDIAGFI(ngrid,"rdust","rdust","m",1,
	567	& rdust)
	568	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
	569	& rsedcloud)
	570	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
	571	& rhocloud)
	572	ENDIF
	573	c------------------------------------------------------------------
	574	return
	575	end

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

Download in other formats:

Original Format