Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

callsedim2q.F @ 1242

Last change on this file since 1242 was 57, checked in by aslmd, 14 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

Line
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