Context Navigation

leapfrog_p.F @ 5427

Last change on this file since 5427 was 1717, checked in by Ehouarn Millour, 12 years ago
Added "arch" files for Ada (using dynamic libraries for NetCDF, you must have export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/smplocal/pub/NetCDF/4.1.3/lib:/smplocal/pub/HDF5/1.8.9/seq/lib in your .bashrc or .bash_login or in your job to run). Also updated some sources so that gcm bench runs in "debug" mode (note that all these changes are minor and have already been implemented in LMDZ5 trunk). EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 49.3 KB

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

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

Download in other formats:

Original Format