Context Navigation

callsedim.F @ 414

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

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

Download in other formats:

Original Format