Context Navigation

leapfrog_p.F @ 1243

Last change on this file since 1243 was 1190, checked in by emillour, 11 years ago

Common dynamics:

bug fix: setting and updating halos for surface pressure (when it is modified by the physics) was missing.

File size: 63.6 KB

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

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

Download in other formats:

Original Format