Context Navigation

callsedim.F @ 459

Last change on this file since 459 was 420, checked in by tnavarro, 14 years ago

24/11/11 == TN

corrected minor bug in updatereffrad.F : reffdust was not saved

ccn_factor as not correctly used in sedimentation.

It is now initialized in inifis.F, declared in tracer.h and
used in both simpleclouds.F & callsedim.F to update ice radius.

Commented diagfi outputs in aeropacity.F & improvedclouds.F for non scavenging users.

File size: 15.2 KB

Rev	Line
[38]	1	SUBROUTINE callsedim(ngrid,nlay, ptimestep,
[411]	2	& pplev,zlev, zlay, pt, rdust, rice,
[358]	3	& rsedcloud,rhocloud,
[411]	4	& pq, pdqfi, pdqsed,pdqs_sed,nq,
	5	& tau, tauscaling)
[38]	6	IMPLICIT NONE
	7
	8	c=======================================================================
	9	c Sedimentation of the Martian aerosols
	10	c depending on their density and radius
	11	c
	12	c F.Forget 1999
	13	c
	14	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
	15	c technique in order to have only one call to callsedim in
	16	c physiq.F.
	17	c
	18	c=======================================================================
	19
	20	c-----------------------------------------------------------------------
	21	c declarations:
	22	c -------------
	23
	24	#include "dimensions.h"
	25	#include "dimphys.h"
	26	#include "comcstfi.h"
	27	#include "tracer.h"
	28	#include "callkeys.h"
	29
	30	c
	31	c arguments:
	32	c ----------
	33
	34	INTEGER ngrid ! number of horizontal grid points
	35	INTEGER nlay ! number of atmospheric layers
	36	REAL ptimestep ! physics time step (s)
	37	REAL pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
	38	REAL pt(ngrid,nlay) ! temperature at mid-layer (K)
	39	REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries
	40	c Aerosol radius provided by the water ice microphysical scheme:
	41	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
	42	REAL rice(ngrid,nlay) ! Ice geometric mean radius (m)
	43
	44	c Traceurs :
	45	integer nq ! number of tracers
	46	real pq(ngrid,nlay,nq) ! tracers (kg/kg)
	47	real pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
	48	real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
	49	real pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
	50
	51	c local:
	52	c ------
	53
	54	REAL CBRT
	55	EXTERNAL CBRT
	56
	57	INTEGER l,ig, iq
	58
	59	real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
	60	real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
	61	real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
	62	real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
[358]	63	real r0(ngridmx,nlayermx) ! geometric mean radius used for
	64	! sedimentation (m)
	65	real r0dust(ngridmx,nlayermx) ! geometric mean radius used for
	66	! dust (m)
	67	real r0ccn(ngridmx,nlayermx) ! geometric mean radius used for
	68	! CCNs (m)
[38]	69	c Sedimentation radius of water ice
[358]	70	real rsedcloud(ngridmx,nlayermx)
	71	c Cloud density (kg.m-3)
	72	real rhocloud(ngridmx,nlayermx)
[411]	73
	74	c for ice radius computation
	75	REAL Mo,No
	76	REAL tau(ngrid,nlay), tauscaling(ngrid)
	77	REAL zlay(ngrid,nlay) ! altitude at the middle of the layers
	78	REAl ccntyp
[38]	79
[358]	80
[411]	81
[38]	82	c Discrete size distributions (doubleq)
	83	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	84	c 1) Parameters used to represent the changes in fall
	85	c velocity as a function of particle size;
	86	integer nr,ir
	87	parameter (nr=12) !(nr=7) ! number of bins
	88	real rd(nr),qr(ngridmx,nlayermx,nr)
	89	real rdi(nr+1) ! extreme and intermediate radii
	90	real Sq(ngridmx,nlayermx)
	91	real rdmin,rdmax,rdimin,rdimax
	92	data rdmin/1.e-8/ !/1.e-7/
	93	data rdmax/30.e-6/
	94	data rdimin/1.e-10/
	95	data rdimax/1e-4/
	96	save rd, rdi
	97
	98	c 2) Second size distribution for the log-normal integration
	99	c (the mass mixing ratio is computed for each radius)
	100
	101	integer ninter, iint
	102	parameter (ninter=4) ! nombre de point entre chaque rayon rdi
	103	real rr(ninter,nr)
	104	save rr
	105	integer radpower
