Context Navigation

leapfrog_p.F @ 154

Last change on this file since 154 was 130, checked in by slebonnois, 14 years ago

Sebastien:

correction de bugs dans phytitan suite a compil avec -debug sur gnome
elimination de nbetat* dans gcm.F et leapfrog.F (et dans le parallele)
correction rday -> rjour dans sortvarc.F
ajustement des arch-GNOME*

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
74	c variables dynamiques
75	REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
76	REAL :: teta(ip1jmp1,llm) ! temperature potentielle
77	REAL :: q(ip1jmp1,llm,nqtot) ! champs advectes
78	REAL :: ps(ip1jmp1) ! pression au sol
79	REAL,SAVE :: p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
80	REAL,SAVE :: pks(ip1jmp1) ! exner au sol
81	REAL,SAVE :: pk(ip1jmp1,llm) ! exner au milieu des couches
82	REAL,SAVE :: pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
83	REAL :: masse(ip1jmp1,llm) ! masse d'air
84	REAL :: phis(ip1jmp1) ! geopotentiel au sol
85	REAL,SAVE :: phi(ip1jmp1,llm) ! geopotentiel
86	REAL,SAVE :: w(ip1jmp1,llm) ! vitesse verticale
87	! ADAPTATION GCM POUR CP(T)
88	REAL,SAVE :: temp(ip1jmp1,llm) ! temperature
89	REAL,SAVE :: tsurpk(ip1jmp1,llm) ! cpp*T/pk
90
91	c variables dynamiques intermediaire pour le transport
92	REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
93
94	c variables dynamiques au pas -1
95	REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
96	REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
97	REAL,SAVE :: massem1(ip1jmp1,llm)
98
99	c tendances dynamiques
100	REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
101	REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
102	REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
103
104	c tendances de la dissipation
105	REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
106	REAL,SAVE :: dtetadis(ip1jmp1,llm)
107
108	c tendances physiques
109	REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
110	REAL,SAVE :: dtetafi(ip1jmp1,llm)
111	REAL,SAVE :: dpfi(ip1jmp1)
112	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
113
114	c tendances top_bound (sponge layer)
115	REAL,SAVE :: dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
116	REAL,SAVE :: dtetatop(ip1jmp1,llm)
117	REAL,SAVE :: dptop(ip1jmp1)
118	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqtop
119
120	c variables pour le fichier histoire
121	REAL dtav ! intervalle de temps elementaire
122
123	REAL tppn(iim),tpps(iim),tpn,tps
124	c
125	INTEGER itau,itaufinp1,iav
126	! INTEGER iday ! jour julien
127	REAL time
128
129	REAL SSUM
130	REAL time_0
131	REAL,SAVE :: finvmaold(ip1jmp1,llm)
132
133	cym LOGICAL lafin
134	LOGICAL :: lafin
135	INTEGER ij,iq,l
136	INTEGER ik
137
138	real time_step, t_wrt, t_ops
139
140	! jD_cur: jour julien courant
141	! jH_cur: heure julienne courante
142	REAL :: jD_cur, jH_cur
143	INTEGER :: an, mois, jour
144	REAL :: secondes
145
146	LOGICAL first,callinigrads
147
148	data callinigrads/.true./
149	character*10 string10
150
151	REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
152	REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
153
154	c+jld variables test conservation energie
155	REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
156	C Tendance de la temp. potentiel d (theta)/ d t due a la
157	C tansformation d'energie cinetique en energie thermique
158	C cree par la dissipation
159	REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
160	REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
161	REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
162	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
163	CHARACTER*15 ztit
164	! INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
165	! SAVE ip_ebil_dyn
166	! DATA ip_ebil_dyn/0/
167	c-jld
168
169	character*80 dynhist_file, dynhistave_file
170	character*20 modname
171	character*80 abort_message
172
173
174	logical,PARAMETER :: dissip_conservative=.TRUE.
175
176	INTEGER testita
177	PARAMETER (testita = 9)
178
179	logical , parameter :: flag_verif = .false.
180
181	! for CP(T) -- Aymeric
182	real :: dtec
183	real,external :: cpdet
184	real,save :: ztetaec(ip1jmp1,llm) !!SAVE ???
185
186	c declaration liees au parallelisme
187	INTEGER :: ierr
188	LOGICAL :: FirstCaldyn
189	LOGICAL :: FirstPhysic
190	INTEGER :: ijb,ije,j,i
191	type(Request) :: TestRequest
192	type(Request) :: Request_Dissip
193	type(Request) :: Request_physic
194	REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
195	REAL,SAVE :: dtetafi_tmp(iip1,llm)
196	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
197	REAL,SAVE :: dpfi_tmp(iip1)
198
199	INTEGER :: true_itau
200	LOGICAL :: verbose=.true.
201	INTEGER :: iapptrac
202	INTEGER :: AdjustCount
203	! INTEGER :: var_time
204	LOGICAL :: ok_start_timer=.FALSE.
205	LOGICAL, SAVE :: firstcall=.TRUE.
206
207	c dummy: sinon cette routine n'est jamais compilee...
208	if(1.eq.0) then
209	CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
210	endif
211
212	c$OMP MASTER
213	ItCount=0
214	c$OMP END MASTER
215	true_itau=0
216	FirstCaldyn=.TRUE.
217	FirstPhysic=.TRUE.
218	iapptrac=0
219	AdjustCount = 0
220	lafin=.false.
221
222	itaufin = nday*day_step
223	if (less1day) then
224	c MODIF VENUS: to run less than one day:
225	itaufin = int(fractday*day_step)
226	endif
227	itaufinp1 = itaufin +1
228	modname="leapfrog_p"
229
230	itau = 0
231	! iday = day_ini+itau/day_step
232	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
233	! IF(time.GT.1.) THEN
234	! time = time-1.
235	! iday = iday+1
236	! ENDIF
237
238	c Allocate variables depending on dynamic variable nqtot
239	c$OMP MASTER
240	IF (firstcall) THEN
241	firstcall=.FALSE.
242	ALLOCATE(dq(ip1jmp1,llm,nqtot))
243	ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
244	ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
245	ALLOCATE(dqtop(ip1jmp1,llm,nqtot))
246	END IF
247	c$OMP END MASTER
248	c$OMP BARRIER
249
250	c-----------------------------------------------------------------------
251	c On initialise la pression et la fonction d'Exner :
252	c --------------------------------------------------
253
254	c$OMP MASTER
255	c INITIALISATIONS
256	dudis(:,:) =0.
257	dvdis(:,:) =0.
258	dtetadis(:,:)=0.
259	dutop(:,:) =0.
260	dvtop(:,:) =0.
261	dtetatop(:,:)=0.
262	dqtop(:,:,:) =0.
263	dptop(:) =0.
264	dufi(:,:) =0.
265	dvfi(:,:) =0.
266	dtetafi(:,:)=0.
267	dqfi(:,:,:) =0.
268	dpfi(:) =0.
269	dq(:,:,:)=0.
270
271	CALL pression ( ip1jmp1, ap, bp, ps, p )
272	if (disvert_type==1) then
273	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
274	else ! we assume that we are in the disvert_type==2 case
275	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
276	endif
277	c$OMP END MASTER
278	c-----------------------------------------------------------------------
279	c Debut de l'integration temporelle:
280	c ----------------------------------
281	c et du parallelisme !!
282
283	1 CONTINUE
284
285	jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
286	jH_cur = jH_ref + &
287	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
288
289
290	#ifdef CPP_IOIPSL
291	if (ok_guide) then
292	!$OMP MASTER
293	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
294	!$OMP END MASTER
295	!$OMP BARRIER
296	endif
297	#endif
298
299	c
300	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
301	c CALL test_period ( ucov,vcov,teta,q,p,phis )
302	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
303	c ENDIF
304	c
305	cym CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
306	cym CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
307	cym CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
308	cym CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
309	cym CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
310
311	if (FirstCaldyn) then
312	c$OMP MASTER
313	ucovm1=ucov
314	vcovm1=vcov
315	tetam1= teta
316	massem1= masse
317	psm1= ps
318
319	finvmaold = masse
320	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
321	c$OMP END MASTER
322	c$OMP BARRIER
323	else
324	! Save fields obtained at previous time step as '...m1'
325	ijb=ij_begin
326	ije=ij_end
327
328	c$OMP MASTER
329	psm1 (ijb:ije) = ps (ijb:ije)
330	c$OMP END MASTER
331
332	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
333	DO l=1,llm
334	ije=ij_end
335	ucovm1 (ijb:ije,l) = ucov (ijb:ije,l)
336	tetam1 (ijb:ije,l) = teta (ijb:ije,l)
337	massem1 (ijb:ije,l) = masse (ijb:ije,l)
338	finvmaold(ijb:ije,l)=masse(ijb:ije,l)
339
340	if (pole_sud) ije=ij_end-iip1
341	vcovm1(ijb:ije,l) = vcov (ijb:ije,l)
342
343
344	ENDDO
345	c$OMP ENDDO
346
347
348	CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
349	. llm, -2,2, .TRUE., 1 )
350
351	endif ! of if (FirstCaldyn)
352
353	forward = .TRUE.
354	leapf = .FALSE.
355	dt = dtvr
356
357	c ... P.Le Van .26/04/94 ....
358
359	cym CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
360	cym CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
361
362	cym ne sert a rien
363	cym call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
364
365	2 CONTINUE
366
367	c$OMP MASTER
368	ItCount=ItCount+1
369	if (MOD(ItCount,1)==1) then
370	debug=.true.
371	else
372	debug=.false.
373	endif
374	c$OMP END MASTER
375	c-----------------------------------------------------------------------
376
377	c date:
378	c -----
379
380
381	c gestion des appels de la physique et des dissipations:
382	c ------------------------------------------------------
383	c
384	c ... P.Le Van ( 6/02/95 ) ....
385
386	apphys = .FALSE.
387	statcl = .FALSE.
388	conser = .FALSE.
389	apdiss = .FALSE.
390	c idissip=1
391	IF( purmats ) THEN
392	! Purely Matsuno time stepping
393	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
394	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
395	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
396	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
397	ELSE
398	! Leapfrog/Matsuno time stepping
399	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
400	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
401	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
402	END IF
403
404	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
405	! supress dissipation step
406	if (llm.eq.1) then
407	apdiss=.false.
408	endif
409
410	cym ---> Pour le moment
411	cym apphys = .FALSE.
412	statcl = .FALSE.
413	conser = .FALSE. ! ie: no output of control variables to stdout in //
414
415	if (firstCaldyn) then
416	c$OMP MASTER
417	call SetDistrib(jj_Nb_Caldyn)
418	c$OMP END MASTER
419	c$OMP BARRIER
420	firstCaldyn=.FALSE.
421	cym call InitTime
422	c$OMP MASTER
423	call Init_timer
424	c$OMP END MASTER
425	endif
426
427	c$OMP MASTER
428	IF (ok_start_timer) THEN
429	CALL InitTime
430	ok_start_timer=.FALSE.
431	ENDIF
432	c$OMP END MASTER
433
434	if (Adjust) then
435	c$OMP MASTER
436	AdjustCount=AdjustCount+1
437	if (iapptrac==iapp_tracvl .and. (forward. OR . leapf)
438	& .and. itau/iphysiq>2 .and. Adjustcount>30) then
439	AdjustCount=0
440	call allgather_timer_average
441
442	if (Verbose) then
443
444	print ,'********************************'
445	print ,'*** TIMER CALDYN ****'
446	do i=0,mpi_size-1
447	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
448	& ' : temps moyen :',
449	& timer_average(jj_nb_caldyn(i),timer_caldyn,i),
450	& '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
451	enddo
452
453	print ,'********************************'
454	print ,'*** TIMER VANLEER ****'
455	do i=0,mpi_size-1
456	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
457	& ' : temps moyen :',
458	& timer_average(jj_nb_vanleer(i),timer_vanleer,i),
459	& '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
460	enddo
461
462	print ,'********************************'
463	print ,'*** TIMER DISSIP ****'
464	do i=0,mpi_size-1
465	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
466	& ' : temps moyen :',
467	& timer_average(jj_nb_dissip(i),timer_dissip,i),
468	& '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
469	enddo
470
471	if (mpi_rank==0) call WriteBands
472
473	endif
474
475	call AdjustBands_caldyn
476	if (mpi_rank==0) call WriteBands
477
478	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
479	& jj_Nb_caldyn,0,0,TestRequest)
480	call Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
481	& jj_Nb_caldyn,0,0,TestRequest)
482	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
483	& jj_Nb_caldyn,0,0,TestRequest)
484	call Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
485	& jj_Nb_caldyn,0,0,TestRequest)
486	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
487	& jj_Nb_caldyn,0,0,TestRequest)
488	call Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
489	& jj_Nb_caldyn,0,0,TestRequest)
490	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
491	& jj_Nb_caldyn,0,0,TestRequest)
492	call Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
493	& jj_Nb_caldyn,0,0,TestRequest)
494	call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
495	& jj_Nb_caldyn,0,0,TestRequest)
496	call Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
497	& jj_Nb_caldyn,0,0,TestRequest)
498	call Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
499	& jj_Nb_caldyn,0,0,TestRequest)
500	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
501	& jj_Nb_caldyn,0,0,TestRequest)
502	call Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
503	& jj_Nb_caldyn,0,0,TestRequest)
504	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
505	& jj_Nb_caldyn,0,0,TestRequest)
506	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
507	& jj_Nb_caldyn,0,0,TestRequest)
508	call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
509	& jj_Nb_caldyn,0,0,TestRequest)
510
511	do j=1,nqtot
512	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
513	& jj_nb_caldyn,0,0,TestRequest)
514	enddo
515	! ADAPTATION GCM POUR CP(T)
516	call Register_SwapFieldHallo(temp,temp,ip1jmp1,llm,
517	& jj_Nb_caldyn,0,0,TestRequest)
518	call Register_SwapFieldHallo(tsurpk,tsurpk,ip1jmp1,llm,
519	& jj_Nb_caldyn,0,0,TestRequest)
520
521	call SetDistrib(jj_nb_caldyn)
522	call SendRequest(TestRequest)
523	call WaitRequest(TestRequest)
524
525	call AdjustBands_dissip
526	call AdjustBands_physic
527
528	endif
529	c$OMP END MASTER
530	endif
531
532
533
534	c-----------------------------------------------------------------------
535	c calcul des tendances dynamiques:
536	c --------------------------------
537	c$OMP BARRIER
538	c$OMP MASTER
539	call VTb(VThallo)
540	c$OMP END MASTER
541
542	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
543	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
544	call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
545	call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
546	call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
547	call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
548	call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
549	call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
550	! ADAPTATION GCM POUR CP(T)
551	call Register_Hallo(temp,ip1jmp1,llm,1,1,1,1,TestRequest)
552	call Register_Hallo(tsurpk,ip1jmp1,llm,1,1,1,1,TestRequest)
553
554	c do j=1,nqtot
555	c call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
556	c * TestRequest)
557	c enddo
558
559	call SendRequest(TestRequest)
560	c$OMP BARRIER
561	call WaitRequest(TestRequest)
562
563	c$OMP MASTER
564	call VTe(VThallo)
565	c$OMP END MASTER
566	c$OMP BARRIER
567
568	if (debug) then
569	!$OMP BARRIER
570	!$OMP MASTER
571	call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
572	call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
573	call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
574	call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
575	call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
576	call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
577	call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
578	call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
579	call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
580	do j=1,nqtot
581	call WriteField_p('q'//trim(int2str(j)),
582	. reshape(q(:,:,j),(/iip1,jmp1,llm/)))
583	enddo
584	!$OMP END MASTER
585	c$OMP BARRIER
586	endif
587
588
589	True_itau=True_itau+1
590
591	! ADAPTATION GCM POUR CP(T)
592	call tpot2t_p(ip1jmp1,llm,teta,temp,pk)
593	ijb=ij_begin
594	ije=ij_end
595	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
596	do l=1,llm
597	tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/pk(ijb:ije,l)
598	enddo
599	!$OMP END DO
600
601	c$OMP MASTER
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 (disvert_type==1) then
753	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
754	else ! we assume that we are in the disvert_type==2 case
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 (disvert_type==1) then
1107	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
1108	else ! we assume that we are in the disvert_type==2 case
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
1521	! ADAPTATION GCM POUR CP(T)
1522	call tpot2t_p(ip1jmp1,llm,teta,temp,pk)
1523	ijb=ij_begin
1524	ije=ij_end
1525	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1526	do l=1,llm
1527	tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/pk(ijb:ije,l)
1528	enddo
1529	!$OMP END DO
1530
1531	!$OMP MASTER
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	!$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
1745	! ADAPTATION GCM POUR CP(T)
1746	call tpot2t_p(ip1jmp1,llm,teta,temp,pk)
1747	ijb=ij_begin
1748	ije=ij_end
1749	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1750	do l=1,llm
1751	tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/
1752	& pk(ijb:ije,l)
1753	enddo
1754	!$OMP END DO
1755
1756	!$OMP MASTER
1757	! CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1758	CALL geopot_p(ip1jmp1,tsurpk,pk,pks,phis,phi)
1759
1760	cym unat=0.
1761	ijb=ij_begin
1762	ije=ij_end
1763
1764	if (pole_nord) then
1765	ijb=ij_begin+iip1
1766	unat(1:iip1,:)=0.
1767	endif
1768
1769	if (pole_sud) then
1770	ije=ij_end-iip1
1771	unat(ij_end-iip1+1:ij_end,:)=0.
1772	endif
1773
1774	do l=1,llm
1775	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1776	enddo
1777
1778	ijb=ij_begin
1779	ije=ij_end
1780	if (pole_sud) ije=ij_end-iip1
1781
1782	do l=1,llm
1783	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1784	enddo
1785
1786	#ifdef CPP_IOIPSL
1787	if (ok_dyn_ins) then
1788	! Ehouarn: Gather fields and make master send to output
1789	call Gather_Field(vcov,ip1jm,llm,0)
1790	call Gather_Field(ucov,ip1jmp1,llm,0)
1791	call Gather_Field(teta,ip1jmp1,llm,0)
1792	call Gather_Field(phi,ip1jmp1,llm,0)
1793	do iq=1,nqtot
1794	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1795	enddo
1796	call Gather_Field(masse,ip1jmp1,llm,0)
1797	call Gather_Field(ps,ip1jmp1,1,0)
1798	call Gather_Field(phis,ip1jmp1,1,0)
1799	if (mpi_rank==0) then
1800	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1801	endif
1802	! CALL writehist_p(histid, histvid, itau,vcov ,
1803	! & ucov,teta,phi,q,masse,ps,phis)
1804	endif ! of if (ok_dyn_ins)
1805	#endif
1806	! For some Grads output (but does it work?)
1807	if (output_grads_dyn) then
1808	call Gather_Field(unat,ip1jmp1,llm,0)
1809	call Gather_Field(vnat,ip1jm,llm,0)
1810	call Gather_Field(teta,ip1jmp1,llm,0)
1811	call Gather_Field(ps,ip1jmp1,1,0)
1812	do iq=1,nqtot
1813	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1814	enddo
1815	c
1816	if (mpi_rank==0) then
1817	#include "write_grads_dyn.h"
1818	endif
1819	endif ! of if (output_grads_dyn)
1820
1821	!$OMP END MASTER
1822	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1823
1824	IF(itau.EQ.itaufin) THEN
1825	if (planet_type.eq."mars") then
1826	! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
1827	abort_message = 'dynredem1_mars A FAIRE'
1828	call abort_gcm(modname,abort_message,0)
1829	else
1830	c$OMP MASTER
1831	CALL dynredem1_p("restart.nc",0.0,
1832	. vcov,ucov,teta,q,masse,ps)
1833	c$OMP END MASTER
1834	endif ! of if (planet_type.eq."mars")
1835	ENDIF ! of IF(itau.EQ.itaufin)
1836
1837	forward = .TRUE.
1838	GO TO 1
1839
1840	ENDIF ! of IF (forward)
1841
1842	END IF ! of IF(.not.purmats)
1843	c$OMP MASTER
1844	call fin_getparam
1845	call finalize_parallel
1846	c$OMP END MASTER
1847	RETURN
1848	END

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

Download in other formats:

Original Format