Context Navigation

callsedim.F @ 706

Last change on this file since 706 was 635, checked in by emillour, 13 years ago

Mars GCM: Update of the chemistry package, including:

93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11)
computation of photoabsorption using raytracing
improved time stepping in the photochemistry
updated parameters (cross-sections); with this new version input files

are in 'EUV/param_v5' of "datafile" directory.

transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90)
Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h
Upadted surfacearea.F
Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers).
Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense.

FGG+FL+EM

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

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

Download in other formats:

Original Format