Context Navigation

leapfrog_p.F @ 123

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

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

Download in other formats:

Original Format