Context Navigation

leapfrog_p.F @ 127

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

Ehouarn: suite de l'implementation de la discretisation verticale,
avec quelques mises a jour pour concorder avec la version terrestre.
-> Finalement, on met un flag "disvert_type" pour fixer la discretisation

disvert_type==1 (defaut si planet_type=="earth") pour cas terrestre
disvert_type==2 (defaut si planet_type!="earth") pour cas planeto (z2sig.def)

-> au passage, pour rester en phase avec modele terrestre on renomme

disvert_terre en disvert (le disvert "alternatif" demeure 'disvert_noterre')

File size: 59.7 KB

Line
1	!
2	! $Id: leapfrog_p.F 1446 2010-10-22 09:27:25Z emillour $
3	!
4	c
5	c
6
7	SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,
8	& time_0)
9
10	USE misc_mod
11	USE parallel
12	USE times
13	USE mod_hallo
14	USE Bands
15	USE Write_Field
16	USE Write_Field_p
17	USE vampir
18	USE timer_filtre, ONLY : print_filtre_timer
19	USE infotrac
20	USE guide_p_mod, ONLY : guide_main
21	USE getparam
22	USE control_mod
23
24	IMPLICIT NONE
25
26	c ...... Version du 10/01/98 ..........
27
28	c avec coordonnees verticales hybrides
29	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
30
31	c=======================================================================
32	c
33	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
34	c -------
35	c
36	c Objet:
37	c ------
38	c
39	c GCM LMD nouvelle grille
40	c
41	c=======================================================================
42	c
43	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
44	c et possibilite d'appeler une fonction f(y) a derivee tangente
45	c hyperbolique a la place de la fonction a derivee sinusoidale.
46
47	c ... Possibilite de choisir le shema pour l'advection de
48	c q , en modifiant iadv dans traceur.def (10/02) .
49	c
50	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
51	c Pour Van-Leer iadv=10
52	c
53	c-----------------------------------------------------------------------
54	c Declarations:
55	c -------------
56
57	#include "dimensions.h"
58	#include "paramet.h"
59	#include "comconst.h"
60	#include "comdissnew.h"
61	#include "comvert.h"
62	#include "comgeom.h"
63	#include "logic.h"
64	#include "temps.h"
65	#include "ener.h"
66	#include "description.h"
67	#include "serre.h"
68	!#include "com_io_dyn.h"
69	#include "iniprint.h"
70	#include "academic.h"
71
72	real zqmin,zqmax
73	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 (disvert_type==1) then
274	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
275	else ! we assume that we are in the disvert_type==2 case
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	c$OMP MASTER
603	IF (prt_level>9) THEN
604	WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
605	ENDIF
606
607
608	call start_timer(timer_caldyn)
609
610	! ADAPTATION GCM POUR CP(T)
611	! CALL geopot_p ( ip1jmp1, teta , pk , pks, phis , phi )
612	CALL geopot_p ( ip1jmp1, tsurpk , pk , pks, phis , phi )
613
614	call VTb(VTcaldyn)
615	c$OMP END MASTER
616	! var_time=time+iday-day_ini
617
618	c$OMP BARRIER
619	! CALL FTRACE_REGION_BEGIN("caldyn")
620	time = jD_cur + jH_cur
621	! ADAPTATION GCM POUR CP(T)
622	! CALL caldyn_p
623	! $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
624	! $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
625	CALL caldyn_p
626	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,tsurpk,phis,
627	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
628
629	! CALL FTRACE_REGION_END("caldyn")
630
631	c$OMP MASTER
632	call VTe(VTcaldyn)
633	c$OMP END MASTER
634
635	cc$OMP BARRIER
636	cc$OMP MASTER
637	! call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
638	! call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
639	! call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
640	! call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
641	! call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
642	! call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
643	! call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
644	! call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
645	! call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
646	! call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
647	cc$OMP END MASTER
648
649	c-----------------------------------------------------------------------
650	c calcul des tendances advection des traceurs (dont l'humidite)
651	c -------------------------------------------------------------
652
653	IF( forward. OR . leapf ) THEN
654	cc$OMP PARALLEL DEFAULT(SHARED)
655	c
656	CALL caladvtrac_p(q,pbaru,pbarv,
657	* p, masse, dq, teta,
658	. flxw,pk, iapptrac)
659
660	C Stokage du flux de masse pour traceurs OFF-LINE
661	IF (offline .AND. .NOT. adjust) THEN
662	CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
663	. dtvr, itau)
664	ENDIF
665
666	ENDIF ! of IF( forward. OR . leapf )
667	cc$OMP END PARALLEL
668
669	c-----------------------------------------------------------------------
670	c integrations dynamique et traceurs:
671	c ----------------------------------
672
673	c$OMP MASTER
674	call VTb(VTintegre)
675	c$OMP END MASTER
676	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
677	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
678	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
679	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
680	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
681	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
682	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
683	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
684	cc$OMP PARALLEL DEFAULT(SHARED)
685	c$OMP BARRIER
686	! CALL FTRACE_REGION_BEGIN("integrd")
687
688	CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
689	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
690	$ finvmaold )
691
692	! CALL FTRACE_REGION_END("integrd")
693	c$OMP BARRIER
694	cc$OMP MASTER
695	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
696	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
697	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
698	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
699	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
700	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
701	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
702	c call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
703	c
704	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
705	c do j=1,nqtot
706	c call WriteField_p('q'//trim(int2str(j)),
707	c . reshape(q(:,:,j),(/iip1,jmp1,llm/)))
708	c call WriteField_p('dq'//trim(int2str(j)),
709	c . reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
710	c enddo
711	cc$OMP END MASTER
712
713
714	c$OMP MASTER
715	call VTe(VTintegre)
716	c$OMP END MASTER
717	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
718	c
719	c-----------------------------------------------------------------------
720	c calcul des tendances physiques:
721	c -------------------------------
722	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
723	c
724	IF( purmats ) THEN
725	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
726	ELSE
727	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
728	ENDIF
729
730	cc$OMP END PARALLEL
731
732	c
733	c
734	IF( apphys ) THEN
735	c
736	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
737	c
738	cc$OMP PARALLEL DEFAULT(SHARED)
739	cc$OMP+ PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
740
741	c$OMP MASTER
742	call suspend_timer(timer_caldyn)
743
744	if (prt_level >= 10) then
745	write(lunout,*)
746	& 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
747	endif
748	c$OMP END MASTER
749
750	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
751
752	c$OMP BARRIER
753	if (disvert_type==1) then
754	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
755	else ! we assume that we are in the disvert_type==2 case
756	CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
757	endif
758	c$OMP BARRIER
759	jD_cur = jD_ref + day_ini - day_ref
760	$ + int (itau * dtvr / daysec)
761	jH_cur = jH_ref + &
762	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
763	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
764
765	c rajout debug
766	c lafin = .true.
767
768
769	c Interface avec les routines de phylmd (phymars ... )
770	c -----------------------------------------------------
771
772	c+jld
773
774	c Diagnostique de conservation de l'energie : initialisation
775	IF (ip_ebil_dyn.ge.1 ) THEN
776	ztit='bil dyn'
777	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
778	IF (planet_type.eq."earth") THEN
779	CALL diagedyn(ztit,2,1,1,dtphys
780	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
781	ENDIF
782	ENDIF
783	c-jld
784	c$OMP BARRIER
785	c$OMP MASTER
786	call VTb(VThallo)
787	c$OMP END MASTER
788
789	call SetTag(Request_physic,800)
790
791	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
792	* jj_Nb_physic,2,2,Request_physic)
793
794	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
795	* jj_Nb_physic,2,2,Request_physic)
796
797	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
798	* jj_Nb_physic,2,2,Request_physic)
799
800	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
801	* jj_Nb_physic,1,2,Request_physic)
802
803	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
804	* jj_Nb_physic,2,2,Request_physic)
805
806	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
807	* jj_Nb_physic,2,2,Request_physic)
808
809	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
810	* jj_Nb_physic,2,2,Request_physic)
811
812	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
813	* jj_Nb_physic,2,2,Request_physic)
814
815	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
816	* jj_Nb_physic,2,2,Request_physic)
817
818	c call SetDistrib(jj_nb_vanleer)
819	do j=1,nqtot
820
821	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
822	* jj_Nb_physic,2,2,Request_physic)
823	enddo
824
825	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
826	* jj_Nb_physic,2,2,Request_physic)
827
828	call SendRequest(Request_Physic)
829	c$OMP BARRIER
830	call WaitRequest(Request_Physic)
831
832	c$OMP BARRIER
833	c$OMP MASTER
834	call SetDistrib(jj_nb_Physic)
835	call VTe(VThallo)
836
837	call VTb(VTphysiq)
838	c$OMP END MASTER
839	c$OMP BARRIER
840
841	cc$OMP MASTER
842	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
843	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
844	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
845	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
846	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
847	cc$OMP END MASTER
848	cc$OMP BARRIER
849	! CALL FTRACE_REGION_BEGIN("calfis")
850	CALL calfis_p(lafin ,jD_cur, jH_cur,
851	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
852	$ du,dv,dteta,dq,
853	$ flxw,
854	$ dufi,dvfi,dtetafi,dqfi,dpfi )
855	! CALL FTRACE_REGION_END("calfis")
856	ijb=ij_begin
857	ije=ij_end
858	if ( .not. pole_nord) then
859	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
860	DO l=1,llm
861	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
862	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
863	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
864	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
865	ENDDO
866	c$OMP END DO NOWAIT
867
868	c$OMP MASTER
869	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
870	c$OMP END MASTER
871	endif ! of if ( .not. pole_nord)
872
873	c$OMP BARRIER
874	c$OMP MASTER
875	call SetDistrib(jj_nb_Physic_bis)
876
877	call VTb(VThallo)
878	c$OMP END MASTER
879	c$OMP BARRIER
880
881	call Register_Hallo(dufi,ip1jmp1,llm,
882	* 1,0,0,1,Request_physic)
883
884	call Register_Hallo(dvfi,ip1jm,llm,
885	* 1,0,0,1,Request_physic)
886
887	call Register_Hallo(dtetafi,ip1jmp1,llm,
888	* 1,0,0,1,Request_physic)
889
890	call Register_Hallo(dpfi,ip1jmp1,1,
891	* 1,0,0,1,Request_physic)
892
893	do j=1,nqtot
894	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
895	* 1,0,0,1,Request_physic)
896	enddo
897
898	call SendRequest(Request_Physic)
899	c$OMP BARRIER
900	call WaitRequest(Request_Physic)
901
902	c$OMP BARRIER
903	c$OMP MASTER
904	call VTe(VThallo)
905
906	call SetDistrib(jj_nb_Physic)
907	c$OMP END MASTER
908	c$OMP BARRIER
909	ijb=ij_begin
910	if (.not. pole_nord) then
911
912	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
913	DO l=1,llm
914	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
915	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
916	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
917	& +dtetafi_tmp(1:iip1,l)
918	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
919	& + dqfi_tmp(1:iip1,l,:)
920	ENDDO
921	c$OMP END DO NOWAIT
922
923	c$OMP MASTER
924	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
925	c$OMP END MASTER
926
927	endif ! of if (.not. pole_nord)
928	c$OMP BARRIER
929	cc$OMP MASTER
930	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
931	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
932	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
933	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
934	cc$OMP END MASTER
935	c
936	c do j=1,nqtot
937	c call WriteField_p('dqfi'//trim(int2str(j)),
938	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
939	c enddo
940
941	c ajout des tendances physiques:
942	c ------------------------------
943	CALL addfi_p( dtphys, leapf, forward ,
944	$ ucov, vcov, teta , q ,ps ,
945	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
946
947	c Couche superieure :
948	c -------------------
949	IF (ok_strato) THEN
950	CALL top_bound_p( vcov,ucov,teta,phi,masse,
951	$ dutop,dvtop,dtetatop)
952	CALL addfi_p( dtphys, leapf, forward ,
953	$ ucov, vcov, teta , q ,ps ,
954	$ dutop, dvtop, dtetatop , dqtop ,dptop )
955
956	ENDIF
957
958	c$OMP BARRIER
959	c$OMP MASTER
960	call VTe(VTphysiq)
961
962	call VTb(VThallo)
963	c$OMP END MASTER
964
965	call SetTag(Request_physic,800)
966	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
967	* jj_Nb_caldyn,Request_physic)
968
969	call Register_SwapField(vcov,vcov,ip1jm,llm,
970	* jj_Nb_caldyn,Request_physic)
971
972	call Register_SwapField(teta,teta,ip1jmp1,llm,
973	* jj_Nb_caldyn,Request_physic)
974
975	call Register_SwapField(masse,masse,ip1jmp1,llm,
976	* jj_Nb_caldyn,Request_physic)
977
978	call Register_SwapField(p,p,ip1jmp1,llmp1,
979	* jj_Nb_caldyn,Request_physic)
980
981	call Register_SwapField(pk,pk,ip1jmp1,llm,
982	* jj_Nb_caldyn,Request_physic)
983
984	call Register_SwapField(phis,phis,ip1jmp1,1,
985	* jj_Nb_caldyn,Request_physic)
986
987	call Register_SwapField(phi,phi,ip1jmp1,llm,
988	* jj_Nb_caldyn,Request_physic)
989
990	call Register_SwapField(w,w,ip1jmp1,llm,
991	* jj_Nb_caldyn,Request_physic)
992
993	do j=1,nqtot
994
995	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
996	* jj_Nb_caldyn,Request_physic)
997
998	enddo
999
1000	call SendRequest(Request_Physic)
1001	c$OMP BARRIER
1002	call WaitRequest(Request_Physic)
1003
1004	c$OMP BARRIER
1005	c$OMP MASTER
1006	call VTe(VThallo)
1007	call SetDistrib(jj_Nb_caldyn)
1008	c$OMP END MASTER
1009	c$OMP BARRIER
1010	c
1011	c Diagnostique de conservation de l'energie : difference
1012	IF (ip_ebil_dyn.ge.1 ) THEN
1013	ztit='bil phys'
1014	CALL diagedyn(ztit,2,1,1,dtphys
1015	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
1016	ENDIF
1017
1018	cc$OMP MASTER
1019	c if (debug) then
1020	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
1021	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
1022	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
1023	c endif
1024	cc$OMP END MASTER
1025
1026
1027	c-jld
1028	c$OMP MASTER
1029	call resume_timer(timer_caldyn)
1030	if (FirstPhysic) then
1031	ok_start_timer=.TRUE.
1032	FirstPhysic=.false.
1033	endif
1034	c$OMP END MASTER
1035	ENDIF ! of IF( apphys )
1036
1037	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
1038	! Academic case : Simple friction and Newtonan relaxation
1039	! -------------------------------------------------------
1040	ijb=ij_begin
1041	ije=ij_end
1042	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1043	do l=1,llm
1044	teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
1045	& (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
1046	& (knewt_g+knewt_t(l)*clat4(ijb:ije))
1047	enddo ! of do l=1,llm
1048	!$OMP END DO
1049
1050	if (planet_type.eq."giant") then
1051	! Intrinsic heat flux
1052	! Aymeric -- for giant planets
1053	if (ihf .gt. 1.e-6) then
1054	!print , '* INTRINSIC HEAT FLUX **', ihf
1055	teta(ijb:ije,1) = teta(ijb:ije,1)
1056	& + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
1057	!print , '*** d teta '
1058	!print , dtvr aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
1059	endif
1060	endif
1061
1062	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
1063	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
1064	call SendRequest(Request_Physic)
1065	c$OMP BARRIER
1066	call WaitRequest(Request_Physic)
1067	c$OMP BARRIER
1068	call friction_p(ucov,vcov,dtvr)
1069	!$OMP BARRIER
1070
1071	! Sponge layer (if any)
1072	IF (ok_strato) THEN
1073	! set dutop,dvtop,... to zero
1074	ijb=ij_begin
1075	ije=ij_end
1076	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1077	do l=1,llm
1078	dutop(ijb:ije,l)=0
1079	dtetatop(ijb:ije,l)=0
1080	dqtop(ijb:ije,l,1:nqtot)=0
1081	enddo
1082	!$OMP END DO
1083	!$OMP SINGLE
1084	dptop(ijb:ije)=0
1085	!$OMP END SINGLE
1086	ijb=ij_begin
1087	ije=ij_end
1088	if (pole_sud) ije=ije-iip1
1089	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1090	do l=1,llm
1091	dvtop(ijb:ije,l)=0
1092	enddo
1093	!$OMP END DO
1094
1095	CALL top_bound_p(vcov,ucov,teta,phi,masse,
1096	$ dutop,dvtop,dtetatop)
1097	CALL addfi_p( dtvr, leapf, forward ,
1098	$ ucov, vcov, teta , q ,ps ,
1099	$ dutop, dvtop, dtetatop , dqtop ,dptop )
1100	!$OMP BARRIER
1101	ENDIF ! of IF (ok_strato)
1102	ENDIF ! of IF(iflag_phys.EQ.2)
1103
1104
1105	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
1106	c$OMP BARRIER
1107	if (disvert_type==1) then
1108	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
1109	else ! we assume that we are in the disvert_type==2 case
1110	CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
1111	endif
1112	c$OMP BARRIER
1113
1114	cc$OMP END PARALLEL
1115
1116	c-----------------------------------------------------------------------
1117	c dissipation horizontale et verticale des petites echelles:
1118	c ----------------------------------------------------------
1119
1120	IF(apdiss) THEN
1121	cc$OMP PARALLEL DEFAULT(SHARED)
1122	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1123	c$OMP MASTER
1124	call suspend_timer(timer_caldyn)
1125
1126	c print*,'Entree dans la dissipation : Iteration No ',true_itau
1127	c calcul de l'energie cinetique avant dissipation
1128	c print *,'Passage dans la dissipation'
1129
1130	call VTb(VThallo)
1131	c$OMP END MASTER
1132
1133	c$OMP BARRIER
1134
1135	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1136	* jj_Nb_dissip,1,1,Request_dissip)
1137
1138	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1139	* jj_Nb_dissip,1,1,Request_dissip)
1140
1141	call Register_SwapField(teta,teta,ip1jmp1,llm,
1142	* jj_Nb_dissip,Request_dissip)
1143
1144	call Register_SwapField(p,p,ip1jmp1,llmp1,
1145	* jj_Nb_dissip,Request_dissip)
1146
1147	call Register_SwapField(pk,pk,ip1jmp1,llm,
1148	* jj_Nb_dissip,Request_dissip)
1149
1150	call SendRequest(Request_dissip)
1151	c$OMP BARRIER
1152	call WaitRequest(Request_dissip)
1153
1154	c$OMP BARRIER
1155	c$OMP MASTER
1156	call SetDistrib(jj_Nb_dissip)
1157	call VTe(VThallo)
1158	call VTb(VTdissipation)
1159	call start_timer(timer_dissip)
1160	c$OMP END MASTER
1161	c$OMP BARRIER
1162
1163	call covcont_p(llm,ucov,vcov,ucont,vcont)
1164	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1165
1166	c dissipation
1167	! ADAPTATION GCM POUR CP(T)
1168	call tpot2t_p(ip1jmp1,llm,teta,temp,pk)
1169
1170	! CALL FTRACE_REGION_BEGIN("dissip")
1171	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1172	! CALL FTRACE_REGION_END("dissip")
1173
1174	ijb=ij_begin
1175	ije=ij_end
1176	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1177	DO l=1,llm
1178	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1179	dudis(ijb:ije,l)=dudis(ijb:ije,l)/dtdiss ! passage en (m/s)/s
1180	ENDDO
1181	c$OMP END DO NOWAIT
1182	if (pole_sud) ije=ije-iip1
1183	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1184	DO l=1,llm
1185	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1186	dvdis(ijb:ije,l)=dvdis(ijb:ije,l)/dtdiss ! passage en (m/s)/s
1187	ENDDO
1188	c$OMP END DO NOWAIT
1189
1190
1191	c------------------------------------------------------------------------
1192	if (dissip_conservative) then
1193	C On rajoute la tendance due a la transform. Ec -> E therm. cree
1194	C lors de la dissipation
1195	c$OMP BARRIER
1196	c$OMP MASTER
1197	call suspend_timer(timer_dissip)
1198	call VTb(VThallo)
1199	c$OMP END MASTER
1200	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1201	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1202	call SendRequest(Request_Dissip)
1203	c$OMP BARRIER
1204	call WaitRequest(Request_Dissip)
1205	c$OMP MASTER
1206	call VTe(VThallo)
1207	call resume_timer(timer_dissip)
1208	c$OMP END MASTER
1209	c$OMP BARRIER
1210	call covcont_p(llm,ucov,vcov,ucont,vcont)
1211	call enercin_p(vcov,ucov,vcont,ucont,ecin)
1212
1213	ijb=ij_begin
1214	ije=ij_end
1215	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1216	do l=1,llm
1217	do ij=ijb,ije
1218	! ADAPTATION GCM POUR CP(T)
1219	! dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1220	! dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1221	temp(ij,l)=temp(ij,l) +
1222	& (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
1223	enddo
1224	enddo
1225	c$OMP END DO
1226	call t2tpot_p(ip1jmp1,llm,temp,ztetaec,pk)
1227	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1228	do l=1,llm
1229	do ij=ijb,ije
1230	dtetaecdt(ij,l)=ztetaec(ij,l)-teta(ij,l)
1231	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1232	enddo
1233	enddo
1234	c$OMP END DO NOWAIT
1235	endif ! of if (dissip_conservative)
1236
1237	ijb=ij_begin
1238	ije=ij_end
1239	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1240	do l=1,llm
1241	do ij=ijb,ije
1242	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1243	dtetadis(ij,l)=dtetadis(ij,l)/dtdiss ! passage en K/s
1244	enddo
1245	enddo
1246	c$OMP END DO NOWAIT
1247	c------------------------------------------------------------------------
1248
1249
1250	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1251	c ... Calcul de la valeur moyenne, unique de h aux poles .....
1252	c
1253
1254	ijb=ij_begin
1255	ije=ij_end
1256
1257	if (pole_nord) then
1258	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1259	DO l = 1, llm
1260	DO ij = 1,iim
1261	tppn(ij) = aire( ij ) * teta( ij ,l)
1262	ENDDO
1263	tpn = SSUM(iim,tppn,1)/apoln
1264
1265	DO ij = 1, iip1
1266	teta( ij ,l) = tpn
1267	ENDDO
1268	ENDDO
1269	c$OMP END DO NOWAIT
1270
1271	c$OMP MASTER
1272	DO ij = 1,iim
1273	tppn(ij) = aire( ij ) * ps ( ij )
1274	ENDDO
1275	tpn = SSUM(iim,tppn,1)/apoln
1276
1277	DO ij = 1, iip1
1278	ps( ij ) = tpn
1279	ENDDO
1280	c$OMP END MASTER
1281	endif
1282
1283	if (pole_sud) then
1284	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1285	DO l = 1, llm
1286	DO ij = 1,iim
1287	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1288	ENDDO
1289	tps = SSUM(iim,tpps,1)/apols
1290
1291	DO ij = 1, iip1
1292	teta(ij+ip1jm,l) = tps
1293	ENDDO
1294	ENDDO
1295	c$OMP END DO NOWAIT
1296
1297	c$OMP MASTER
1298	DO ij = 1,iim
1299	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1300	ENDDO
1301	tps = SSUM(iim,tpps,1)/apols
1302
1303	DO ij = 1, iip1
1304	ps(ij+ip1jm) = tps
1305	ENDDO
1306	c$OMP END MASTER
1307	endif
1308
1309
1310	c$OMP BARRIER
1311	c$OMP MASTER
1312	call VTe(VTdissipation)
1313
1314	call stop_timer(timer_dissip)
1315
1316	call VTb(VThallo)
1317	c$OMP END MASTER
1318	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1319	* jj_Nb_caldyn,Request_dissip)
1320
1321	call Register_SwapField(vcov,vcov,ip1jm,llm,
1322	* jj_Nb_caldyn,Request_dissip)
1323
1324	call Register_SwapField(teta,teta,ip1jmp1,llm,
1325	* jj_Nb_caldyn,Request_dissip)
1326
1327	call Register_SwapField(p,p,ip1jmp1,llmp1,
1328	* jj_Nb_caldyn,Request_dissip)
1329
1330	call Register_SwapField(pk,pk,ip1jmp1,llm,
1331	* jj_Nb_caldyn,Request_dissip)
1332
1333	call SendRequest(Request_dissip)
1334	c$OMP BARRIER
1335	call WaitRequest(Request_dissip)
1336
1337	c$OMP BARRIER
1338	c$OMP MASTER
1339	call SetDistrib(jj_Nb_caldyn)
1340	call VTe(VThallo)
1341	call resume_timer(timer_caldyn)
1342	c print *,'fin dissipation'
1343	c$OMP END MASTER
1344	c$OMP BARRIER
1345	END IF ! of IF(apdiss)
1346
1347	cc$OMP END PARALLEL
1348
1349	c ajout debug
1350	c IF( lafin ) then
1351	c abort_message = 'Simulation finished'
1352	c call abort_gcm(modname,abort_message,0)
1353	c ENDIF
1354
1355	c ********************************************************************
1356	c ********************************************************************
1357	c .... fin de l'integration dynamique et physique pour le pas itau ..
1358	c ********************************************************************
1359	c ********************************************************************
1360
1361	c preparation du pas d'integration suivant ......
1362	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1363	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1364	c$OMP MASTER
1365	call stop_timer(timer_caldyn)
1366	c$OMP END MASTER
1367	IF (itau==itaumax) then
1368	c$OMP MASTER
1369	call allgather_timer_average
1370
1371	if (mpi_rank==0) then
1372
1373	print ,'********************************'
1374	print ,'*** TIMER CALDYN ****'
1375	do i=0,mpi_size-1
1376	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1377	& ' : temps moyen :',
1378	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1379	enddo
1380
1381	print ,'********************************'
1382	print ,'*** TIMER VANLEER ****'
1383	do i=0,mpi_size-1
1384	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1385	& ' : temps moyen :',
1386	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1387	enddo
1388
1389	print ,'********************************'
1390	print ,'*** TIMER DISSIP ****'
1391	do i=0,mpi_size-1
1392	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1393	& ' : temps moyen :',
1394	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1395	enddo
1396
1397	print ,'********************************'
1398	print ,'*** TIMER PHYSIC ****'
1399	do i=0,mpi_size-1
1400	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1401	& ' : temps moyen :',
1402	& timer_average(jj_nb_physic(i),timer_physic,i)
1403	enddo
1404
1405	endif
1406
1407	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1408	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1409	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1410	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1411	CALL print_filtre_timer
1412	call fin_getparam
1413	call finalize_parallel
1414	c$OMP END MASTER
1415	c$OMP BARRIER
1416	RETURN
1417	ENDIF
1418
1419	IF ( .NOT.purmats ) THEN
1420	c ........................................................
1421	c .............. schema matsuno + leapfrog ..............
1422	c ........................................................
1423
1424	IF(forward. OR. leapf) THEN
1425	itau= itau + 1
1426	! iday= day_ini+itau/day_step
1427	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1428	! IF(time.GT.1.) THEN
1429	! time = time-1.
1430	! iday = iday+1
1431	! ENDIF
1432	ENDIF
1433
1434
1435	IF( itau. EQ. itaufinp1 ) then
1436
1437	if (flag_verif) then
1438	write(79,*) 'ucov',ucov
1439	write(80,*) 'vcov',vcov
1440	write(81,*) 'teta',teta
1441	write(82,*) 'ps',ps
1442	write(83,*) 'q',q
1443	WRITE(85,*) 'q1 = ',q(:,:,1)
1444	WRITE(86,*) 'q3 = ',q(:,:,3)
1445	endif
1446
1447
1448	c$OMP MASTER
1449	call fin_getparam
1450	call finalize_parallel
1451	c$OMP END MASTER
1452	abort_message = 'Simulation finished'
1453	call abort_gcm(modname,abort_message,0)
1454	RETURN
1455	ENDIF
1456	c-----------------------------------------------------------------------
1457	c ecriture du fichier histoire moyenne:
1458	c -------------------------------------
1459
1460	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1461	c$OMP BARRIER
1462	IF(itau.EQ.itaufin) THEN
1463	iav=1
1464	ELSE
1465	iav=0
1466	ENDIF
1467	#ifdef CPP_IOIPSL
1468	IF (ok_dynzon) THEN
1469	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1470	call SendRequest(TestRequest)
1471	c$OMP BARRIER
1472	call WaitRequest(TestRequest)
1473	c$OMP BARRIER
1474	c$OMP MASTER
1475	! CALL writedynav_p(histaveid, itau,vcov ,
1476	! , ucov,teta,pk,phi,q,masse,ps,phis)
1477
1478	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
1479	! CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1480	! , ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1481	c les traceurs ne sont pas sortis, trop lourd.
1482	c Peut changer eventuellement si besoin.
1483	CALL bilan_dyn_p(dtvriperiod,dtvrday_step*periodav,
1484	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
1485	& du,dudis,dutop,dufi)
1486	c$OMP END MASTER
1487	ENDIF !ok_dynzon
1488	#endif
1489	IF (ok_dyn_ave) THEN
1490	!$OMP MASTER
1491	#ifdef CPP_IOIPSL
1492	! Ehouarn: Gather fields and make master send to output
1493	call Gather_Field(vcov,ip1jm,llm,0)
1494	call Gather_Field(ucov,ip1jmp1,llm,0)
1495	call Gather_Field(teta,ip1jmp1,llm,0)
1496	call Gather_Field(pk,ip1jmp1,llm,0)
1497	call Gather_Field(phi,ip1jmp1,llm,0)
1498	do iq=1,nqtot
1499	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1500	enddo
1501	call Gather_Field(masse,ip1jmp1,llm,0)
1502	call Gather_Field(ps,ip1jmp1,1,0)
1503	call Gather_Field(phis,ip1jmp1,1,0)
1504	if (mpi_rank==0) then
1505	CALL writedynav(itau,vcov,
1506	& ucov,teta,pk,phi,q,masse,ps,phis)
1507	endif
1508	#endif
1509	!$OMP END MASTER
1510	ENDIF ! of IF (ok_dyn_ave)
1511	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1512
1513	c-----------------------------------------------------------------------
1514	c ecriture de la bande histoire:
1515	c ------------------------------
1516
1517	IF( MOD(itau,iecri).EQ.0) THEN
1518	! Ehouarn: output only during LF or Backward Matsuno
1519	if (leapf.or.(.not.leapf.and.(.not.forward))) then
1520	c$OMP BARRIER
1521	c$OMP MASTER
1522	nbetat = nbetatdem
1523	c$OMP END MASTER
1524	c$OMP BARRIER
1525
1526	! ADAPTATION GCM POUR CP(T)
1527	call tpot2t_p(ip1jmp1,llm,teta,temp,pk)
1528	ijb=ij_begin
1529	ije=ij_end
1530	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1531	do l=1,llm
1532	tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/pk(ijb:ije,l)
1533	enddo
1534	!$OMP END DO
1535
1536	!$OMP MASTER
1537	! CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1538	CALL geopot_p ( ip1jmp1, tsurpk , pk , pks, phis , phi )
1539
1540	cym unat=0.
1541
1542	ijb=ij_begin
1543	ije=ij_end
1544
1545	if (pole_nord) then
1546	ijb=ij_begin+iip1
1547	unat(1:iip1,:)=0.
1548	endif
1549
1550	if (pole_sud) then
1551	ije=ij_end-iip1
1552	unat(ij_end-iip1+1:ij_end,:)=0.
1553	endif
1554
1555	do l=1,llm
1556	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1557	enddo
1558
1559	ijb=ij_begin
1560	ije=ij_end
1561	if (pole_sud) ije=ij_end-iip1
1562
1563	do l=1,llm
1564	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1565	enddo
1566
1567	#ifdef CPP_IOIPSL
1568	if (ok_dyn_ins) then
1569	! Ehouarn: Gather fields and make master write to output
1570	call Gather_Field(vcov,ip1jm,llm,0)
1571	call Gather_Field(ucov,ip1jmp1,llm,0)
1572	call Gather_Field(teta,ip1jmp1,llm,0)
1573	call Gather_Field(phi,ip1jmp1,llm,0)
1574	do iq=1,nqtot
1575	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1576	enddo
1577	call Gather_Field(masse,ip1jmp1,llm,0)
1578	call Gather_Field(ps,ip1jmp1,1,0)
1579	call Gather_Field(phis,ip1jmp1,1,0)
1580	if (mpi_rank==0) then
1581	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1582	endif
1583	! CALL writehist_p(histid,histvid, itau,vcov,
1584	! & ucov,teta,phi,q,masse,ps,phis)
1585	! or use writefield_p
1586	! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1587	! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1588	! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1589	! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1590	endif ! of if (ok_dyn_ins)
1591	#endif
1592	! For some Grads outputs of fields
1593	if (output_grads_dyn) then
1594	! Ehouarn: hope this works the way I think it does:
1595	call Gather_Field(unat,ip1jmp1,llm,0)
1596	call Gather_Field(vnat,ip1jm,llm,0)
1597	call Gather_Field(teta,ip1jmp1,llm,0)
1598	call Gather_Field(ps,ip1jmp1,1,0)
1599	do iq=1,nqtot
1600	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1601	enddo
1602	if (mpi_rank==0) then
1603	#include "write_grads_dyn.h"
1604	endif
1605	endif ! of if (output_grads_dyn)
1606	!$OMP END MASTER
1607	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1608	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1609
1610	IF(itau.EQ.itaufin) THEN
1611
1612	c$OMP BARRIER
1613	c$OMP MASTER
1614
1615	if (planet_type.eq."mars") then
1616	! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
1617	abort_message = 'dynredem1_mars A FAIRE'
1618	call abort_gcm(modname,abort_message,0)
1619	else
1620	! Write an Earth-format restart file
1621	CALL dynredem1_p("restart.nc",0.0,
1622	& vcov,ucov,teta,q,masse,ps)
1623	endif ! of if (planet_type.eq."mars")
1624
1625	! CLOSE(99)
1626	c$OMP END MASTER
1627	ENDIF ! of IF (itau.EQ.itaufin)
1628
1629	c-----------------------------------------------------------------------
1630	c gestion de l'integration temporelle:
1631	c ------------------------------------
1632
1633	IF( MOD(itau,iperiod).EQ.0 ) THEN
1634	GO TO 1
1635	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1636
1637	IF( forward ) THEN
1638	c fin du pas forward et debut du pas backward
1639
1640	forward = .FALSE.
1641	leapf = .FALSE.
1642	GO TO 2
1643
1644	ELSE
1645	c fin du pas backward et debut du premier pas leapfrog
1646
1647	leapf = .TRUE.
1648	dt = 2.*dtvr
1649	GO TO 2
1650	END IF
1651	ELSE
1652
1653	c ...... pas leapfrog .....
1654
1655	leapf = .TRUE.
1656	dt = 2.*dtvr
1657	GO TO 2
1658	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1659	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1660
1661
1662	ELSE ! of IF (.not.purmats)
1663
1664	c ........................................................
1665	c .............. schema matsuno ...............
1666	c ........................................................
1667	IF( forward ) THEN
1668
1669	itau = itau + 1
1670	! iday = day_ini+itau/day_step
1671	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1672	!
1673	! IF(time.GT.1.) THEN
1674	! time = time-1.
1675	! iday = iday+1
1676	! ENDIF
1677
1678	forward = .FALSE.
1679	IF( itau. EQ. itaufinp1 ) then
1680	c$OMP MASTER
1681	call fin_getparam
1682	call finalize_parallel
1683	c$OMP END MASTER
1684	abort_message = 'Simulation finished'
1685	call abort_gcm(modname,abort_message,0)
1686	RETURN
1687	ENDIF
1688	GO TO 2
1689
1690	ELSE ! of IF(forward) i.e. backward step
1691
1692	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1693	IF(itau.EQ.itaufin) THEN
1694	iav=1
1695	ELSE
1696	iav=0
1697	ENDIF
1698	#ifdef CPP_IOIPSL
1699	IF (ok_dynzon) THEN
1700	c$OMP BARRIER
1701	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1702	call SendRequest(TestRequest)
1703	c$OMP BARRIER
1704	call WaitRequest(TestRequest)
1705	c$OMP BARRIER
1706	c$OMP MASTER
1707	! CALL writedynav_p(histaveid, itau,vcov ,
1708	! , ucov,teta,pk,phi,q,masse,ps,phis)
1709	! CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1710	! , ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1711	c les traceurs ne sont pas sortis, trop lourd.
1712	c Peut changer eventuellement si besoin.
1713	CALL bilan_dyn_p(dtvriperiod,dtvrday_step*periodav,
1714	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
1715	& du,dudis,dutop,dufi)
1716
1717	c$OMP END MASTER
1718	END IF !ok_dynzon
1719	#endif
1720	IF (ok_dyn_ave) THEN
1721	!$OMP MASTER
1722	#ifdef CPP_IOIPSL
1723	! Ehouarn: Gather fields and make master send to output
1724	call Gather_Field(vcov,ip1jm,llm,0)
1725	call Gather_Field(ucov,ip1jmp1,llm,0)
1726	call Gather_Field(teta,ip1jmp1,llm,0)
1727	call Gather_Field(pk,ip1jmp1,llm,0)
1728	call Gather_Field(phi,ip1jmp1,llm,0)
1729	do iq=1,nqtot
1730	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1731	enddo
1732	call Gather_Field(masse,ip1jmp1,llm,0)
1733	call Gather_Field(ps,ip1jmp1,1,0)
1734	call Gather_Field(phis,ip1jmp1,1,0)
1735	if (mpi_rank==0) then
1736	CALL writedynav(itau,vcov,
1737	& ucov,teta,pk,phi,q,masse,ps,phis)
1738	endif
1739	#endif
1740	!$OMP END MASTER
1741	ENDIF ! of IF (ok_dyn_ave)
1742
1743	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1744
1745
1746	IF(MOD(itau,iecri ).EQ.0) THEN
1747	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1748	c$OMP BARRIER
1749	c$OMP MASTER
1750	nbetat = nbetatdem
1751	c$OMP END MASTER
1752	c$OMP BARRIER
1753
1754	! ADAPTATION GCM POUR CP(T)
1755	call tpot2t_p(ip1jmp1,llm,teta,temp,pk)
1756	ijb=ij_begin
1757	ije=ij_end
1758	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1759	do l=1,llm
1760	tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/
1761	& pk(ijb:ije,l)
1762	enddo
1763	!$OMP END DO
1764
1765	!$OMP MASTER
1766	! CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1767	CALL geopot_p(ip1jmp1,tsurpk,pk,pks,phis,phi)
1768
1769	cym unat=0.
1770	ijb=ij_begin
1771	ije=ij_end
1772
1773	if (pole_nord) then
1774	ijb=ij_begin+iip1
1775	unat(1:iip1,:)=0.
1776	endif
1777
1778	if (pole_sud) then
1779	ije=ij_end-iip1
1780	unat(ij_end-iip1+1:ij_end,:)=0.
1781	endif
1782
1783	do l=1,llm
1784	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1785	enddo
1786
1787	ijb=ij_begin
1788	ije=ij_end
1789	if (pole_sud) ije=ij_end-iip1
1790
1791	do l=1,llm
1792	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1793	enddo
1794
1795	#ifdef CPP_IOIPSL
1796	if (ok_dyn_ins) then
1797	! Ehouarn: Gather fields and make master send to output
1798	call Gather_Field(vcov,ip1jm,llm,0)
1799	call Gather_Field(ucov,ip1jmp1,llm,0)
1800	call Gather_Field(teta,ip1jmp1,llm,0)
1801	call Gather_Field(phi,ip1jmp1,llm,0)
1802	do iq=1,nqtot
1803	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1804	enddo
1805	call Gather_Field(masse,ip1jmp1,llm,0)
1806	call Gather_Field(ps,ip1jmp1,1,0)
1807	call Gather_Field(phis,ip1jmp1,1,0)
1808	if (mpi_rank==0) then
1809	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1810	endif
1811	! CALL writehist_p(histid, histvid, itau,vcov ,
1812	! & ucov,teta,phi,q,masse,ps,phis)
1813	endif ! of if (ok_dyn_ins)
1814	#endif
1815	! For some Grads output (but does it work?)
1816	if (output_grads_dyn) then
1817	call Gather_Field(unat,ip1jmp1,llm,0)
1818	call Gather_Field(vnat,ip1jm,llm,0)
1819	call Gather_Field(teta,ip1jmp1,llm,0)
1820	call Gather_Field(ps,ip1jmp1,1,0)
1821	do iq=1,nqtot
1822	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1823	enddo
1824	c
1825	if (mpi_rank==0) then
1826	#include "write_grads_dyn.h"
1827	endif
1828	endif ! of if (output_grads_dyn)
1829
1830	!$OMP END MASTER
1831	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1832
1833	IF(itau.EQ.itaufin) THEN
1834	if (planet_type.eq."mars") then
1835	! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
1836	abort_message = 'dynredem1_mars A FAIRE'
1837	call abort_gcm(modname,abort_message,0)
1838	else
1839	c$OMP MASTER
1840	CALL dynredem1_p("restart.nc",0.0,
1841	. vcov,ucov,teta,q,masse,ps)
1842	c$OMP END MASTER
1843	endif ! of if (planet_type.eq."mars")
1844	ENDIF ! of IF(itau.EQ.itaufin)
1845
1846	forward = .TRUE.
1847	GO TO 1
1848
1849	ENDIF ! of IF (forward)
1850
1851	END IF ! of IF(.not.purmats)
1852	c$OMP MASTER
1853	call fin_getparam
1854	call finalize_parallel
1855	c$OMP END MASTER
1856	RETURN
1857	END

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

Download in other formats:

Original Format