[358]	106	real sigma0
[38]	107
	108	c 3) Other local variables used in doubleq
	109
	110	INTEGER idust_mass ! index of tracer containing dust mass
	111	! mix. ratio
	112	INTEGER idust_number ! index of tracer containing dust number
	113	! mix. ratio
[358]	114	INTEGER iccn_mass ! index of tracer containing CCN mass
	115	! mix. ratio
	116	INTEGER iccn_number ! index of tracer containing CCN number
	117	! mix. ratio
[38]	118	SAVE idust_mass,idust_number
[358]	119	SAVE iccn_mass,iccn_number
[38]	120
	121	c Firstcall:
	122
	123	LOGICAL firstcall
	124	SAVE firstcall
	125	DATA firstcall/.true./
	126
	127	c ** un petit test de coherence
	128	c --------------------------
	129
	130	IF (firstcall) THEN
[411]	131
[38]	132	IF(ngrid.NE.ngridmx) THEN
	133	PRINT*,'STOP dans callsedim'
	134	PRINT*,'probleme de dimensions :'
	135	PRINT*,'ngrid =',ngrid
	136	PRINT*,'ngridmx =',ngridmx
	137	STOP
	138	ENDIF
	139
	140	c Doubleq: initialization
	141	IF (doubleq) THEN
	142	do ir=1,nr
	143	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
	144	end do
	145	rdi(1)=rdimin
	146	do ir=2,nr
	147	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
	148	end do
	149	rdi(nr+1)=rdimax
	150
	151	do ir=1,nr
	152	do iint=1,ninter
	153	rr(iint,ir)=
	154	& rdi(ir)*
	155	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
	156	c write(,) rr(iint,ir)
	157	end do
	158	end do
	159
	160	! identify tracers corresponding to mass mixing ratio and
	161	! number mixing ratio
	162	idust_mass=0 ! dummy initialization
	163	idust_number=0 ! dummy initialization
	164
	165	do iq=1,nq
	166	if (noms(iq).eq."dust_mass") then
	167	idust_mass=iq
	168	endif
	169	if (noms(iq).eq."dust_number") then
	170	idust_number=iq
	171	endif
	172	enddo
	173
	174	! check that we did find the tracers
	175	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
	176	write(,) 'callsedim: error! could not identify'
	177	write(,) ' tracers for dust mass and number mixing'
	178	write(,) ' ratio and doubleq is activated!'
	179	stop
	180	endif
	181	ENDIF !of if (doubleq)
	182
[358]	183	IF (scavenging) THEN
	184	iccn_mass=0
	185	iccn_number=0
	186	do iq=1,nq
	187	if (noms(iq).eq."ccn_mass") then
	188	iccn_mass=iq
	189	endif
	190	if (noms(iq).eq."ccn_number") then
	191	iccn_number=iq
	192	endif
	193	enddo
	194	! check that we did find the tracers
	195	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
	196	write(,) 'callsedim: error! could not identify'
	197	write(,) ' tracers for ccn mass and number mixing'
	198	write(,) ' ratio and scavenging is activated!'
	199	stop
	200	endif
	201	ENDIF !of if (scavenging)
	202
