Context Navigation

leapfrog_p.F @ 497

Last change on this file since 497 was 497, checked in by slebonnois, 14 years ago
Ajout de scripts NCL dans UTIL. Pas encore tres developpe mais dites moi ce que vous en pensez. + correction bug pour marees grav Titan dans leapfrog
File size: 60.5 KB

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

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

Download in other formats:

Original Format