Context Navigation

leapfrog_p.F @ 1164

Last change on this file since 1164 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 48.3 KB

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

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

Download in other formats:

Original Format