Context Navigation

leapfrog_p.F @ 1672

Last change on this file since 1672 was 1564, checked in by emillour, 9 years ago

All GCMS:

correction wrt previous commit: inigeom is also the name of a routine

in dyn3d_common! To avoid confusion rename inigeom (in dynphy_lonlat)
inigeomphy.

cosmetic cleanup in leapfrog ('fake' calls to init_phys_lmdz,

which no longer exists, removed).
EM

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

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

Download in other formats:

Original Format