Context Navigation

leapfrog_p.F @ 1533

Last change on this file since 1533 was 1272, checked in by Ehouarn Millour, 15 years ago

Bug fix (problem only occured when running Newtonian case with both MPI and OpenMP turned on): call SendRequest?(Request_Physic) must be followed by an OMP barrier.

EM & YM

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 48.9 KB

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

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

Download in other formats:

Original Format