Context Navigation

callsedim.F @ 73

Last change on this file since 73 was 38, checked in by emillour, 15 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 11.7 KB

Rev	Line
[38]	1	SUBROUTINE callsedim(ngrid,nlay, ptimestep,
	2	$ pplev,zlev, pt, rdust, rice,
	3	& pq, pdqfi, pdqsed,pdqs_sed,nq)
	4	IMPLICIT NONE
	5
	6	c=======================================================================
	7	c Sedimentation of the Martian aerosols
	8	c depending on their density and radius
	9	c
	10	c F.Forget 1999
	11	c
	12	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
	13	c technique in order to have only one call to callsedim in
	14	c physiq.F.
	15	c
	16	c=======================================================================
	17
	18	c-----------------------------------------------------------------------
	19	c declarations:
	20	c -------------
	21
	22	#include "dimensions.h"
	23	#include "dimphys.h"
	24	#include "comcstfi.h"
	25	#include "tracer.h"
	26	#include "callkeys.h"
	27
	28	c
	29	c arguments:
	30	c ----------
	31
	32	INTEGER ngrid ! number of horizontal grid points
	33	INTEGER nlay ! number of atmospheric layers
	34	REAL ptimestep ! physics time step (s)
	35	REAL pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
	36	REAL pt(ngrid,nlay) ! temperature at mid-layer (K)
	37	REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries
	38	c Aerosol radius provided by the water ice microphysical scheme:
	39	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
	40	REAL rice(ngrid,nlay) ! Ice geometric mean radius (m)
	41
	42	c Traceurs :
	43	integer nq ! number of tracers
	44	real pq(ngrid,nlay,nq) ! tracers (kg/kg)
	45	real pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
	46	real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
	47	real pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
	48
	49	c local:
	50	c ------
	51
	52	REAL CBRT
	53	EXTERNAL CBRT
	54
	55	INTEGER l,ig, iq
	56
	57	real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
	58	real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
	59	real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
	60	real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
	61	real r0(ngridmx,nlayermx) ! geometric mean doubleq radius (m)
	62	c Sedimentation radius of water ice
	63	real rfall(ngridmx,nlayermx)
	64	c Sedimentation effective variance of water ice
	65	REAL, PARAMETER :: nuice_sed = 0.1
	66
	67	c Discrete size distributions (doubleq)
	68	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	69	c 1) Parameters used to represent the changes in fall
	70	c velocity as a function of particle size;
	71	integer nr,ir
	72	parameter (nr=12) !(nr=7) ! number of bins
	73	real rd(nr),qr(ngridmx,nlayermx,nr)
	74	real rdi(nr+1) ! extreme and intermediate radii
	75	real Sq(ngridmx,nlayermx)
	76	real rdmin,rdmax,rdimin,rdimax
	77	data rdmin/1.e-8/ !/1.e-7/
	78	data rdmax/30.e-6/
	79	data rdimin/1.e-10/
	80	data rdimax/1e-4/
	81	save rd, rdi
	82
	83	c 2) Second size distribution for the log-normal integration
	84	c (the mass mixing ratio is computed for each radius)
	85
	86	integer ninter, iint
	87	parameter (ninter=4) ! nombre de point entre chaque rayon rdi
	88	real rr(ninter,nr)
	89	save rr
	90	integer radpower
	91
	92	c 3) Other local variables used in doubleq
	93
	94	real reff(ngridmx,nlayermx,2) ! for diagnostic only
	95	INTEGER idust_mass ! index of tracer containing dust mass
	96	! mix. ratio
	97	INTEGER idust_number ! index of tracer containing dust number
	98	! mix. ratio
	99	SAVE idust_mass,idust_number
	100
	101	c Firstcall:
	102
	103	LOGICAL firstcall
	104	SAVE firstcall
	105	DATA firstcall/.true./
	106
	107	c ** un petit test de coherence
	108	c --------------------------
	109
	110	IF (firstcall) THEN
	111	IF(ngrid.NE.ngridmx) THEN
	112	PRINT*,'STOP dans callsedim'
	113	PRINT*,'probleme de dimensions :'
	114	PRINT*,'ngrid =',ngrid
	115	PRINT*,'ngridmx =',ngridmx
	116	STOP
	117	ENDIF
	118
	119	c Doubleq: initialization
	120	IF (doubleq) THEN
	121	do ir=1,nr
	122	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
	123	end do
	124	rdi(1)=rdimin
	125	do ir=2,nr
	126	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
	127	end do
	128	rdi(nr+1)=rdimax
	129
	130	do ir=1,nr
	131	do iint=1,ninter
	132	rr(iint,ir)=
	133	& rdi(ir)*
	134	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
	135	c write(,) rr(iint,ir)
	136	end do
	137	end do
	138
	139	! identify tracers corresponding to mass mixing ratio and
	140	! number mixing ratio
	141	idust_mass=0 ! dummy initialization
	142	idust_number=0 ! dummy initialization
	143
	144	do iq=1,nq
	145	if (noms(iq).eq."dust_mass") then
	146	idust_mass=iq
	147	endif
	148	if (noms(iq).eq."dust_number") then
	149	idust_number=iq
	150	endif
	151	enddo
	152
	153	! check that we did find the tracers
	154	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
	155	write(,) 'callsedim: error! could not identify'
	156	write(,) ' tracers for dust mass and number mixing'
	157	write(,) ' ratio and doubleq is activated!'
	158	stop
	159	endif
	160	ENDIF !of if (doubleq)
	161
	162	IF (water) THEN
	163	write(,) "water_param nueff Sedimentation:", nuice_sed
	164	IF (activice) THEN
	165	write(,) "water_param nueff Radiative:", nuice_ref
	166	ENDIF
	167	ENDIF
	168
	169	firstcall=.false.
	170	ENDIF ! of IF (firstcall)
	171
	172	c-----------------------------------------------------------------------
	173	c 1. Initialization
	174	c -----------------
	175
	176	zqi(1:ngrid,1:nlay,1:nqmx) = 0.
	177	c Updating the mass mixing ratio with the tendencies coming
	178	c from other parameterizations:
	179	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	180
	181	do iq=1,nq
	182	do l=1,nlay
	183	do ig=1,ngrid
	184	zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
	185	enddo
	186	enddo
	187	enddo
	188
	189	c Computing the different layer properties
	190	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	191	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
	192
	193	do l=1,nlay
	194	do ig=1, ngrid
	195	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	196	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
	197	end do
	198	end do
	199
	200	c =================================================================
	201	do iq=1,nq
	202	if(radius(iq).gt.1.e-9) then ! no sedim for gaz
	203
	204	c -----------------------------------------------------------------
	205	c DOUBLEQ CASE
	206	c -----------------------------------------------------------------
	207	if (doubleq.and.
	208	& ((iq.eq.idust_mass).or.
	209	& (iq.eq.idust_number))) then
	210
	211	c Computing size distribution:
	212	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	213
	214	do l=1,nlay
	215	do ig=1, ngrid
	216	r0(ig,l)=
	217	& CBRT(r3n_q*zqi(ig,l,idust_mass)/
	218	& max(zqi(ig,l,idust_number),0.01))
	219	r0(ig,l)=min(max(r0(ig,l),1.e-10),500.e-6)
	220	end do
	221	end do
	222
	223	c Computing mass mixing ratio for each particle size
	224	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	225	IF (iq.EQ.idust_mass) then
	226	radpower = 2
	227	ELSE
	228	radpower = -1
	229	ENDIF
	230	Sq(1:ngrid,1:nlay) = 0.
	231	do ir=1,nr
	232	do l=1,nlay
	233	do ig=1,ngrid
	234	c ****************
	235	c Size distribution integration
	236	c (Trapezoid Integration Method)
	237	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
	238	& (rr(1,ir)*radpower)
	239	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2varian**2))
	240	do iint=2,ninter-1
	241	qr(ig,l,ir)=qr(ig,l,ir) +
	242	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
	243	& (rr(iint,ir)*radpower)
	244	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
	245	& (2varian*2))
	246	end do
	247	qr(ig,l,ir)=qr(ig,l,ir) +
	248	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
	249	& (rr(ninter,ir)*radpower)
	250	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
	251	& (2varian*2))
	252
	253	c **************** old method (not recommended!)
	254	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
	255	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2varian**2) )
	256	c ******************************
	257
	258	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
	259	enddo
	260	enddo
	261	enddo
	262
	263	do ir=1,nr
	264	do l=1,nlay
	265	do ig=1,ngrid
	266	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
	267	enddo
	268	enddo
	269	enddo
	270
	271	c Computing sedimentation for each tracer
	272	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	273
	274	c call zerophys(ngridmx*nlayermx,zqi(1,1,iq))
	275	zqi(1:ngrid,1:nlay,iq) = 0.
	276	c call zerophys(ngridmx,pdqs_sed(1,iq))
	277	pdqs_sed(1:ngrid,iq) = 0.
	278
	279	do ir=1,nr
	280	call newsedim(ngrid,nlay,1,ptimestep,
	281	& pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
	282	& wq,0.5)
	283
	284	c Tendencies
	285	c ~~~~~~~~~~
	286	do ig=1,ngrid
	287	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
	288	& + wq(ig,1)/ptimestep
	289	end do
	290	DO l = 1, nlay
	291	DO ig=1,ngrid
	292	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
	293	ENDDO
	294	ENDDO
	295	enddo ! of do ir=1,nr
	296	c -----------------------------------------------------------------
	297	c WATER CYCLE CASE
	298	c -----------------------------------------------------------------
	299	else if (water.and.(iq.eq.igcm_h2o_ice)) then
	300	c if (doubleq) then
	301	c do l=1,nlay
	302	c do ig=1,ngrid
	303	c rfall(ig,l)=max( rice(ig,l),rdust(ig,l) )
	304	c rfall(ig,l)=min(rfall(ig,l),1.e-4)
	305	c enddo
	306	c enddo ! of do l=1,nlay
	307	c else
	308	do l=1,nlay
	309	do ig=1,ngrid
	310	c For water cycle, a typical dust size is assumed:
	311	c r(z)=r0*exp(-z/H) with r0=0.8 micron and H=18 km.
	312	c rfall(ig,l)=max( rice(ig,l)*1.5,rdust(ig,l) )
	313	rfall(ig,l)=max( rice(ig,l)(1.+nuice_sed)*3.,
	314	& rdust(ig,l) )
	315	c modif FranckMM pour ameliorer cycle H2O: rfall= 20 microns
	316	c mars commente pour l'instant rfall(ig,l)=20.e-6
	317	rfall(ig,l)=min(rfall(ig,l),1.e-4)
	318	enddo
	319	enddo ! of do l=1,nlay
	320	c endif
	321	call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
	322	& pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi(1,1,iq),
	323	& wq,1.0)
	324	c Tendencies
	325	c ~~~~~~~~~~
	326	do ig=1,ngrid
	327	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	328	end do
	329	c -----------------------------------------------------------------
	330	c GENERAL CASE
	331	c -----------------------------------------------------------------
	332	else
	333	call newsedim(ngrid,nlay,1,ptimestep,
	334	& pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
	335	& zqi(1,1,iq),wq,1.0)
	336	c Tendencies
	337	c ~~~~~~~~~~
	338	do ig=1,ngrid
	339	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
	340	end do
	341	endif ! of if doubleq and if water
	342	c -----------------------------------------------------------------
	343
	344	DO l = 1, nlay
	345	DO ig=1,ngrid
	346	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
	347	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
	348	ENDDO
	349	ENDDO
	350
	351	endif ! of if(radius(iq).gt.1.e-9)
	352	c =================================================================
	353	enddo ! of do iq=1,nq
	354
	355	RETURN
	356	END
	357

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

Download in other formats:

Original Format