Context Navigation

advx.F @ 802

Last change on this file since 802 was 644, checked in by Laurent Fairhead, 20 years ago
Synchronisation avec tous les diagnostiques de Ionela IM Inclusion du slab ocean IM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.1 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE advx(limit,dtx,pbaru,sm,s0,
5	$ sx,sy,sz,lati,latf)
6	IMPLICIT NONE
7
8	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9	C C
10	C first-order moments (FOM) advection of tracer in X direction C
11	C C
12	C Source : Pascal Simon (Meteo,CNRM) C
13	C Adaptation : A.Armengaud (LGGE) juin 94 C
14	C C
15	C limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz C
16	C sont des arguments d'entree pour le s-pg... C
17	C C
18	C sm,s0,sx,sy,sz C
19	C sont les arguments de sortie pour le s-pg C
20	C C
21	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
22	C
23	C parametres principaux du modele
24	C
25	#include "dimensions.h"
26	#include "paramet.h"
27	#include "comconst.h"
28	#include "comvert.h"
29
30	C Arguments :
31	C -----------
32	C dtx : frequence fictive d'appel du transport
33	C pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
34
35	INTEGER ntra
36	PARAMETER (ntra = 1)
37
38	C ATTENTION partout ou on trouve ntra, insertion de boucle
39	C possible dans l'avenir.
40
41	REAL dtx
42	REAL pbaru ( iip1,jjp1,llm )
43
44	C moments: SM total mass in each grid box
45	C S0 mass of tracer in each grid box
46	C Si 1rst order moment in i direction
47	C
48	REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
49	REAL sx(iip1,jjp1,llm,ntra)
50	$ ,sy(iip1,jjp1,llm,ntra)
51	REAL sz(iip1,jjp1,llm,ntra)
52
53	C Local :
54	C -------
55
56	C mass fluxes across the boundaries (UGRI,VGRI,WGRI)
57	C mass fluxes in kg
58	C declaration :
59
60	REAL UGRI(iip1,jjp1,llm)
61
62	C Rem : VGRI et WGRI ne sont pas utilises dans
63	C cette subroutine ( advection en x uniquement )
64	C
65	C Ti are the moments for the current latitude and level
66	C
67	REAL TM(iim)
68	REAL T0(iim,ntra),TX(iim,ntra)
69	REAL TY(iim,ntra),TZ(iim,ntra)
70	REAL TEMPTM ! just a temporary variable
71	C
72	C the moments F are similarly defined and used as temporary
73	C storage for portions of the grid boxes in transit
74	C
75	REAL FM(iim)
76	REAL F0(iim,ntra),FX(iim,ntra)
77	REAL FY(iim,ntra),FZ(iim,ntra)
78	C
79	C work arrays
80	C
81	REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
82	C
83	REAL SMNEW(iim),UEXT(iim)
84	C
85	REAL sqi,sqf
86
87	LOGICAL LIMIT
88	INTEGER NUM(jjp1),LONK,NUMK
89	INTEGER lon,lati,latf,niv
90	INTEGER i,i2,i3,j,jv,l,k,itrac
91
92	lon = iim
93	niv = llm
94
95	C * Test de passage d'arguments ****
96
97
98	C -------------------------------------
99	DO 300 j = 1,jjp1
100	NUM(j) = 1
101	300 CONTINUE
102	sqi = 0.
103	sqf = 0.
104
105	DO l = 1,llm
106	DO j = 1,jjp1
107	DO i = 1,iim
108	cIM 240305 sqi = sqi + S0(i,j,l,9)
109	sqi = sqi + S0(i,j,l,ntra)
110	ENDDO
111	ENDDO
112	ENDDO
113	PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
114	PRINT*,'sqi=',sqi
115
116
117	C Interface : adaptation nouveau modele
118	C -------------------------------------
119	C
120	C ---------------------------------------------------------
121	C Conversion des flux de masses en kg/s
122	C pbaru est en N/s d'ou :
123	C ugri est en kg/s
124
125	DO 500 l = 1,llm
126	DO 500 j = 1,jjm+1
127	DO 500 i = 1,iip1
128	C ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
129	ugri (i,j,llm+1-l) = pbaru (i,j,l)
130	500 CONTINUE
131
132
133	C ---------------------------------------------------------
134	C ---------------------------------------------------------
135	C ---------------------------------------------------------
136
137	C start here
138	C
139	C boucle principale sur les niveaux et les latitudes
140	C
141	DO 1 L=1,NIV
142	DO 1 K=lati,latf
143	C
144	C initialisation
145	C
146	C program assumes periodic boundaries in X
147	C
148	DO 10 I=2,LON
149	SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
150	10 CONTINUE
151	SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
152	C
153	C modifications for extended polar zones
154	C
155	NUMK=NUM(K)
156	LONK=LON/NUMK
157	C
158	IF(NUMK.GT.1) THEN
159	C
160	DO 111 I=1,LON
161	TM(I)=0.
162	111 CONTINUE
163	DO 112 JV=1,NTRA
164	DO 1120 I=1,LON
165	T0(I,JV)=0.
166	TX(I,JV)=0.
167	TY(I,JV)=0.
168	TZ(I,JV)=0.
169	1120 CONTINUE
170	112 CONTINUE
171	C
172	DO 11 I2=1,NUMK
173	C
174	DO 113 I=1,LONK
175	I3=(I-1)*NUMK+I2
176	TM(I)=TM(I)+SM(I3,K,L)
177	ALF(I)=SM(I3,K,L)/TM(I)
178	ALF1(I)=1.-ALF(I)
179	113 CONTINUE
180	C
181	DO JV=1,NTRA
182	DO I=1,LONK
183	I3=(I-1)*NUMK+I2
184	TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)
185	$ *S0(I3,K,L,JV)
186	T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
187	TX(I,JV)=ALF(I) *sx(I3,K,L,JV)+
188	$ ALF1(I)TX(I,JV) +3.TEMPTM
189	TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
190	TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
191	ENDDO
192	ENDDO
193	C
194	11 CONTINUE
195	C
196	ELSE
197	C
198	DO 115 I=1,LON
199	TM(I)=SM(I,K,L)
200	115 CONTINUE
201	DO 116 JV=1,NTRA
202	DO 1160 I=1,LON
203	T0(I,JV)=S0(I,K,L,JV)
204	TX(I,JV)=sx(I,K,L,JV)
205	TY(I,JV)=sy(I,K,L,JV)
206	TZ(I,JV)=sz(I,K,L,JV)
207	1160 CONTINUE
208	116 CONTINUE
209	C
210	ENDIF
211	C
212	DO 117 I=1,LONK
213	UEXT(I)=UGRI(I*NUMK,K,L)
214	117 CONTINUE
215	C
216	C place limits on appropriate moments before transport
217	C (if flux-limiting is to be applied)
218	C
219	IF(.NOT.LIMIT) GO TO 13
220	C
221	DO 12 JV=1,NTRA
222	DO 120 I=1,LONK
223	TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
224	120 CONTINUE
225	12 CONTINUE
226	C
227	13 CONTINUE
228	C
229	C calculate flux and moments between adjacent boxes
230	C 1- create temporary moments/masses for partial boxes in transit
231	C 2- reajusts moments remaining in the box
232	C
233	C flux from IP to I if U(I).lt.0
234	C
235	DO 140 I=1,LONK-1
236	IF(UEXT(I).LT.0.) THEN
237	FM(I)=-UEXT(I)*DTX
238	ALF(I)=FM(I)/TM(I+1)
239	TM(I+1)=TM(I+1)-FM(I)
240	ENDIF
241	140 CONTINUE
242	C
243	I=LONK
244	IF(UEXT(I).LT.0.) THEN
245	FM(I)=-UEXT(I)*DTX
246	ALF(I)=FM(I)/TM(1)
247	TM(1)=TM(1)-FM(I)
248	ENDIF
249	C
250	C flux from I to IP if U(I).gt.0
251	C
252	DO 141 I=1,LONK
253	IF(UEXT(I).GE.0.) THEN
254	FM(I)=UEXT(I)*DTX
255	ALF(I)=FM(I)/TM(I)
256	TM(I)=TM(I)-FM(I)
257	ENDIF
258	141 CONTINUE
259	C
260	DO 142 I=1,LONK
261	ALFQ(I)=ALF(I)*ALF(I)
262	ALF1(I)=1.-ALF(I)
263	ALF1Q(I)=ALF1(I)*ALF1(I)
264	142 CONTINUE
265	C
266	DO 150 JV=1,NTRA
267	DO 1500 I=1,LONK-1
268	C
269	IF(UEXT(I).LT.0.) THEN
270	C
271	F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
272	FX(I,JV)=ALFQ(I)*TX(I+1,JV)
273	FY(I,JV)=ALF (I)*TY(I+1,JV)
274	FZ(I,JV)=ALF (I)*TZ(I+1,JV)
275	C
276	T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
277	TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
278	TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
279	TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
280	C
281	ENDIF
282	C
283	1500 CONTINUE
284	150 CONTINUE
285	C
286	I=LONK
287	IF(UEXT(I).LT.0.) THEN
288	C
289	DO 151 JV=1,NTRA
290	C
291	F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
292	FX (I,JV)=ALFQ(I)*TX(1,JV)
293	FY (I,JV)=ALF (I)*TY(1,JV)
294	FZ (I,JV)=ALF (I)*TZ(1,JV)
295	C
296	T0(1,JV)=T0(1,JV)-F0(I,JV)
297	TX(1,JV)=ALF1Q(I)*TX(1,JV)
298	TY(1,JV)=TY(1,JV)-FY(I,JV)
299	TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
300	C
301	151 CONTINUE
302	C
303	ENDIF
304	C
305	DO 152 JV=1,NTRA
306	DO 1520 I=1,LONK
307	C
308	IF(UEXT(I).GE.0.) THEN
309	C
310	F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
311	FX(I,JV)=ALFQ(I)*TX(I,JV)
312	FY(I,JV)=ALF (I)*TY(I,JV)
313	FZ(I,JV)=ALF (I)*TZ(I,JV)
314	C
315	T0(I,JV)=T0(I,JV)-F0(I,JV)
316	TX(I,JV)=ALF1Q(I)*TX(I,JV)
317	TY(I,JV)=TY(I,JV)-FY(I,JV)
318	TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
319	C
320	ENDIF
321	C
322	1520 CONTINUE
323	152 CONTINUE
324	C
325	C puts the temporary moments Fi into appropriate neighboring boxes
326	C
327	DO 160 I=1,LONK
328	IF(UEXT(I).LT.0.) THEN
329	TM(I)=TM(I)+FM(I)
330	ALF(I)=FM(I)/TM(I)
331	ENDIF
332	160 CONTINUE
333	C
334	DO 161 I=1,LONK-1
335	IF(UEXT(I).GE.0.) THEN
336	TM(I+1)=TM(I+1)+FM(I)
337	ALF(I)=FM(I)/TM(I+1)
338	ENDIF
339	161 CONTINUE
340	C
341	I=LONK
342	IF(UEXT(I).GE.0.) THEN
343	TM(1)=TM(1)+FM(I)
344	ALF(I)=FM(I)/TM(1)
345	ENDIF
346	C
347	DO 162 I=1,LONK
348	ALF1(I)=1.-ALF(I)
349	162 CONTINUE
350	C
351	DO 170 JV=1,NTRA
352	DO 1700 I=1,LONK
353	C
354	IF(UEXT(I).LT.0.) THEN
355	C
356	TEMPTM=-ALF(I)T0(I,JV)+ALF1(I)F0(I,JV)
357	T0(I,JV)=T0(I,JV)+F0(I,JV)
358	TX(I,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(I,JV)+3.*TEMPTM
359	TY(I,JV)=TY(I,JV)+FY(I,JV)
360	TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
361	C
362	ENDIF
363	C
364	1700 CONTINUE
365	170 CONTINUE
366	C
367	DO 171 JV=1,NTRA
368	DO 1710 I=1,LONK-1
369	C
370	IF(UEXT(I).GE.0.) THEN
371	C
372	TEMPTM=ALF(I)T0(I+1,JV)-ALF1(I)F0(I,JV)
373	T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
374	TX(I+1,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(I+1,JV)+3.*TEMPTM
375	TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
376	TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
377	C
378	ENDIF
379	C
380	1710 CONTINUE
381	171 CONTINUE
382	C
383	I=LONK
384	IF(UEXT(I).GE.0.) THEN
385	DO 172 JV=1,NTRA
386	TEMPTM=ALF(I)T0(1,JV)-ALF1(I)F0(I,JV)
387	T0(1,JV)=T0(1,JV)+F0(I,JV)
388	TX(1,JV)=ALF(I)FX(I,JV)+ALF1(I)TX(1,JV)+3.*TEMPTM
389	TY(1,JV)=TY(1,JV)+FY(I,JV)
390	TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
391	172 CONTINUE
392	ENDIF
393	C
394	C retour aux mailles d'origine (passage des Tij aux Sij)
395	C
396	IF(NUMK.GT.1) THEN
397	C
398	DO 180 I2=1,NUMK
399	C
400	DO 180 I=1,LONK
401	C
402	I3=I2+(I-1)*NUMK
403	SM(I3,K,L)=SMNEW(I3)
404	ALF(I)=SMNEW(I3)/TM(I)
405	TM(I)=TM(I)-SMNEW(I3)
406	C
407	ALFQ(I)=ALF(I)*ALF(I)
408	ALF1(I)=1.-ALF(I)
409	ALF1Q(I)=ALF1(I)*ALF1(I)
410	C
411	180 CONTINUE
412	C
413	DO JV=1,NTRA
414	DO I=1,LONK
415	C
416	I3=I2+(I-1)*NUMK
417	S0(I3,K,L,JV)=ALF (I)
418	$ * (T0(I,JV)-ALF1(I)*TX(I,JV))
419	sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
420	sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
421	sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
422	C
423	C reajusts moments remaining in the box
424	C
425	T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
426	TX(I,JV)=ALF1Q(I)*TX(I,JV)
427	TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
428	TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
429	ENDDO
430	ENDDO
431	C
432	C
433	ELSE
434	C
435	DO 190 I=1,LON
436	SM(I,K,L)=TM(I)
437	190 CONTINUE
438	DO 191 JV=1,NTRA
439	DO 1910 I=1,LON
440	S0(I,K,L,JV)=T0(I,JV)
441	sx(I,K,L,JV)=TX(I,JV)
442	sy(I,K,L,JV)=TY(I,JV)
443	sz(I,K,L,JV)=TZ(I,JV)
444	1910 CONTINUE
445	191 CONTINUE
446	C
447	ENDIF
448	C
449	1 CONTINUE
450	C
451	C ----------- AA Test en fin de ADVX ------ Controle des S*
452	c OK
453	c DO 9998 l = 1, llm
454	c DO 9998 j = 1, jjp1
455	c DO 9998 i = 1, iip1
456	c IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
457	c PRINT*, '-------------------'
458	c PRINT*, 'En fin de ADVX'
459	c PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
460	c PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
461	c print*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
462	c print*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
463	c print*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
464	c WRITE (,) 'On arrete !! - pbl en fin de ADVX1'
465	cc STOP
466	c ENDIF
467	c 9998 CONTINUE
468	c
469	C ---------- bouclage cyclique
470	DO itrac=1,ntra
471	DO l = 1,llm
472	DO j = lati,latf
473	SM(iip1,j,l) = SM(1,j,l)
474	S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
475	sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
476	sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
477	sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
478	END DO
479	END DO
480	ENDDO
481
482	c ----------- qqtite totale de traceur dans tte l'atmosphere
483	DO l = 1, llm
484	DO j = 1, jjp1
485	DO i = 1, iim
486	cIM 240405 sqf = sqf + S0(i,j,l,9)
487	sqf = sqf + S0(i,j,l,ntra)
488	END DO
489	END DO
490	END DO
491	c
492	PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
493	PRINT*,'sqf=',sqf
494	c-------------
495
496	RETURN
497	END
498	C_________________________________________________________________
499	C_________________________________________________________________

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

Download in other formats:

Original Format