[38]	203	IF (water) THEN
	204	write(,) "water_param nueff Sedimentation:", nuice_sed
	205	IF (activice) THEN
	206	write(,) "water_param nueff Radiative:", nuice_ref
	207	ENDIF
	208	ENDIF
	209
	210	firstcall=.false.
	211	ENDIF ! of IF (firstcall)
	212
	213	c-----------------------------------------------------------------------
	214	c 1. Initialization
	215	c -----------------
	216
	217	zqi(1:ngrid,1:nlay,1:nqmx) = 0.
	218	c Updating the mass mixing ratio with the tendencies coming
	219	c from other parameterizations:
	220	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	221
	222	do iq=1,nq
	223	do l=1,nlay
	224	do ig=1,ngrid
	225	zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
	226	enddo
	227	enddo
	228	enddo
	229
	230	c Computing the different layer properties
	231	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	232	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
	233
	234	do l=1,nlay
	235	do ig=1, ngrid
	236	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	237	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
	238	end do
	239	end do
	240
	241	c =================================================================
[358]	242	c Compute the geometric mean radius used for sedimentation
	243
	244	if (doubleq) then
	245	do l=1,nlay
	246	do ig=1, ngrid
	247	r0dust(ig,l) =
	248	& CBRT(r3n_q*zqi(ig,l,idust_mass)/
	249	& max(zqi(ig,l,idust_number),0.01))
	250	r0dust(ig,l)=min(max(r0dust(ig,l),1.e-10),500.e-6)
	251	end do
	252	end do
	253	endif
	254	if (scavenging) then
	255	do l=1,nlay
	256	do ig=1, ngrid
	257	r0ccn(ig,l) = rsedcloud(ig,l)/(1.+nuice_sed)**4.5
	258	end do
	259	end do
	260	endif
	261
	262	c =================================================================
[38]	263	do iq=1,nq
	264	if(radius(iq).gt.1.e-9) then ! no sedim for gaz
	265
	266	c -----------------------------------------------------------------
	267	c DOUBLEQ CASE
	268	c -----------------------------------------------------------------
[358]	269	if ((doubleq.and.
[38]	270	& ((iq.eq.idust_mass).or.
[358]	271	& (iq.eq.idust_number))).or.
	272	& (scavenging.and.
	273	& ((iq.eq.iccn_mass).or.
	274	& (iq.eq.iccn_number)))) then
[38]	275
	276	c Computing size distribution:
	277	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	278
[358]	279	if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
	280	do l=1,nlay
	281	do ig=1, ngrid
	282	r0(ig,l)=r0dust(ig,l)
	283	end do
[38]	284	end do
[358]	285	sigma0 = varian
	286	else
	287	do l=1,nlay
	288	do ig=1, ngrid
	289	r0(ig,l)=r0ccn(ig,l)
	290	end do
	291	end do
	292	sigma0 = sqrt(log(1.+nuice_sed))
	293	endif
[38]	294
	295	c Computing mass mixing ratio for each particle size
	296	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[358]	297	IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
[38]	298	radpower = 2
	299	ELSE
	300	radpower = -1
	301	ENDIF
	302	Sq(1:ngrid,1:nlay) = 0.
	303	do ir=1,nr
	304	do l=1,nlay
	305	do ig=1,ngrid
	306	c ****************
	307	c Size distribution integration
	308	c (Trapezoid Integration Method)
	309	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
	310	& (rr(1,ir)*radpower)
[358]	311	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
[38]	312	do iint=2,ninter-1
	313	qr(ig,l,ir)=qr(ig,l,ir) +
	314	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
	315	& (rr(iint,ir)*radpower)
	316	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
[358]	317	& (2sigma0*2))
[38]	318	end do
	319	qr(ig,l,ir)=qr(ig,l,ir) +
	320	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
	321	& (rr(ninter,ir)*radpower)
	322	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
[358]	323	& (2sigma0*2))
[38]	324
	325	c **************** old method (not recommended!)
	326	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
[358]	327	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
[38]	328	c ******************************
	329
	330	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
	331	enddo
	332	enddo
	333	enddo
	334
	335	do ir=1,nr
	336	do l=1,nlay
	337	do ig=1,ngrid
	338	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
	339	enddo
	340	enddo
	341	enddo
	342
	343	c Computing sedimentation for each tracer
	344	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	345
	346	zqi(1:ngrid,1:nlay,iq) = 0.
	347	pdqs_sed(1:ngrid,iq) = 0.
	348
	349	do ir=1,nr
