Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

callsedim2q.F @ 105

Last change on this file since 105 was 57, checked in by aslmd, 15 years ago
mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)
File size: 8.7 KB

Rev	Line
[57]	1	SUBROUTINE callsedim2q(ngrid,nlay, ptimestep,
	2	$ pplev,zlev, pt,
	3	& pq, pdqfi, pdqsed,pdqs_sed,nq)
	4	IMPLICIT NONE
	5
	6	c=======================================================================
	7	c
	8	c Sedimentation of the Martian dust
	9	c using method with mass (iq=1) and number(iq=2) mixing ratio
	10	c
	11	c F.Forget 1999
	12	c
	13	c
	14	c=======================================================================
	15
	16	c-----------------------------------------------------------------------
	17	c declarations:
	18	c -------------
	19
	20	#include "dimensions.h"
	21	#include "dimphys.h"
	22	#include "comcstfi.h"
	23	#include "tracer.h"
	24	c
	25	c arguments:
	26	c ----------
	27
	28	INTEGER ngrid,nlay
	29	REAL ptimestep ! pas de temps physique (s)
	30	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
	31	REAL pt(ngrid,nlay) ! temperature au centre des couches (K)
	32	REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries
	33
	34
	35	c Traceurs :
	36	integer nq ! nombre de traceurs
	37	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
	38	real pdqfi(ngrid,nlay,nq) ! tendance avant sedimentation (kg/kg.s-1)
	39	real pdqsed(ngrid,nlay,nq) ! tendande due a la sedimentation (kg/kg.s-1)
	40	real pdqs_sed(ngrid,nq) ! flux en surface (kg.m-2.s-1)
	41
	42	c local:
	43	c ------
	44
	45	INTEGER l,ig
	46
	47	LOGICAL firstcall
	48	SAVE firstcall
	49
	50	c Traceurs :
	51	c ~~~~~~~~
	52	INTEGER iq
	53	real zq(ngridmx,nlayermx,2)
	54	real masse(ngridmx,nlayermx) ! Layer mass (kg.m-2)
	55	real epaisseur(ngridmx,nlayermx) ! Layer thickness (m)
	56	real wq(ngridmx,nlayermx+1) ! moved tracer mass (kg.m-2)
	57	real r0(ngridmx,nlayermx) ! geometric mean radius (m)
	58
	59
	60	c Discretisation en taille:
	61	c ~~~~~~~~~~~~~~
	62	c 1) Discretisation pour representer la variation de Vitesse de chute
	63	integer nr,ir
	64	parameter (nr=7) ! nombre de rayon pour le calcul du transport
	65	real rd(nr),qr(ngridmx,nlayermx,nr)
	66	real rdi(nr+1) ! extreme and intermediate radii
	67	real Sq(ngridmx,nlayermx)
	68	real rdmin,rdmax,rdimin,rdimax
	69	data rdmin/1.e-7/
	70	data rdmax/30.e-6/
	71	data rdimin/1.e-10/
	72	data rdimax/1e-4/
	73	save rd, rdi
	74
	75	c 2) Sous discretisation pour bien integrer la loi lognormale
	76	c (calcul de q pour chaque gamme de rayon rd)
	77
	78	integer nint, iint
	79	parameter (nint=4) ! nombre de point entre chaque rayon rdi
	80	real rr(nint,nr)
	81	save rr
	82
	83	real reff(ngridmx,nlayermx,2) ! for diagnostic only
	84
	85
	86
	87
	88	DATA firstcall/.true./
	89
	90	c ** un petit test de coherence
	91	c --------------------------
	92
	93	IF (firstcall) THEN
	94	IF(ngrid.NE.ngridmx) THEN
	95	Print*,'STOP dans coefdifv'
	96	Print*,'probleme de dimensions :'
	97	Print*,'ngrid =',ngrid
	98	Print*,'ngridmx =',ngridmx
	99	STOP
	100	ENDIF
	101	firstcall=.false.
	102
	103	do ir=1,nr
	104	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
	105	c rd(ir) = rdmin + (float(ir-1)/float(nr-1))*(rdmax-rdmin)
	106	end do
	107	rdi(1)=rdimin
	108	do ir=2,nr
	109	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
	110	end do
	111	rdi(nr+1)=rdimax
	112
	113	do ir=1,nr
	114	do iint=1,nint
	115	rr(iint,ir)=
	116	& rdi(ir)(rdi(ir+1)/rdi(ir))*(float(iint-1)/float(nint-1))
	117	c write(,) rr(iint,ir)
	118	end do
	119	end do
	120	ENDIF
	121
	122
	123
	124
	125
	126	c-----------------------------------------------------------------------
	127	c 1. initialisation
	128	c -----------------
	129
	130	c On "update" la valeur de q apres les autres parametrisations
	131	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	132	do iq=1,2
	133	do l=1,nlay
	134	do ig=1,ngrid
	135	zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
	136	enddo
	137	enddo
	138	enddo
	139
	140
	141	c Calcul preliminaires de caracteristiques des couches
	142	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	143	c Masse (kg.m-2), epaisseur(m), temps de traversee (s) etc...
	144
	145	do l=1,nlay
	146	do ig=1, ngrid
	147	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	148	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
	149	end do
	150	end do
	151
	152
	153	c Calcul de la distribution de taille
	154	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	155
	156	do l=1,nlay
	157	do ig=1, ngrid
	158	r0(ig,l)=
	159	& (r3n_qzq(ig,l,1)/max(zq(ig,l,min(2,nq)),0.01))*(1./3.)
	160	r0(ig,l)=min(max(r0(ig,l),1.e-10),500.e-6)
	161	end do
	162	end do
	163
	164	DO iq=1,2 ! LOOP on Mass (iq=1) then Number(iq=2) mixing ratio
	165
	166	c Calcul du rapport de melange pour les nir tailles considerees
	167	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	168	call zerophys(ngridmx*nlayermx, Sq)
	169	do ir=1,nr
	170	do l=1,nlay
	171	do ig=1,ngrid
	172	c **************** nouvelle methode :
	173	c Calcul de q pour chaque gamme de rayon rd
	174	c Sous discretisation pour bien integrer la loi lognormale
	175	c (Methode des trapezes)
	176	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
	177	& (rr(1,ir)*(5-3iq))*
	178	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2varian**2))
	179	do iint=2,nint-1
	180	qr(ig,l,ir)=qr(ig,l,ir) +
	181	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
	182	& (rr(iint,ir)*(5-3iq))*
	183	& exp(-(log(rr(iint,ir)/r0(ig,l)))*2/(2varian**2))
	184	end do
	185	qr(ig,l,ir)=qr(ig,l,ir) +
	186	& 0.5(rr(nint,ir)-rr(nint-1,ir))
	187	& (rr(nint,ir)*(5-3iq))*
	188	& exp(-(log(rr(nint,ir)/r0(ig,l)))*2/(2varian**2))
	189
	190	c **************** Methode simple BUGGEE
	191	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
	192	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2varian**2) )
	193	c ******************************
	194	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
	195	enddo
	196	enddo
	197	enddo
	198	do ir=1,nr
	199	do l=1,nlay
	200	do ig=1,ngrid
	201	qr(ig,l,ir)= zq(ig,l,iq)* qr(ig,l,ir)/Sq(ig,l)
	202	enddo
	203	enddo
	204	c write(123,) rd(ir)1.e6, qr(1,10,ir)
	205	enddo
	206
	207	c Calcul du transport par sedimentation pour chaque traceur
	208	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	209
	210	call zerophys(ngridmx*nlayermx,zq(1,1,iq))
	211	call zerophys(ngridmx*nlayermx,pdqs_sed(1,iq))
	212
	213	c write(,)
	214	c if(iq.eq.1) write(,) '5 fois wq(6), wq(7), q(6)'
	215
	216	do ir=1,nr
	217	call newsedim(ngrid,nlay,1,ptimestep,
	218	& pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),wq)
	219
	220	c Calcul des tendances
	221	c ~~~~~~~~~~~~~~~~~~~~
	222	do ig=1,ngrid
	223	pdqs_sed(ig,iq) = pdqs_sed(ig,iq) + wq(ig,1)/ptimestep
	224	end do
	225	DO l = 1, nlay
	226	DO ig=1,ngrid
	227	zq(ig,l,iq)=zq(ig,l,iq)+qr(ig,l,ir)
	228	ENDDO
	229	ENDDO
	230	c if(iq.eq.1) write(,) 'callsed: wq(6), wq(7), q(6)',
	231	c & wq(1,6),wq(1,7),qr(1,6,ir)
	232
	233	c if(iq.eq.1)write(,)'R=', rd(ir)*1.e6
	234	c write(,)
	235	enddo
	236
	237	do l = 1, nlay
	238	do ig=1,ngrid
	239	pdqsed(ig,l,iq)=(zq(ig,l,iq)-
	240	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
	241	enddo
	242	enddo
	243
	244	c *** Diagnostique ***************
	245	c call zerophys(ngridmx*nlayermx,reff(1,1,iq))
	246	c do ir=1,nr
	247	c do l = 1, nlay
	248	c do ig=1,ngrid
	249	c if(iq.eq.1) then
	250	c reff(ig,l,iq)=reff(ig,l,iq) + qr(ig,l,ir)/rd(ir)
	251	c else
	252	c reff(ig,l,iq)=reff(ig,l,iq) + qr(ig,l,ir)*log(rd(ir))
	253	c endif
	254	c enddo
	255	c enddo
	256	c enddo
	257	c do l = 1, nlay
	258	c do ig=1,ngrid
	259	c if(iq.eq.1) then
	260	c reff(ig,l,iq)=zq(ig,l,iq)/reff(ig,l,iq)
	261	c else
	262	c reff(ig,l,iq)=exp(reff(ig,l,iq)/zq(ig,l,iq))
	263	c endif
	264	c enddo
	265	c enddo
	266	c ** FIN Diagnostique ***************
	267
	268	enddo ! end loop on iq=1,2
	269
	270
	271
	272	c ********************************************
	273	c Diagnostique
	274	c if(ngrid.eq.1) then
	275	c do l=1,nlay
	276	c do ig=1, ngrid
	277	c r0(ig,l)=
	278	c & (r3n_qzq(ig,l,1)/max(zq(ig,l,min(2,nq)),0.01))*(1./3.)
	279	c r0(ig,l)=min(max(r0(ig,l),1.e-10),500.e-6)
	280	c
	281	c write(,)'l,r0, r02 ', l, r0(ig,l),reff(ig,l,1)*0.3626
	282	c & ,reff(ig,l,2)
	283	c & ,r0(ig,l)/(reff(ig,l,1)*0.3626)
	284	c
	285	c end do
	286	c end do
	287	c CALL writeg1d(ngrid,nlay,r0,'r0','m')
	288	c end if
	289	c ********************************************
	290
	291
	292	RETURN
	293	END
	294

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

Download in other formats:

Original Format