Context Navigation

leapfrog_p.F @ 2125

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

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

Download in other formats:

Original Format