Context Navigation

leapfrog_p.F @ 907

Last change on this file since 907 was 847, checked in by aslmd, 13 years ago
LMDZ.COMMON. Corrected bugs with using variable cp in parallel. Corrected bugs related to problems when no tracers are used. Updated makelmdz_fcm with the latest version and added a few lines necessary for generic physics. Added a build_gcm script in csh. Updated AMD64_CICLAD compilation settings. Uploaded arch files to make the model work with ifort on ciclad.
File size: 61.7 KB

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

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

Download in other formats:

Original Format