Context Navigation

leapfrog_p.F @ 2241

Last change on this file since 2241 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 55.7 KB

Line
1	!
2	! $Id: leapfrog_p.F 2239 2015-03-23 07:27:30Z fhourdin $
3	!
4	c
5	c
6
7	SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,time_0)
8
9	use exner_hyb_m, only: exner_hyb
10	use exner_milieu_m, only: exner_milieu
11	use exner_hyb_p_m, only: exner_hyb_p
12	use exner_milieu_p_m, only: exner_milieu_p
13	USE misc_mod
14	USE parallel_lmdz
15	USE times
16	USE mod_hallo
17	USE Bands
18	USE Write_Field
19	USE Write_Field_p
20	USE vampir
21	USE timer_filtre, ONLY : print_filtre_timer
22	USE infotrac, ONLY: nqtot
23	USE guide_p_mod, ONLY : guide_main
24	USE getparam
25	USE control_mod, ONLY: nday, day_step, planet_type, offline,
26	& iconser, iphysiq, iperiod, dissip_period,
27	& iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
28	& periodav, ok_dyn_ave, output_grads_dyn,
29	& iapp_tracvl
30	IMPLICIT NONE
31
32	c ...... Version du 10/01/98 ..........
33
34	c avec coordonnees verticales hybrides
35	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
36
37	c=======================================================================
38	c
39	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
40	c -------
41	c
42	c Objet:
43	c ------
44	c
45	c GCM LMD nouvelle grille
46	c
47	c=======================================================================
48	c
49	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
50	c et possibilite d'appeler une fonction f(y) a derivee tangente
51	c hyperbolique a la place de la fonction a derivee sinusoidale.
52
53	c ... Possibilite de choisir le shema pour l'advection de
54	c q , en modifiant iadv dans traceur.def (10/02) .
55	c
56	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
57	c Pour Van-Leer iadv=10
58	c
59	c-----------------------------------------------------------------------
60	c Declarations:
61	c -------------
62
63	#include "dimensions.h"
64	#include "paramet.h"
65	#include "comconst.h"
66	#include "comdissnew.h"
67	#include "comvert.h"
68	#include "comgeom.h"
69	#include "logic.h"
70	#include "temps.h"
71	#include "ener.h"
72	#include "description.h"
73	#include "serre.h"
74	!#include "com_io_dyn.h"
75	#include "iniprint.h"
76	#include "academic.h"
77
78	REAL,INTENT(IN) :: time_0 ! not used
79
80	c dynamical variables:
81	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
82	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
83	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
84	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
85	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
86	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
87	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
88
89	REAL,SAVE :: p (ip1jmp1,llmp1 ) ! interlayer pressure
90	REAL,SAVE :: pks(ip1jmp1) ! exner at the surface
91	REAL,SAVE :: pk(ip1jmp1,llm) ! exner at mid-layer
92	REAL,SAVE :: pkf(ip1jmp1,llm) ! filtered exner at mid-layer
93	REAL,SAVE :: phi(ip1jmp1,llm) ! geopotential
94	REAL,SAVE :: w(ip1jmp1,llm) ! vertical velocity
95
96	real zqmin,zqmax
97
98	c variables dynamiques intermediaire pour le transport
99	REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
100
101	c variables dynamiques au pas -1
102	REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
103	REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
104	REAL,SAVE :: massem1(ip1jmp1,llm)
105
106	c tendances dynamiques
107	REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
108	REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
109	REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
110
111	c tendances de la dissipation
112	REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
113	REAL,SAVE :: dtetadis(ip1jmp1,llm)
114
115	c tendances physiques
116	REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
117	REAL,SAVE :: dtetafi(ip1jmp1,llm)
118	REAL,SAVE :: dpfi(ip1jmp1)
119	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
120
121	c variables pour le fichier histoire
122	REAL dtav ! intervalle de temps elementaire
123
124	REAL tppn(iim),tpps(iim),tpn,tps
125	c
126	INTEGER itau,itaufinp1,iav
127	! INTEGER iday ! jour julien
128	REAL time
129
130	REAL SSUM
131	! REAL,SAVE :: finvmaold(ip1jmp1,llm)
132
133	cym LOGICAL lafin
134	LOGICAL :: lafin
135	INTEGER ij,iq,l
136	INTEGER ik
137
138	real time_step, t_wrt, t_ops
139
140	! jD_cur: jour julien courant
141	! jH_cur: heure julienne courante
142	REAL :: jD_cur, jH_cur
143	INTEGER :: an, mois, jour
144	REAL :: secondes
145
146	logical :: physic
147	LOGICAL first,callinigrads
148
149	data callinigrads/.true./
150	character*10 string10
151
152	REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
153
154	c+jld variables test conservation energie
155	REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
156	C Tendance de la temp. potentiel d (theta)/ d t due a la
157	C tansformation d'energie cinetique en energie thermique
158	C cree par la dissipation
159	REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
160	REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
161	REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
162	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
163	CHARACTER*15 ztit
164	! INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
165	! SAVE ip_ebil_dyn
166	! DATA ip_ebil_dyn/0/
167	c-jld
168
169	character*80 dynhist_file, dynhistave_file
170	character(len=*),parameter :: modname="leapfrog"
171	character*80 abort_message
172
173
174	logical,PARAMETER :: dissip_conservative=.TRUE.
175
176	INTEGER testita
177	PARAMETER (testita = 9)
178
179	logical , parameter :: flag_verif = .false.
180
181	c declaration liees au parallelisme
182	INTEGER :: ierr
183	LOGICAL :: FirstCaldyn
184	LOGICAL :: FirstPhysic
185	INTEGER :: ijb,ije,j,i
186	type(Request) :: TestRequest
187	type(Request) :: Request_Dissip
188	type(Request) :: Request_physic
189	REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
190	REAL,SAVE :: dtetafi_tmp(iip1,llm)
191	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
192	REAL,SAVE :: dpfi_tmp(iip1)
193
194	INTEGER :: true_itau
195	INTEGER :: iapptrac
196	INTEGER :: AdjustCount
197	! INTEGER :: var_time
198	LOGICAL :: ok_start_timer=.FALSE.
199	LOGICAL, SAVE :: firstcall=.TRUE.
200
201	c$OMP MASTER
202	ItCount=0
203	c$OMP END MASTER
204	true_itau=0
205	FirstCaldyn=.TRUE.
206	FirstPhysic=.TRUE.
207	iapptrac=0
208	AdjustCount = 0
209	lafin=.false.
210
211	if (nday>=0) then
212	itaufin = nday*day_step
213	else
214	itaufin = -nday
215	endif
216
217	itaufinp1 = itaufin +1
218
219	itau = 0
220	physic=.true.
221	if (iflag_phys==0.or.iflag_phys==2) physic=.false.
222	! iday = day_ini+itau/day_step
223	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
224	! IF(time.GT.1.) THEN
225	! time = time-1.
226	! iday = iday+1
227	! ENDIF
228
229	c Allocate variables depending on dynamic variable nqtot
230	c$OMP MASTER
231	IF (firstcall) THEN
232	firstcall=.FALSE.
233	ALLOCATE(dq(ip1jmp1,llm,nqtot))
234	ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
235	ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
236	END IF
237	c$OMP END MASTER
238	c$OMP BARRIER
239
240	c-----------------------------------------------------------------------
241	c On initialise la pression et la fonction d'Exner :
242	c --------------------------------------------------
243
244	c$OMP MASTER
245	dq(:,:,:)=0.
246	CALL pression ( ip1jmp1, ap, bp, ps, p )
247	if (pressure_exner) then
248	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
249	else
250	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
251	endif
252	c$OMP END MASTER
253	c-----------------------------------------------------------------------
254	c Debut de l'integration temporelle:
255	c ----------------------------------
256	c et du parallelisme !!
257
258	1 CONTINUE ! Matsuno Forward step begins here
259
260	jD_cur = jD_ref + day_ini - day_ref + &
261	& itau/day_step
262	jH_cur = jH_ref + start_time + &
263	& mod(itau,day_step)/float(day_step)
264	if (jH_cur > 1.0 ) then
265	jD_cur = jD_cur +1.
266	jH_cur = jH_cur -1.
267	endif
268
269
270	#ifdef CPP_IOIPSL
271	if (ok_guide) then
272	!$OMP MASTER
273	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
274	!$OMP END MASTER
275	!$OMP BARRIER
276	endif
277	#endif
278
279	c
280	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
281	c CALL test_period ( ucov,vcov,teta,q,p,phis )
282	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
283	c ENDIF
284	c
285	cym CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
286	cym CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
287	cym CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
288	cym CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
289	cym CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
290
291	if (FirstCaldyn) then
292	c$OMP MASTER
293	ucovm1=ucov
294	vcovm1=vcov
295	tetam1= teta
296	massem1= masse
297	psm1= ps
298
299	! Ehouarn: finvmaold is actually not used
300	! finvmaold = masse
301	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
302	c$OMP END MASTER
303	c$OMP BARRIER
304	else
305	! Save fields obtained at previous time step as '...m1'
306	ijb=ij_begin
307	ije=ij_end
308
309	c$OMP MASTER
310	psm1 (ijb:ije) = ps (ijb:ije)
311	c$OMP END MASTER
312
313	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
314	DO l=1,llm
315	ije=ij_end
316	ucovm1 (ijb:ije,l) = ucov (ijb:ije,l)
317	tetam1 (ijb:ije,l) = teta (ijb:ije,l)
318	massem1 (ijb:ije,l) = masse (ijb:ije,l)
319	! finvmaold(ijb:ije,l)=masse(ijb:ije,l)
320
321	if (pole_sud) ije=ij_end-iip1
322	vcovm1(ijb:ije,l) = vcov (ijb:ije,l)
323
324
325	ENDDO
326	c$OMP ENDDO
327
328	! Ehouarn: finvmaold not used
329	! CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
330	! . llm, -2,2, .TRUE., 1 )
331
332	endif ! of if (FirstCaldyn)
333
334	forward = .TRUE.
335	leapf = .FALSE.
336	dt = dtvr
337
338	c ... P.Le Van .26/04/94 ....
339
340	cym CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
341	cym CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
342
343	cym ne sert a rien
344	cym call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
345
346	2 CONTINUE ! Matsuno backward or leapfrog step begins here
347
348	c$OMP MASTER
349	ItCount=ItCount+1
350	if (MOD(ItCount,1)==1) then
351	debug=.true.
352	else
353	debug=.false.
354	endif
355	c$OMP END MASTER
356	c-----------------------------------------------------------------------
357
358	c date:
359	c -----
360
361
362	c gestion des appels de la physique et des dissipations:
363	c ------------------------------------------------------
364	c
365	c ... P.Le Van ( 6/02/95 ) ....
366
367	apphys = .FALSE.
368	statcl = .FALSE.
369	conser = .FALSE.
370	apdiss = .FALSE.
371
372	IF( purmats ) THEN
373	! Purely Matsuno time stepping
374	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
375	IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
376	s apdiss = .TRUE.
377	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
378	s .and. physic ) apphys = .TRUE.
379	ELSE
380	! Leapfrog/Matsuno time stepping
381	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
382	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
383	s apdiss = .TRUE.
384	IF( MOD(itau+1,iphysiq).EQ.0.AND.physic) apphys=.TRUE.
385	END IF
386
387	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
388	! supress dissipation step
389	if (llm.eq.1) then
390	apdiss=.false.
391	endif
392
393	cym ---> Pour le moment
394	cym apphys = .FALSE.
395	statcl = .FALSE.
396	conser = .FALSE. ! ie: no output of control variables to stdout in //
397
398	if (firstCaldyn) then
399	c$OMP MASTER
400	call SetDistrib(jj_Nb_Caldyn)
401	c$OMP END MASTER
402	c$OMP BARRIER
403	firstCaldyn=.FALSE.
404	cym call InitTime
405	c$OMP MASTER
406	call Init_timer
407	c$OMP END MASTER
408	endif
409
410	c$OMP MASTER
411	IF (ok_start_timer) THEN
412	CALL InitTime
413	ok_start_timer=.FALSE.
414	ENDIF
415	c$OMP END MASTER
416
417	if (Adjust) then
418	c$OMP MASTER
419	AdjustCount=AdjustCount+1
420	if (iapptrac==iapp_tracvl .and. (forward. OR . leapf)
421	& .and. itau/iphysiq>2 .and. Adjustcount>30) then
422	AdjustCount=0
423	call allgather_timer_average
424
425	if (prt_level > 9) then
426
427	print ,'********************************'
428	print ,'*** TIMER CALDYN ****'
429	do i=0,mpi_size-1
430	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
431	& ' : temps moyen :',
432	& timer_average(jj_nb_caldyn(i),timer_caldyn,i),
433	& '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
434	enddo
435
436	print ,'********************************'
437	print ,'*** TIMER VANLEER ****'
438	do i=0,mpi_size-1
439	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
440	& ' : temps moyen :',
441	& timer_average(jj_nb_vanleer(i),timer_vanleer,i),
442	& '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
443	enddo
444
445	print ,'********************************'
446	print ,'*** TIMER DISSIP ****'
447	do i=0,mpi_size-1
448	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
449	& ' : temps moyen :',
450	& timer_average(jj_nb_dissip(i),timer_dissip,i),
451	& '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
452	enddo
453
454	if (mpi_rank==0) call WriteBands
455
456	endif
457
458	call AdjustBands_caldyn
459	if (mpi_rank==0) call WriteBands
460
461	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
462	& jj_Nb_caldyn,0,0,TestRequest)
463	call Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
464	& jj_Nb_caldyn,0,0,TestRequest)
465	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
466	& jj_Nb_caldyn,0,0,TestRequest)
467	call Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
468	& jj_Nb_caldyn,0,0,TestRequest)
469	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
470	& jj_Nb_caldyn,0,0,TestRequest)
471	call Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
472	& jj_Nb_caldyn,0,0,TestRequest)
473	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
474	& jj_Nb_caldyn,0,0,TestRequest)
475	call Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
476	& jj_Nb_caldyn,0,0,TestRequest)
477	call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
478	& jj_Nb_caldyn,0,0,TestRequest)
479	call Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
480	& jj_Nb_caldyn,0,0,TestRequest)
481	call Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
482	& jj_Nb_caldyn,0,0,TestRequest)
483	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
484	& jj_Nb_caldyn,0,0,TestRequest)
485	call Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
486	& jj_Nb_caldyn,0,0,TestRequest)
487	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
488	& jj_Nb_caldyn,0,0,TestRequest)
489	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
490	& jj_Nb_caldyn,0,0,TestRequest)
491	! call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
492	! & jj_Nb_caldyn,0,0,TestRequest)
493
494	do j=1,nqtot
495	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
496	& jj_nb_caldyn,0,0,TestRequest)
497	enddo
498
499	call SetDistrib(jj_nb_caldyn)
500	call SendRequest(TestRequest)
501	call WaitRequest(TestRequest)
502
503	call AdjustBands_dissip
504	call AdjustBands_physic
505
506	endif
507	c$OMP END MASTER
508	endif
509
510
511
512	c-----------------------------------------------------------------------
513	c calcul des tendances dynamiques:
514	c --------------------------------
515	c$OMP BARRIER
516	c$OMP MASTER
517	call VTb(VThallo)
518	c$OMP END MASTER
519
520	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
521	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
522	call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
523	call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
524	call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
525	call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
526	call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
527	call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
528
529	c do j=1,nqtot
530	c call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
531	c * TestRequest)
532	c enddo
533
534	call SendRequest(TestRequest)
535	c$OMP BARRIER
536	call WaitRequest(TestRequest)
537
538	c$OMP MASTER
539	call VTe(VThallo)
540	c$OMP END MASTER
541	c$OMP BARRIER
542
543	if (debug) then
544	!$OMP BARRIER
545	!$OMP MASTER
546	call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
547	call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
548	call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
549	call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
550	call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
551	call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
552	call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
553	call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
554	call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
555	do j=1,nqtot
556	call WriteField_p('q'//trim(int2str(j)),
557	. reshape(q(:,:,j),(/iip1,jmp1,llm/)))
558	enddo
559	!$OMP END MASTER
560	c$OMP BARRIER
561	endif
562
563
564	True_itau=True_itau+1
565
566	c$OMP MASTER
567	IF (prt_level>9) THEN
568	WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
569	ENDIF
570
571
572	call start_timer(timer_caldyn)
573
574	! compute geopotential phi()
575	CALL geopot_p ( ip1jmp1, teta , pk , pks, phis , phi )
576
577
578	call VTb(VTcaldyn)
579	c$OMP END MASTER
580	! var_time=time+iday-day_ini
581
582	c$OMP BARRIER
583	! CALL FTRACE_REGION_BEGIN("caldyn")
584	time = jD_cur + jH_cur
585	CALL caldyn_p
586	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
587	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
588
589	! CALL FTRACE_REGION_END("caldyn")
590
591	c$OMP MASTER
592	call VTe(VTcaldyn)
593	c$OMP END MASTER
594
595	cc$OMP BARRIER
596	cc$OMP MASTER
597	! call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
598	! call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
599	! call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
600	! call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
601	! call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
602	! call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
603	! call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
604	! call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
605	! call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
606	! call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
607	cc$OMP END MASTER
608
609	c-----------------------------------------------------------------------
610	c calcul des tendances advection des traceurs (dont l'humidite)
611	c -------------------------------------------------------------
612
613	IF( forward. OR . leapf ) THEN
614	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
615	CALL caladvtrac_p(q,pbaru,pbarv,
616	* p, masse, dq, teta,
617	. flxw,pk, iapptrac)
618
619	C Stokage du flux de masse pour traceurs OFF-LINE
620	IF (offline .AND. .NOT. adjust) THEN
621	CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
622	. dtvr, itau)
623	ENDIF
624
625	ENDIF ! of IF( forward. OR . leapf )
626
627	c-----------------------------------------------------------------------
628	c integrations dynamique et traceurs:
629	c ----------------------------------
630
631	c$OMP MASTER
632	call VTb(VTintegre)
633	c$OMP END MASTER
634	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
635	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
636	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
637	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
638	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
639	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
640	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
641	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
642	cc$OMP PARALLEL DEFAULT(SHARED)
643	c$OMP BARRIER
644	! CALL FTRACE_REGION_BEGIN("integrd")
645
646	CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
647	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
648	! $ finvmaold )
649
650	! CALL FTRACE_REGION_END("integrd")
651	c$OMP BARRIER
652	cc$OMP MASTER
653	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
654	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
655	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
656	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
657	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
658	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
659	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
660	c call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
661	c
662	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
663	c do j=1,nqtot
664	c call WriteField_p('q'//trim(int2str(j)),
665	c . reshape(q(:,:,j),(/iip1,jmp1,llm/)))
666	c call WriteField_p('dq'//trim(int2str(j)),
667	c . reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
668	c enddo
669	cc$OMP END MASTER
670
671
672	c$OMP MASTER
673	call VTe(VTintegre)
674	c$OMP END MASTER
675	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
676	c
677	c-----------------------------------------------------------------------
678	c calcul des tendances physiques:
679	c -------------------------------
680	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
681	c
682	IF( purmats ) THEN
683	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
684	ELSE
685	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
686	ENDIF
687
688	cc$OMP END PARALLEL
689
690	c
691	c
692	IF( apphys ) THEN
693	c
694	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
695	c
696	cc$OMP PARALLEL DEFAULT(SHARED)
697	cc$OMP+ PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
698
699	c$OMP MASTER
700	call suspend_timer(timer_caldyn)
701
702	if (prt_level >= 10) then
703	write(lunout,*)
704	& 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
705	endif
706	c$OMP END MASTER
707
708	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
709
710	c$OMP BARRIER
711	if (pressure_exner) then
712	CALL exner_hyb_p( ip1jmp1, ps, p,pks, pk, pkf )
713	else
714	CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
715	endif
716	c$OMP BARRIER
717	! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
718	! avec dyn3dmem
719	CALL geopot_p ( ip1jmp1, teta , pk , pks, phis , phi )
720
721	jD_cur = jD_ref + day_ini - day_ref
722	$ + itau/day_step
723
724	IF (planet_type .eq."generic") THEN
725	! AS: we make jD_cur to be pday
726	jD_cur = int(day_ini + itau/day_step)
727	ENDIF
728
729	jH_cur = jH_ref + start_time + &
730	& mod(itau,day_step)/float(day_step)
731	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
732	if (jH_cur > 1.0 ) then
733	jD_cur = jD_cur +1.
734	jH_cur = jH_cur -1.
735	endif
736
737	c rajout debug
738	c lafin = .true.
739
740
741	c Inbterface avec les routines de phylmd (phymars ... )
742	c -----------------------------------------------------
743
744	c+jld
745
746	c Diagnostique de conservation de l'energie : initialisation
747	IF (ip_ebil_dyn.ge.1 ) THEN
748	ztit='bil dyn'
749	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
750	IF (planet_type.eq."earth") THEN
751	#ifdef CPP_EARTH
752	CALL diagedyn(ztit,2,1,1,dtphys
753	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
754	#endif
755	ENDIF
756	ENDIF
757	c-jld
758	c$OMP BARRIER
759	c$OMP MASTER
760	call VTb(VThallo)
761	c$OMP END MASTER
762
763	call SetTag(Request_physic,800)
764
765	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
766	* jj_Nb_physic,2,2,Request_physic)
767
768	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
769	* jj_Nb_physic,2,2,Request_physic)
770
771	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
772	* jj_Nb_physic,2,2,Request_physic)
773
774	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
775	* jj_Nb_physic,1,2,Request_physic)
776
777	call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
778	* jj_Nb_physic,2,2,Request_physic)
779
780	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
781	* jj_Nb_physic,2,2,Request_physic)
782
783	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
784	* jj_Nb_physic,2,2,Request_physic)
785
786	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
787	* jj_Nb_physic,2,2,Request_physic)
788
789	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
790	* jj_Nb_physic,2,2,Request_physic)
791
792	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
793	* jj_Nb_physic,2,2,Request_physic)
794
795	c call SetDistrib(jj_nb_vanleer)
796	do j=1,nqtot
797
798	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
799	* jj_Nb_physic,2,2,Request_physic)
800	enddo
801
802	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
803	* jj_Nb_physic,2,2,Request_physic)
804
805	call SendRequest(Request_Physic)
806	c$OMP BARRIER
807	call WaitRequest(Request_Physic)
808
809	c$OMP BARRIER
810	c$OMP MASTER
811	call SetDistrib(jj_nb_Physic)
812	call VTe(VThallo)
813
814	call VTb(VTphysiq)
815	c$OMP END MASTER
816	c$OMP BARRIER
817
818	cc$OMP MASTER
819	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
820	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
821	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
822	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
823	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
824	cc$OMP END MASTER
825	cc$OMP BARRIER
826	! CALL FTRACE_REGION_BEGIN("calfis")
827	#ifdef CPP_PHYS
828	CALL calfis_p(lafin ,jD_cur, jH_cur,
829	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
830	$ du,dv,dteta,dq,
831	$ flxw, dufi,dvfi,dtetafi,dqfi,dpfi )
832	#endif
833	! CALL FTRACE_REGION_END("calfis")
834	ijb=ij_begin
835	ije=ij_end
836	if ( .not. pole_nord) then
837	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
838	DO l=1,llm
839	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
840	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
841	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
842	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
843	ENDDO
844	c$OMP END DO NOWAIT
845
846	c$OMP MASTER
847	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
848	c$OMP END MASTER
849	endif ! of if ( .not. pole_nord)
850
851	c$OMP BARRIER
852	c$OMP MASTER
853	call SetDistrib(jj_nb_Physic_bis)
854
855	call VTb(VThallo)
856	c$OMP END MASTER
857	c$OMP BARRIER
858
859	call Register_Hallo(dufi,ip1jmp1,llm,
860	* 1,0,0,1,Request_physic)
861
862	call Register_Hallo(dvfi,ip1jm,llm,
863	* 1,0,0,1,Request_physic)
864
865	call Register_Hallo(dtetafi,ip1jmp1,llm,
866	* 1,0,0,1,Request_physic)
867
868	call Register_Hallo(dpfi,ip1jmp1,1,
869	* 1,0,0,1,Request_physic)
870
871	do j=1,nqtot
872	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
873	* 1,0,0,1,Request_physic)
874	enddo
875
876	call SendRequest(Request_Physic)
877	c$OMP BARRIER
878	call WaitRequest(Request_Physic)
879
880	c$OMP BARRIER
881	c$OMP MASTER
882	call VTe(VThallo)
883
884	call SetDistrib(jj_nb_Physic)
885	c$OMP END MASTER
886	c$OMP BARRIER
887	ijb=ij_begin
888	if (.not. pole_nord) then
889
890	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
891	DO l=1,llm
892	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
893	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
894	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
895	& +dtetafi_tmp(1:iip1,l)
896	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
897	& + dqfi_tmp(1:iip1,l,:)
898	ENDDO
899	c$OMP END DO NOWAIT
900
901	c$OMP MASTER
902	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
903	c$OMP END MASTER
904
905	endif ! of if (.not. pole_nord)
906	c$OMP BARRIER
907	cc$OMP MASTER
908	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
909	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
910	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
911	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
912	cc$OMP END MASTER
913	c
914	c do j=1,nqtot
915	c call WriteField_p('dqfi'//trim(int2str(j)),
916	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
917	c enddo
918
919	c ajout des tendances physiques:
920	c ------------------------------
921	CALL addfi_p( dtphys, leapf, forward ,
922	$ ucov, vcov, teta , q ,ps ,
923	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
924	! since addfi updates ps(), also update p(), masse() and pk()
925	CALL pression_p(ip1jmp1,ap,bp,ps,p)
926	c$OMP BARRIER
927	CALL massdair_p(p,masse)
928	c$OMP BARRIER
929	if (pressure_exner) then
930	CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
931	else
932	CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
933	endif
934	c$OMP BARRIER
935
936	IF (ok_strato) THEN
937	CALL top_bound_p(vcov,ucov,teta,masse,dtphys)
938	ENDIF
939
940	c$OMP BARRIER
941	c$OMP MASTER
942	call VTe(VTphysiq)
943
944	call VTb(VThallo)
945	c$OMP END MASTER
946
947	call SetTag(Request_physic,800)
948	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
949	* jj_Nb_caldyn,Request_physic)
950
951	call Register_SwapField(vcov,vcov,ip1jm,llm,
952	* jj_Nb_caldyn,Request_physic)
953
954	call Register_SwapField(teta,teta,ip1jmp1,llm,
955	* jj_Nb_caldyn,Request_physic)
956
957	call Register_SwapField(masse,masse,ip1jmp1,llm,
958	* jj_Nb_caldyn,Request_physic)
959
960	call Register_SwapField(ps,ps,ip1jmp1,1,
961	* jj_Nb_caldyn,Request_physic)
962
963	call Register_SwapField(p,p,ip1jmp1,llmp1,
964	* jj_Nb_caldyn,Request_physic)
965
966	call Register_SwapField(pk,pk,ip1jmp1,llm,
967	* jj_Nb_caldyn,Request_physic)
968
969	call Register_SwapField(phis,phis,ip1jmp1,1,
970	* jj_Nb_caldyn,Request_physic)
971
972	call Register_SwapField(phi,phi,ip1jmp1,llm,
973	* jj_Nb_caldyn,Request_physic)
974
975	call Register_SwapField(w,w,ip1jmp1,llm,
976	* jj_Nb_caldyn,Request_physic)
977
978	do j=1,nqtot
979
980	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
981	* jj_Nb_caldyn,Request_physic)
982
983	enddo
984
985	call SendRequest(Request_Physic)
986	c$OMP BARRIER
987	call WaitRequest(Request_Physic)
988
989	c$OMP BARRIER
990	c$OMP MASTER
991	call VTe(VThallo)
992	call SetDistrib(jj_Nb_caldyn)
993	c$OMP END MASTER
994	c$OMP BARRIER
995	c
996	c Diagnostique de conservation de l'energie : difference
997	IF (ip_ebil_dyn.ge.1 ) THEN
998	ztit='bil phys'
999	CALL diagedyn(ztit,2,1,1,dtphys
1000	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
1001	ENDIF
1002
1003	cc$OMP MASTER
1004	c if (debug) then
1005	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
1006	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
1007	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
1008	c endif
1009	cc$OMP END MASTER
1010
1011
1012	c-jld
1013	c$OMP MASTER
1014	call resume_timer(timer_caldyn)
1015	if (FirstPhysic) then
1016	ok_start_timer=.TRUE.
1017	FirstPhysic=.false.
1018	endif
1019	c$OMP END MASTER
1020	ENDIF ! of IF( apphys )
1021
1022	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
1023	! Academic case : Simple friction and Newtonan relaxation
1024	! -------------------------------------------------------
1025	c$OMP MASTER
1026	if (FirstPhysic) then
1027	ok_start_timer=.TRUE.
1028	FirstPhysic=.false.
1029	endif
1030	c$OMP END MASTER
1031
1032	ijb=ij_begin
1033	ije=ij_end
1034	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1035	do l=1,llm
1036	teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
1037	& (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
1038	& (knewt_g+knewt_t(l)*clat4(ijb:ije))
1039	enddo ! of do l=1,llm
1040	!$OMP END DO
1041
1042	!$OMP MASTER
1043	if (planet_type.eq."giant") then
1044	! add an intrinsic heat flux at the base of the atmosphere
1045	teta(ijb:ije,1) = teta(ijb:ije,1)
1046	& + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
1047	endif
1048	!$OMP END MASTER
1049	!$OMP BARRIER
1050
1051	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
1052	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
1053	call SendRequest(Request_Physic)
1054	c$OMP BARRIER
1055	call WaitRequest(Request_Physic)
1056	c$OMP BARRIER
1057	call friction_p(ucov,vcov,dtvr)
1058	!$OMP BARRIER
1059
1060	! Sponge layer (if any)
1061	IF (ok_strato) THEN
1062	CALL top_bound_p(vcov,ucov,teta,masse,dtvr)
1063	!$OMP BARRIER
1064	ENDIF ! of IF (ok_strato)
1065	ENDIF ! of IF(iflag_phys.EQ.2)
1066
1067
1068	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
1069	c$OMP BARRIER
1070	if (pressure_exner) then
1071	CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
1072	else
1073	CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
1074	endif
1075	c$OMP BARRIER
1076	CALL massdair_p(p,masse)
1077	c$OMP BARRIER
1078
1079	cc$OMP END PARALLEL
1080
1081	c-----------------------------------------------------------------------
1082	c dissipation horizontale et verticale des petites echelles:
1083	c ----------------------------------------------------------
1084
1085	IF(apdiss) THEN
1086	cc$OMP PARALLEL DEFAULT(SHARED)
1087	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1088	c$OMP MASTER
1089	call suspend_timer(timer_caldyn)
1090
1091	c print*,'Entree dans la dissipation : Iteration No ',true_itau
1092	c calcul de l'energie cinetique avant dissipation
1093	c print *,'Passage dans la dissipation'
1094
1095	call VTb(VThallo)
1096	c$OMP END MASTER
1097
1098	c$OMP BARRIER
1099
1100	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1101	* jj_Nb_dissip,1,1,Request_dissip)
1102
1103	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1104	* jj_Nb_dissip,1,1,Request_dissip)
1105
1106	call Register_SwapField(teta,teta,ip1jmp1,llm,
1107	* jj_Nb_dissip,Request_dissip)
1108
1109	call Register_SwapField(p,p,ip1jmp1,llmp1,
1110	* jj_Nb_dissip,Request_dissip)
1111
1112	call Register_SwapField(pk,pk,ip1jmp1,llm,
1113	* jj_Nb_dissip,Request_dissip)
1114
1115	call SendRequest(Request_dissip)
1116	c$OMP BARRIER
1117	call WaitRequest(Request_dissip)
1118
1119	c$OMP BARRIER
1120	c$OMP MASTER
1121	call SetDistrib(jj_Nb_dissip)
1122	call VTe(VThallo)
1123	call VTb(VTdissipation)
1124	call start_timer(timer_dissip)
1125	c$OMP END MASTER
1126	c$OMP BARRIER
1127
1128	call covcont_p(llm,ucov,vcov,ucont,vcont)
1129	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1130
1131	c dissipation
1132
1133	! CALL FTRACE_REGION_BEGIN("dissip")
1134	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1135	! CALL FTRACE_REGION_END("dissip")
1136
1137	ijb=ij_begin
1138	ije=ij_end
1139	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1140	DO l=1,llm
1141	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1142	ENDDO
1143	c$OMP END DO NOWAIT
1144	if (pole_sud) ije=ije-iip1
1145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1146	DO l=1,llm
1147	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1148	ENDDO
1149	c$OMP END DO NOWAIT
1150
1151	c teta=teta+dtetadis
1152
1153
1154	c------------------------------------------------------------------------
1155	if (dissip_conservative) then
1156	C On rajoute la tendance due a la transform. Ec -> E therm. cree
1157	C lors de la dissipation
1158	c$OMP BARRIER
1159	c$OMP MASTER
1160	call suspend_timer(timer_dissip)
1161	call VTb(VThallo)
1162	c$OMP END MASTER
1163	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1164	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1165	call SendRequest(Request_Dissip)
1166	c$OMP BARRIER
1167	call WaitRequest(Request_Dissip)
1168	c$OMP MASTER
1169	call VTe(VThallo)
1170	call resume_timer(timer_dissip)
1171	c$OMP END MASTER
1172	c$OMP BARRIER
1173	call covcont_p(llm,ucov,vcov,ucont,vcont)
1174	call enercin_p(vcov,ucov,vcont,ucont,ecin)
1175
1176	ijb=ij_begin
1177	ije=ij_end
1178	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1179	do l=1,llm
1180	do ij=ijb,ije
1181	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1182	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1183	enddo
1184	enddo
1185	c$OMP END DO NOWAIT
1186	endif ! of if (dissip_conservative)
1187
1188	ijb=ij_begin
1189	ije=ij_end
1190	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1191	do l=1,llm
1192	do ij=ijb,ije
1193	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1194	enddo
1195	enddo
1196	c$OMP END DO NOWAIT
1197	c------------------------------------------------------------------------
1198
1199
1200	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1201	c ... Calcul de la valeur moyenne, unique de h aux poles .....
1202	c
1203
1204	ijb=ij_begin
1205	ije=ij_end
1206
1207	if (pole_nord) then
1208	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1209	DO l = 1, llm
1210	DO ij = 1,iim
1211	tppn(ij) = aire( ij ) * teta( ij ,l)
1212	ENDDO
1213	tpn = SSUM(iim,tppn,1)/apoln
1214
1215	DO ij = 1, iip1
1216	teta( ij ,l) = tpn
1217	ENDDO
1218	ENDDO
1219	c$OMP END DO NOWAIT
1220
1221	if (1 == 0) then
1222	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1223	!!! 2) should probably not be here anyway
1224	!!! but are kept for those who would want to revert to previous behaviour
1225	c$OMP MASTER
1226	DO ij = 1,iim
1227	tppn(ij) = aire( ij ) * ps ( ij )
1228	ENDDO
1229	tpn = SSUM(iim,tppn,1)/apoln
1230
1231	DO ij = 1, iip1
1232	ps( ij ) = tpn
1233	ENDDO
1234	c$OMP END MASTER
1235	endif ! of if (1 == 0)
1236	endif ! of of (pole_nord)
1237
1238	if (pole_sud) then
1239	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1240	DO l = 1, llm
1241	DO ij = 1,iim
1242	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1243	ENDDO
1244	tps = SSUM(iim,tpps,1)/apols
1245
1246	DO ij = 1, iip1
1247	teta(ij+ip1jm,l) = tps
1248	ENDDO
1249	ENDDO
1250	c$OMP END DO NOWAIT
1251
1252	if (1 == 0) then
1253	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1254	!!! 2) should probably not be here anyway
1255	!!! but are kept for those who would want to revert to previous behaviour
1256	c$OMP MASTER
1257	DO ij = 1,iim
1258	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1259	ENDDO
1260	tps = SSUM(iim,tpps,1)/apols
1261
1262	DO ij = 1, iip1
1263	ps(ij+ip1jm) = tps
1264	ENDDO
1265	c$OMP END MASTER
1266	endif ! of if (1 == 0)
1267	endif ! of if (pole_sud)
1268
1269
1270	c$OMP BARRIER
1271	c$OMP MASTER
1272	call VTe(VTdissipation)
1273
1274	call stop_timer(timer_dissip)
1275
1276	call VTb(VThallo)
1277	c$OMP END MASTER
1278	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1279	* jj_Nb_caldyn,Request_dissip)
1280
1281	call Register_SwapField(vcov,vcov,ip1jm,llm,
1282	* jj_Nb_caldyn,Request_dissip)
1283
1284	call Register_SwapField(teta,teta,ip1jmp1,llm,
1285	* jj_Nb_caldyn,Request_dissip)
1286
1287	call Register_SwapField(p,p,ip1jmp1,llmp1,
1288	* jj_Nb_caldyn,Request_dissip)
1289
1290	call Register_SwapField(pk,pk,ip1jmp1,llm,
1291	* jj_Nb_caldyn,Request_dissip)
1292
1293	call SendRequest(Request_dissip)
1294	c$OMP BARRIER
1295	call WaitRequest(Request_dissip)
1296
1297	c$OMP BARRIER
1298	c$OMP MASTER
1299	call SetDistrib(jj_Nb_caldyn)
1300	call VTe(VThallo)
1301	call resume_timer(timer_caldyn)
1302	c print *,'fin dissipation'
1303	c$OMP END MASTER
1304	c$OMP BARRIER
1305	END IF ! of IF(apdiss)
1306
1307	cc$OMP END PARALLEL
1308
1309	c ajout debug
1310	c IF( lafin ) then
1311	c abort_message = 'Simulation finished'
1312	c call abort_gcm(modname,abort_message,0)
1313	c ENDIF
1314
1315	c ********************************************************************
1316	c ********************************************************************
1317	c .... fin de l'integration dynamique et physique pour le pas itau ..
1318	c ********************************************************************
1319	c ********************************************************************
1320
1321	c preparation du pas d'integration suivant ......
1322	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1323	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1324	c$OMP MASTER
1325	call stop_timer(timer_caldyn)
1326	c$OMP END MASTER
1327	IF (itau==itaumax) then
1328	c$OMP MASTER
1329	call allgather_timer_average
1330
1331	if (mpi_rank==0) then
1332
1333	print ,'********************************'
1334	print ,'*** TIMER CALDYN ****'
1335	do i=0,mpi_size-1
1336	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1337	& ' : temps moyen :',
1338	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1339	enddo
1340
1341	print ,'********************************'
1342	print ,'*** TIMER VANLEER ****'
1343	do i=0,mpi_size-1
1344	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1345	& ' : temps moyen :',
1346	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1347	enddo
1348
1349	print ,'********************************'
1350	print ,'*** TIMER DISSIP ****'
1351	do i=0,mpi_size-1
1352	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1353	& ' : temps moyen :',
1354	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1355	enddo
1356
1357	print ,'********************************'
1358	print ,'*** TIMER PHYSIC ****'
1359	do i=0,mpi_size-1
1360	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1361	& ' : temps moyen :',
1362	& timer_average(jj_nb_physic(i),timer_physic,i)
1363	enddo
1364
1365	endif
1366
1367	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1368	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1369	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1370	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1371	CALL print_filtre_timer
1372	call fin_getparam
1373	call finalize_parallel
1374	c$OMP END MASTER
1375	c$OMP BARRIER
1376	RETURN
1377	ENDIF ! of IF (itau==itaumax)
1378
1379	IF ( .NOT.purmats ) THEN
1380	c ........................................................
1381	c .............. schema matsuno + leapfrog ..............
1382	c ........................................................
1383
1384	IF(forward. OR. leapf) THEN
1385	itau= itau + 1
1386	! iday= day_ini+itau/day_step
1387	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1388	! IF(time.GT.1.) THEN
1389	! time = time-1.
1390	! iday = iday+1
1391	! ENDIF
1392	ENDIF
1393
1394
1395	IF( itau. EQ. itaufinp1 ) then
1396
1397	if (flag_verif) then
1398	write(79,*) 'ucov',ucov
1399	write(80,*) 'vcov',vcov
1400	write(81,*) 'teta',teta
1401	write(82,*) 'ps',ps
1402	write(83,*) 'q',q
1403	WRITE(85,*) 'q1 = ',q(:,:,1)
1404	WRITE(86,*) 'q3 = ',q(:,:,3)
1405	endif
1406
1407
1408	c$OMP MASTER
1409	call fin_getparam
1410	c$OMP END MASTER
1411	#ifdef INCA
1412	call finalize_inca
1413	#endif
1414	c$OMP MASTER
1415	call finalize_parallel
1416	c$OMP END MASTER
1417	abort_message = 'Simulation finished'
1418	call abort_gcm(modname,abort_message,0)
1419	RETURN
1420	ENDIF
1421	c-----------------------------------------------------------------------
1422	c ecriture du fichier histoire moyenne:
1423	c -------------------------------------
1424
1425	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1426	c$OMP BARRIER
1427	IF(itau.EQ.itaufin) THEN
1428	iav=1
1429	ELSE
1430	iav=0
1431	ENDIF
1432	#ifdef CPP_IOIPSL
1433	IF (ok_dynzon) THEN
1434	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1435	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1436	ENDIF !ok_dynzon
1437	#endif
1438	IF (ok_dyn_ave) THEN
1439	!$OMP MASTER
1440	#ifdef CPP_IOIPSL
1441	! Ehouarn: Gather fields and make master send to output
1442	call Gather_Field(vcov,ip1jm,llm,0)
1443	call Gather_Field(ucov,ip1jmp1,llm,0)
1444	call Gather_Field(teta,ip1jmp1,llm,0)
1445	call Gather_Field(pk,ip1jmp1,llm,0)
1446	call Gather_Field(phi,ip1jmp1,llm,0)
1447	do iq=1,nqtot
1448	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1449	enddo
1450	call Gather_Field(masse,ip1jmp1,llm,0)
1451	call Gather_Field(ps,ip1jmp1,1,0)
1452	call Gather_Field(phis,ip1jmp1,1,0)
1453	if (mpi_rank==0) then
1454	CALL writedynav(itau,vcov,
1455	& ucov,teta,pk,phi,q,masse,ps,phis)
1456	endif
1457	#endif
1458	!$OMP END MASTER
1459	ENDIF ! of IF (ok_dyn_ave)
1460	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1461
1462	c-----------------------------------------------------------------------
1463	c ecriture de la bande histoire:
1464	c ------------------------------
1465
1466	IF( MOD(itau,iecri).EQ.0) THEN
1467	! Ehouarn: output only during LF or Backward Matsuno
1468	if (leapf.or.(.not.leapf.and.(.not.forward))) then
1469	c$OMP BARRIER
1470	c$OMP MASTER
1471	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1472
1473	cym unat=0.
1474
1475	ijb=ij_begin
1476	ije=ij_end
1477
1478	if (pole_nord) then
1479	ijb=ij_begin+iip1
1480	unat(1:iip1,:)=0.
1481	endif
1482
1483	if (pole_sud) then
1484	ije=ij_end-iip1
1485	unat(ij_end-iip1+1:ij_end,:)=0.
1486	endif
1487
1488	do l=1,llm
1489	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1490	enddo
1491
1492	ijb=ij_begin
1493	ije=ij_end
1494	if (pole_sud) ije=ij_end-iip1
1495
1496	do l=1,llm
1497	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1498	enddo
1499
1500	#ifdef CPP_IOIPSL
1501	if (ok_dyn_ins) then
1502	! Ehouarn: Gather fields and make master write to output
1503	call Gather_Field(vcov,ip1jm,llm,0)
1504	call Gather_Field(ucov,ip1jmp1,llm,0)
1505	call Gather_Field(teta,ip1jmp1,llm,0)
1506	call Gather_Field(phi,ip1jmp1,llm,0)
1507	do iq=1,nqtot
1508	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1509	enddo
1510	call Gather_Field(masse,ip1jmp1,llm,0)
1511	call Gather_Field(ps,ip1jmp1,1,0)
1512	call Gather_Field(phis,ip1jmp1,1,0)
1513	if (mpi_rank==0) then
1514	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1515	endif
1516	! CALL writehist_p(histid,histvid, itau,vcov,
1517	! & ucov,teta,phi,q,masse,ps,phis)
1518	! or use writefield_p
1519	! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1520	! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1521	! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1522	! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1523	endif ! of if (ok_dyn_ins)
1524	#endif
1525	! For some Grads outputs of fields
1526	if (output_grads_dyn) then
1527	! Ehouarn: hope this works the way I think it does:
1528	call Gather_Field(unat,ip1jmp1,llm,0)
1529	call Gather_Field(vnat,ip1jm,llm,0)
1530	call Gather_Field(teta,ip1jmp1,llm,0)
1531	call Gather_Field(ps,ip1jmp1,1,0)
1532	do iq=1,nqtot
1533	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1534	enddo
1535	if (mpi_rank==0) then
1536	#include "write_grads_dyn.h"
1537	endif
1538	endif ! of if (output_grads_dyn)
1539	c$OMP END MASTER
1540	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1541	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1542
1543	IF(itau.EQ.itaufin) THEN
1544
1545	c$OMP BARRIER
1546	c$OMP MASTER
1547
1548	! if (planet_type.eq."earth") then
1549	! Write an Earth-format restart file
1550	CALL dynredem1_p("restart.nc",0.0,
1551	& vcov,ucov,teta,q,masse,ps)
1552	! endif ! of if (planet_type.eq."earth")
1553
1554	! CLOSE(99)
1555	c$OMP END MASTER
1556	ENDIF ! of IF (itau.EQ.itaufin)
1557
1558	c-----------------------------------------------------------------------
1559	c gestion de l'integration temporelle:
1560	c ------------------------------------
1561
1562	IF( MOD(itau,iperiod).EQ.0 ) THEN
1563	GO TO 1
1564	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1565
1566	IF( forward ) THEN
1567	c fin du pas forward et debut du pas backward
1568
1569	forward = .FALSE.
1570	leapf = .FALSE.
1571	GO TO 2
1572
1573	ELSE
1574	c fin du pas backward et debut du premier pas leapfrog
1575
1576	leapf = .TRUE.
1577	dt = 2.*dtvr
1578	GO TO 2
1579	END IF
1580	ELSE
1581
1582	c ...... pas leapfrog .....
1583
1584	leapf = .TRUE.
1585	dt = 2.*dtvr
1586	GO TO 2
1587	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1588	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1589
1590
1591	ELSE ! of IF (.not.purmats)
1592
1593	c ........................................................
1594	c .............. schema matsuno ...............
1595	c ........................................................
1596	IF( forward ) THEN
1597
1598	itau = itau + 1
1599	! iday = day_ini+itau/day_step
1600	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1601	!
1602	! IF(time.GT.1.) THEN
1603	! time = time-1.
1604	! iday = iday+1
1605	! ENDIF
1606
1607	forward = .FALSE.
1608	IF( itau. EQ. itaufinp1 ) then
1609	c$OMP MASTER
1610	call fin_getparam
1611	call finalize_parallel
1612	c$OMP END MASTER
1613	abort_message = 'Simulation finished'
1614	call abort_gcm(modname,abort_message,0)
1615	RETURN
1616	ENDIF
1617	GO TO 2
1618
1619	ELSE ! of IF(forward) i.e. backward step
1620
1621	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1622	IF(itau.EQ.itaufin) THEN
1623	iav=1
1624	ELSE
1625	iav=0
1626	ENDIF
1627	#ifdef CPP_IOIPSL
1628	IF (ok_dynzon) THEN
1629	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1630	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1631	END IF !ok_dynzon
1632	#endif
1633	IF (ok_dyn_ave) THEN
1634	!$OMP MASTER
1635	#ifdef CPP_IOIPSL
1636	! Ehouarn: Gather fields and make master send to output
1637	call Gather_Field(vcov,ip1jm,llm,0)
1638	call Gather_Field(ucov,ip1jmp1,llm,0)
1639	call Gather_Field(teta,ip1jmp1,llm,0)
1640	call Gather_Field(pk,ip1jmp1,llm,0)
1641	call Gather_Field(phi,ip1jmp1,llm,0)
1642	do iq=1,nqtot
1643	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1644	enddo
1645	call Gather_Field(masse,ip1jmp1,llm,0)
1646	call Gather_Field(ps,ip1jmp1,1,0)
1647	call Gather_Field(phis,ip1jmp1,1,0)
1648	if (mpi_rank==0) then
1649	CALL writedynav(itau,vcov,
1650	& ucov,teta,pk,phi,q,masse,ps,phis)
1651	endif
1652	#endif
1653	!$OMP END MASTER
1654	ENDIF ! of IF (ok_dyn_ave)
1655
1656	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1657
1658
1659	IF(MOD(itau,iecri ).EQ.0) THEN
1660	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1661	c$OMP BARRIER
1662	c$OMP MASTER
1663	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1664
1665	cym unat=0.
1666	ijb=ij_begin
1667	ije=ij_end
1668
1669	if (pole_nord) then
1670	ijb=ij_begin+iip1
1671	unat(1:iip1,:)=0.
1672	endif
1673
1674	if (pole_sud) then
1675	ije=ij_end-iip1
1676	unat(ij_end-iip1+1:ij_end,:)=0.
1677	endif
1678
1679	do l=1,llm
1680	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1681	enddo
1682
1683	ijb=ij_begin
1684	ije=ij_end
1685	if (pole_sud) ije=ij_end-iip1
1686
1687	do l=1,llm
1688	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1689	enddo
1690
1691	#ifdef CPP_IOIPSL
1692	if (ok_dyn_ins) then
1693	! Ehouarn: Gather fields and make master send to output
1694	call Gather_Field(vcov,ip1jm,llm,0)
1695	call Gather_Field(ucov,ip1jmp1,llm,0)
1696	call Gather_Field(teta,ip1jmp1,llm,0)
1697	call Gather_Field(phi,ip1jmp1,llm,0)
1698	do iq=1,nqtot
1699	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1700	enddo
1701	call Gather_Field(masse,ip1jmp1,llm,0)
1702	call Gather_Field(ps,ip1jmp1,1,0)
1703	call Gather_Field(phis,ip1jmp1,1,0)
1704	if (mpi_rank==0) then
1705	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1706	endif
1707	! CALL writehist_p(histid, histvid, itau,vcov ,
1708	! & ucov,teta,phi,q,masse,ps,phis)
1709	endif ! of if (ok_dyn_ins)
1710	#endif
1711	! For some Grads output (but does it work?)
1712	if (output_grads_dyn) then
1713	call Gather_Field(unat,ip1jmp1,llm,0)
1714	call Gather_Field(vnat,ip1jm,llm,0)
1715	call Gather_Field(teta,ip1jmp1,llm,0)
1716	call Gather_Field(ps,ip1jmp1,1,0)
1717	do iq=1,nqtot
1718	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1719	enddo
1720	c
1721	if (mpi_rank==0) then
1722	#include "write_grads_dyn.h"
1723	endif
1724	endif ! of if (output_grads_dyn)
1725
1726	c$OMP END MASTER
1727	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1728
1729	IF(itau.EQ.itaufin) THEN
1730	! if (planet_type.eq."earth") then
1731	c$OMP MASTER
1732	CALL dynredem1_p("restart.nc",0.0,
1733	. vcov,ucov,teta,q,masse,ps)
1734	c$OMP END MASTER
1735	! endif ! of if (planet_type.eq."earth")
1736	ENDIF ! of IF(itau.EQ.itaufin)
1737
1738	forward = .TRUE.
1739	GO TO 1
1740
1741	ENDIF ! of IF (forward)
1742
1743	END IF ! of IF(.not.purmats)
1744	c$OMP MASTER
1745	call fin_getparam
1746	call finalize_parallel
1747	c$OMP END MASTER
1748	RETURN
1749	END

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

Download in other formats:

Original Format