Context Navigation

leapfrog_p.F @ 272

Last change on this file since 272 was 270, checked in by emillour, 14 years ago

Ehouarn: Mise a jour des dynamiques (seq et ) pour suivre
les changements et developpements de LMDZ5 terrestre
(mise a niveau avec LMDZ5 trunk, rev 1560). Ce qui ne devrais pas changer grand chose au fonctionnement de base du GCM).
Modification importante: correction du bug sur le cumul des flux de masse pour le transport des traceurs (mais dans les faits semble avoir peu d'impact).
(voir "commit_importants.log" pour les details des ajouts et modifications).

File size: 59.7 KB

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

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

Download in other formats:

Original Format