[358]	350	IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
	351	call newsedim(ngrid,nlay,1,1,ptimestep,
	352	& pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
	353	& wq,0.5)
	354	ELSE
	355	call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
	356	& pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
	357	& wq,1.0)
	358	ENDIF
[38]	359
	360	c Tendencies
	361	c ~~~~~~~~~~
	362	do ig=1,ngrid
	363	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
	364	& + wq(ig,1)/ptimestep
	365	end do
	366	DO l = 1, nlay
	367	DO ig=1,ngrid
	368	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
	369	ENDDO
	370	ENDDO
	371	enddo ! of do ir=1,nr
	372	c -----------------------------------------------------------------
	373	c WATER CYCLE CASE
	374	c -----------------------------------------------------------------
	375	else if (water.and.(iq.eq.igcm_h2o_ice)) then
[358]	376	if (microphys) then
[411]	377	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	378	& ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
	379	& zqi(1,1,iq),wq,1.0)
[358]	380	else
[411]	381	call newsedim(ngrid,nlay,ngrid*nlay,1,
	382	& ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
	383	& zqi(1,1,iq),wq,1.0)
[358]	384	endif ! of if (microphys)
[38]	385	c Tendencies
	386	c ~~~~~~~~~~
	387	do ig=1,ngrid
	388	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	389	end do
	390	c -----------------------------------------------------------------
	391	c GENERAL CASE
	392	c -----------------------------------------------------------------
	393	else
[358]	394	call newsedim(ngrid,nlay,1,1,ptimestep,
[38]	395	& pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
	396	& zqi(1,1,iq),wq,1.0)
	397	c Tendencies
	398	c ~~~~~~~~~~
	399	do ig=1,ngrid
	400	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	401	end do
	402	endif ! of if doubleq and if water
	403	c -----------------------------------------------------------------
	404
[358]	405	c Compute the final tendency:
	406	c ---------------------------
[38]	407	DO l = 1, nlay
	408	DO ig=1,ngrid
	409	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
	410	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
	411	ENDDO
	412	ENDDO
	413
	414	endif ! of if(radius(iq).gt.1.e-9)
	415	c =================================================================
	416	enddo ! of do iq=1,nq
	417
[358]	418	c Update the dust particle size "rdust"
	419	c -------------------------------------
	420	DO l = 1, nlay
	421	DO ig=1,ngrid
	422	rdust(ig,l)=
	423	& CBRT(r3n_q*zqi(ig,l,idust_mass)/
	424	& max(zqi(ig,l,idust_number),0.01))
	425	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
	426	ENDDO
	427	ENDDO
[411]	428
	429	c Update the ice particle size "rice"
	430	c -------------------------------------
	431	IF(scavenging) THEN
	432	DO l = 1, nlay
	433	DO ig=1,ngrid
	434	Mo = zqi(ig,l,igcm_h2o_ice) +
	435	& zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
	436	No = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
	437	rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
	438	& +zqi(ig,l,iccn_mass)* tauscaling(ig)
	439	& / Mo * rho_dust
	440	rhocloud(ig,l) =
	441	& min(max(rhocloud(ig,l),rho_ice),rho_dust)
	442	rice(ig,l) =
	443	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
	444	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
	445	c print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
	446	c print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
	447
	448	ENDDO
	449	ENDDO
	450	ELSE
	451	DO l = 1, nlay
	452	DO ig=1,ngrid
	453	ccntyp =
	454	& 1.3e+8max(tau(ig,1),0.001)/0.1exp(-zlay(ig,l)/10000.)
	455	ccntyp = ccntyp /ccn_factor
	456	rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
	457	& +ccntyp(4./3.)pirdust(ig,l)*3.)
[420]	458	& /(ccntyp4./3.pi) ), rdust(ig,l))
[411]	459	ENDDO
	460	ENDDO
	461	ENDIF
[358]	462
[38]	463	RETURN
	464	END
	465

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

Download in other formats:

Original Format