Context Navigation

leapfrog_p.F @ 1322

Last change on this file since 1322 was 1322, checked in by Laurent Fairhead, 14 years ago

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file

Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

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

File size: 49.4 KB

Line
1	!
2	! $Id: leapfrog_p.F 1322 2010-03-12 10:54:11Z fairhead $
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	USE control_mod
23
24	IMPLICIT NONE
25
26	c ...... Version du 10/01/98 ..........
27
28	c avec coordonnees verticales hybrides
29	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
30
31	c=======================================================================
32	c
33	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
34	c -------
35	c
36	c Objet:
37	c ------
38	c
39	c GCM LMD nouvelle grille
40	c
41	c=======================================================================
42	c
43	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
44	c et possibilite d'appeler une fonction f(y) a derivee tangente
45	c hyperbolique a la place de la fonction a derivee sinusoidale.
46
47	c ... Possibilite de choisir le shema pour l'advection de
48	c q , en modifiant iadv dans traceur.def (10/02) .
49	c
50	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
51	c Pour Van-Leer iadv=10
52	c
53	c-----------------------------------------------------------------------
54	c Declarations:
55	c -------------
56
57	#include "dimensions.h"
58	#include "paramet.h"
59	#include "comconst.h"
60	#include "comdissnew.h"
61	#include "comvert.h"
62	#include "comgeom.h"
63	#include "logic.h"
64	#include "temps.h"
65	#include "ener.h"
66	#include "description.h"
67	#include "serre.h"
68	#include "com_io_dyn.h"
69	#include "iniprint.h"
70	#include "academic.h"
71
72	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 = REAL(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) 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	if (prt_level >= 10) then
680	write(lunout,*)
681	& 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
682	endif
683	c$OMP END MASTER
684
685	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
686
687	c$OMP BARRIER
688	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
689	c$OMP BARRIER
690	jD_cur = jD_ref + day_ini - day_ref
691	$ + int (itau * dtvr / daysec)
692	jH_cur = jH_ref + &
693	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
694	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
695
696	c rajout debug
697	c lafin = .true.
698
699
700	c Inbterface avec les routines de phylmd (phymars ... )
701	c -----------------------------------------------------
702
703	c+jld
704
705	c Diagnostique de conservation de l'energie : initialisation
706	IF (ip_ebil_dyn.ge.1 ) THEN
707	ztit='bil dyn'
708	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
709	IF (planet_type.eq."earth") THEN
710	CALL diagedyn(ztit,2,1,1,dtphys
711	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
712	ENDIF
713	ENDIF
714	c-jld
715	c$OMP BARRIER
716	c$OMP MASTER
717	call VTb(VThallo)
718	c$OMP END MASTER
719
720	call SetTag(Request_physic,800)
721
722	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
723	* jj_Nb_physic,2,2,Request_physic)
724
725	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
726	* jj_Nb_physic,2,2,Request_physic)
727
728	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
729	* jj_Nb_physic,2,2,Request_physic)
730
731	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
732	* jj_Nb_physic,1,2,Request_physic)
733
734	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
735	* jj_Nb_physic,2,2,Request_physic)
736
737	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
738	* jj_Nb_physic,2,2,Request_physic)
739
740	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
741	* jj_Nb_physic,2,2,Request_physic)
742
743	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
744	* jj_Nb_physic,2,2,Request_physic)
745
746	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
747	* jj_Nb_physic,2,2,Request_physic)
748
749	c call SetDistrib(jj_nb_vanleer)
750	do j=1,nqtot
751
752	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
753	* jj_Nb_physic,2,2,Request_physic)
754	enddo
755
756	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
757	* jj_Nb_physic,2,2,Request_physic)
758
759	call SendRequest(Request_Physic)
760	c$OMP BARRIER
761	call WaitRequest(Request_Physic)
762
763	c$OMP BARRIER
764	c$OMP MASTER
765	call SetDistrib(jj_nb_Physic)
766	call VTe(VThallo)
767
768	call VTb(VTphysiq)
769	c$OMP END MASTER
770	c$OMP BARRIER
771
772	cc$OMP MASTER
773	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
774	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
775	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
776	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
777	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
778	cc$OMP END MASTER
779	cc$OMP BARRIER
780	! CALL FTRACE_REGION_BEGIN("calfis")
781	CALL calfis_p(lafin ,jD_cur, jH_cur,
782	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
783	$ du,dv,dteta,dq,
784	$ flxw,
785	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
786	! CALL FTRACE_REGION_END("calfis")
787	ijb=ij_begin
788	ije=ij_end
789	if ( .not. pole_nord) then
790	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
791	DO l=1,llm
792	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
793	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
794	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
795	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
796	ENDDO
797	c$OMP END DO NOWAIT
798
799	c$OMP MASTER
800	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
801	c$OMP END MASTER
802	endif ! of if ( .not. pole_nord)
803
804	c$OMP BARRIER
805	c$OMP MASTER
806	call SetDistrib(jj_nb_Physic_bis)
807
808	call VTb(VThallo)
809	c$OMP END MASTER
810	c$OMP BARRIER
811
812	call Register_Hallo(dufi,ip1jmp1,llm,
813	* 1,0,0,1,Request_physic)
814
815	call Register_Hallo(dvfi,ip1jm,llm,
816	* 1,0,0,1,Request_physic)
817
818	call Register_Hallo(dtetafi,ip1jmp1,llm,
819	* 1,0,0,1,Request_physic)
820
821	call Register_Hallo(dpfi,ip1jmp1,1,
822	* 1,0,0,1,Request_physic)
823
824	do j=1,nqtot
825	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
826	* 1,0,0,1,Request_physic)
827	enddo
828
829	call SendRequest(Request_Physic)
830	c$OMP BARRIER
831	call WaitRequest(Request_Physic)
832
833	c$OMP BARRIER
834	c$OMP MASTER
835	call VTe(VThallo)
836
837	call SetDistrib(jj_nb_Physic)
838	c$OMP END MASTER
839	c$OMP BARRIER
840	ijb=ij_begin
841	if (.not. pole_nord) then
842
843	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
844	DO l=1,llm
845	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
846	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
847	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
848	& +dtetafi_tmp(1:iip1,l)
849	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
850	& + dqfi_tmp(1:iip1,l,:)
851	ENDDO
852	c$OMP END DO NOWAIT
853
854	c$OMP MASTER
855	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
856	c$OMP END MASTER
857
858	endif ! of if (.not. pole_nord)
859	c$OMP BARRIER
860	cc$OMP MASTER
861	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
862	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
863	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
864	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
865	cc$OMP END MASTER
866	c
867	c do j=1,nqtot
868	c call WriteField_p('dqfi'//trim(int2str(j)),
869	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
870	c enddo
871
872	c ajout des tendances physiques:
873	c ------------------------------
874	IF (ok_strato) THEN
875	CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
876	ENDIF
877
878	CALL addfi_p( dtphys, leapf, forward ,
879	$ ucov, vcov, teta , q ,ps ,
880	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
881
882	c$OMP BARRIER
883	c$OMP MASTER
884	call VTe(VTphysiq)
885
886	call VTb(VThallo)
887	c$OMP END MASTER
888
889	call SetTag(Request_physic,800)
890	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
891	* jj_Nb_caldyn,Request_physic)
892
893	call Register_SwapField(vcov,vcov,ip1jm,llm,
894	* jj_Nb_caldyn,Request_physic)
895
896	call Register_SwapField(teta,teta,ip1jmp1,llm,
897	* jj_Nb_caldyn,Request_physic)
898
899	call Register_SwapField(masse,masse,ip1jmp1,llm,
900	* jj_Nb_caldyn,Request_physic)
901
902	call Register_SwapField(p,p,ip1jmp1,llmp1,
903	* jj_Nb_caldyn,Request_physic)
904
905	call Register_SwapField(pk,pk,ip1jmp1,llm,
906	* jj_Nb_caldyn,Request_physic)
907
908	call Register_SwapField(phis,phis,ip1jmp1,1,
909	* jj_Nb_caldyn,Request_physic)
910
911	call Register_SwapField(phi,phi,ip1jmp1,llm,
912	* jj_Nb_caldyn,Request_physic)
913
914	call Register_SwapField(w,w,ip1jmp1,llm,
915	* jj_Nb_caldyn,Request_physic)
916
917	do j=1,nqtot
918
919	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
920	* jj_Nb_caldyn,Request_physic)
921
922	enddo
923
924	call SendRequest(Request_Physic)
925	c$OMP BARRIER
926	call WaitRequest(Request_Physic)
927
928	c$OMP BARRIER
929	c$OMP MASTER
930	call VTe(VThallo)
931	call SetDistrib(jj_Nb_caldyn)
932	c$OMP END MASTER
933	c$OMP BARRIER
934	c
935	c Diagnostique de conservation de l'energie : difference
936	IF (ip_ebil_dyn.ge.1 ) THEN
937	ztit='bil phys'
938	CALL diagedyn(ztit,2,1,1,dtphys
939	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
940	ENDIF
941
942	cc$OMP MASTER
943	c if (debug) then
944	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
945	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
946	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
947	c endif
948	cc$OMP END MASTER
949
950
951	c-jld
952	c$OMP MASTER
953	call resume_timer(timer_caldyn)
954	if (FirstPhysic) then
955	ok_start_timer=.TRUE.
956	FirstPhysic=.false.
957	endif
958	c$OMP END MASTER
959	ENDIF ! of IF( apphys )
960
961	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
962	c Calcul academique de la physique = Rappel Newtonien + fritcion
963	c --------------------------------------------------------------
964	cym teta(:,:)=teta(:,:)
965	cym s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
966	ijb=ij_begin
967	ije=ij_end
968	teta(ijb:ije,:)=teta(ijb:ije,:)
969	s -iphysiqdtvr(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
970
971	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
972	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
973	call SendRequest(Request_Physic)
974	c$OMP BARRIER
975	call WaitRequest(Request_Physic)
976
977	call friction_p(ucov,vcov,iphysiq*dtvr)
978	ENDIF ! of IF(iflag_phys.EQ.2)
979
980
981	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
982	c$OMP BARRIER
983	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
984	c$OMP BARRIER
985
986	cc$OMP END PARALLEL
987
988	c-----------------------------------------------------------------------
989	c dissipation horizontale et verticale des petites echelles:
990	c ----------------------------------------------------------
991
992	IF(apdiss) THEN
993	cc$OMP PARALLEL DEFAULT(SHARED)
994	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
995	c$OMP MASTER
996	call suspend_timer(timer_caldyn)
997
998	c print*,'Entree dans la dissipation : Iteration No ',true_itau
999	c calcul de l'energie cinetique avant dissipation
1000	c print *,'Passage dans la dissipation'
1001
1002	call VTb(VThallo)
1003	c$OMP END MASTER
1004
1005	c$OMP BARRIER
1006
1007	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1008	* jj_Nb_dissip,1,1,Request_dissip)
1009
1010	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1011	* jj_Nb_dissip,1,1,Request_dissip)
1012
1013	call Register_SwapField(teta,teta,ip1jmp1,llm,
1014	* jj_Nb_dissip,Request_dissip)
1015
1016	call Register_SwapField(p,p,ip1jmp1,llmp1,
1017	* jj_Nb_dissip,Request_dissip)
1018
1019	call Register_SwapField(pk,pk,ip1jmp1,llm,
1020	* jj_Nb_dissip,Request_dissip)
1021
1022	call SendRequest(Request_dissip)
1023	c$OMP BARRIER
1024	call WaitRequest(Request_dissip)
1025
1026	c$OMP BARRIER
1027	c$OMP MASTER
1028	call SetDistrib(jj_Nb_dissip)
1029	call VTe(VThallo)
1030	call VTb(VTdissipation)
1031	call start_timer(timer_dissip)
1032	c$OMP END MASTER
1033	c$OMP BARRIER
1034
1035	call covcont_p(llm,ucov,vcov,ucont,vcont)
1036	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1037
1038	c dissipation
1039
1040	! CALL FTRACE_REGION_BEGIN("dissip")
1041	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1042	! CALL FTRACE_REGION_END("dissip")
1043
1044	ijb=ij_begin
1045	ije=ij_end
1046	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1047	DO l=1,llm
1048	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1049	ENDDO
1050	c$OMP END DO NOWAIT
1051	if (pole_sud) ije=ije-iip1
1052	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1053	DO l=1,llm
1054	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1055	ENDDO
1056	c$OMP END DO NOWAIT
1057
1058	c teta=teta+dtetadis
1059
1060
1061	c------------------------------------------------------------------------
1062	if (dissip_conservative) then
1063	C On rajoute la tendance due a la transform. Ec -> E therm. cree
1064	C lors de la dissipation
1065	c$OMP BARRIER
1066	c$OMP MASTER
1067	call suspend_timer(timer_dissip)
1068	call VTb(VThallo)
1069	c$OMP END MASTER
1070	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1071	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1072	call SendRequest(Request_Dissip)
1073	c$OMP BARRIER
1074	call WaitRequest(Request_Dissip)
1075	c$OMP MASTER
1076	call VTe(VThallo)
1077	call resume_timer(timer_dissip)
1078	c$OMP END MASTER
1079	c$OMP BARRIER
1080	call covcont_p(llm,ucov,vcov,ucont,vcont)
1081	call enercin_p(vcov,ucov,vcont,ucont,ecin)
1082
1083	ijb=ij_begin
1084	ije=ij_end
1085	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1086	do l=1,llm
1087	do ij=ijb,ije
1088	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1089	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1090	enddo
1091	enddo
1092	c$OMP END DO NOWAIT
1093	endif
1094
1095	ijb=ij_begin
1096	ije=ij_end
1097	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1098	do l=1,llm
1099	do ij=ijb,ije
1100	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1101	enddo
1102	enddo
1103	c$OMP END DO NOWAIT
1104	c------------------------------------------------------------------------
1105
1106
1107	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1108	c ... Calcul de la valeur moyenne, unique de h aux poles .....
1109	c
1110
1111	ijb=ij_begin
1112	ije=ij_end
1113
1114	if (pole_nord) then
1115	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1116	DO l = 1, llm
1117	DO ij = 1,iim
1118	tppn(ij) = aire( ij ) * teta( ij ,l)
1119	ENDDO
1120	tpn = SSUM(iim,tppn,1)/apoln
1121
1122	DO ij = 1, iip1
1123	teta( ij ,l) = tpn
1124	ENDDO
1125	ENDDO
1126	c$OMP END DO NOWAIT
1127
1128	c$OMP MASTER
1129	DO ij = 1,iim
1130	tppn(ij) = aire( ij ) * ps ( ij )
1131	ENDDO
1132	tpn = SSUM(iim,tppn,1)/apoln
1133
1134	DO ij = 1, iip1
1135	ps( ij ) = tpn
1136	ENDDO
1137	c$OMP END MASTER
1138	endif
1139
1140	if (pole_sud) then
1141	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1142	DO l = 1, llm
1143	DO ij = 1,iim
1144	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1145	ENDDO
1146	tps = SSUM(iim,tpps,1)/apols
1147
1148	DO ij = 1, iip1
1149	teta(ij+ip1jm,l) = tps
1150	ENDDO
1151	ENDDO
1152	c$OMP END DO NOWAIT
1153
1154	c$OMP MASTER
1155	DO ij = 1,iim
1156	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1157	ENDDO
1158	tps = SSUM(iim,tpps,1)/apols
1159
1160	DO ij = 1, iip1
1161	ps(ij+ip1jm) = tps
1162	ENDDO
1163	c$OMP END MASTER
1164	endif
1165
1166
1167	c$OMP BARRIER
1168	c$OMP MASTER
1169	call VTe(VTdissipation)
1170
1171	call stop_timer(timer_dissip)
1172
1173	call VTb(VThallo)
1174	c$OMP END MASTER
1175	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1176	* jj_Nb_caldyn,Request_dissip)
1177
1178	call Register_SwapField(vcov,vcov,ip1jm,llm,
1179	* jj_Nb_caldyn,Request_dissip)
1180
1181	call Register_SwapField(teta,teta,ip1jmp1,llm,
1182	* jj_Nb_caldyn,Request_dissip)
1183
1184	call Register_SwapField(p,p,ip1jmp1,llmp1,
1185	* jj_Nb_caldyn,Request_dissip)
1186
1187	call Register_SwapField(pk,pk,ip1jmp1,llm,
1188	* jj_Nb_caldyn,Request_dissip)
1189
1190	call SendRequest(Request_dissip)
1191	c$OMP BARRIER
1192	call WaitRequest(Request_dissip)
1193
1194	c$OMP BARRIER
1195	c$OMP MASTER
1196	call SetDistrib(jj_Nb_caldyn)
1197	call VTe(VThallo)
1198	call resume_timer(timer_caldyn)
1199	c print *,'fin dissipation'
1200	c$OMP END MASTER
1201	c$OMP BARRIER
1202	END IF
1203
1204	cc$OMP END PARALLEL
1205
1206	c ajout debug
1207	c IF( lafin ) then
1208	c abort_message = 'Simulation finished'
1209	c call abort_gcm(modname,abort_message,0)
1210	c ENDIF
1211
1212	c ********************************************************************
1213	c ********************************************************************
1214	c .... fin de l'integration dynamique et physique pour le pas itau ..
1215	c ********************************************************************
1216	c ********************************************************************
1217
1218	c preparation du pas d'integration suivant ......
1219	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1220	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1221	c$OMP MASTER
1222	call stop_timer(timer_caldyn)
1223	c$OMP END MASTER
1224	IF (itau==itaumax) then
1225	c$OMP MASTER
1226	call allgather_timer_average
1227
1228	if (mpi_rank==0) then
1229
1230	print ,'********************************'
1231	print ,'*** TIMER CALDYN ****'
1232	do i=0,mpi_size-1
1233	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1234	& ' : temps moyen :',
1235	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1236	enddo
1237
1238	print ,'********************************'
1239	print ,'*** TIMER VANLEER ****'
1240	do i=0,mpi_size-1
1241	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1242	& ' : temps moyen :',
1243	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1244	enddo
1245
1246	print ,'********************************'
1247	print ,'*** TIMER DISSIP ****'
1248	do i=0,mpi_size-1
1249	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1250	& ' : temps moyen :',
1251	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1252	enddo
1253
1254	print ,'********************************'
1255	print ,'*** TIMER PHYSIC ****'
1256	do i=0,mpi_size-1
1257	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1258	& ' : temps moyen :',
1259	& timer_average(jj_nb_physic(i),timer_physic,i)
1260	enddo
1261
1262	endif
1263
1264	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1265	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1266	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1267	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1268	CALL print_filtre_timer
1269	call fin_getparam
1270	call finalize_parallel
1271	c$OMP END MASTER
1272	c$OMP BARRIER
1273	RETURN
1274	ENDIF
1275
1276	IF ( .NOT.purmats ) THEN
1277	c ........................................................
1278	c .............. schema matsuno + leapfrog ..............
1279	c ........................................................
1280
1281	IF(forward. OR. leapf) THEN
1282	itau= itau + 1
1283	! iday= day_ini+itau/day_step
1284	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1285	! IF(time.GT.1.) THEN
1286	! time = time-1.
1287	! iday = iday+1
1288	! ENDIF
1289	ENDIF
1290
1291
1292	IF( itau. EQ. itaufinp1 ) then
1293
1294	if (flag_verif) then
1295	write(79,*) 'ucov',ucov
1296	write(80,*) 'vcov',vcov
1297	write(81,*) 'teta',teta
1298	write(82,*) 'ps',ps
1299	write(83,*) 'q',q
1300	WRITE(85,*) 'q1 = ',q(:,:,1)
1301	WRITE(86,*) 'q3 = ',q(:,:,3)
1302	endif
1303
1304
1305	c$OMP MASTER
1306	call fin_getparam
1307	call finalize_parallel
1308	c$OMP END MASTER
1309	abort_message = 'Simulation finished'
1310	call abort_gcm(modname,abort_message,0)
1311	RETURN
1312	ENDIF
1313	c-----------------------------------------------------------------------
1314	c ecriture du fichier histoire moyenne:
1315	c -------------------------------------
1316
1317	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1318	c$OMP BARRIER
1319	IF(itau.EQ.itaufin) THEN
1320	iav=1
1321	ELSE
1322	iav=0
1323	ENDIF
1324	#ifdef CPP_IOIPSL
1325	IF (ok_dynzon) THEN
1326	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1327	call SendRequest(TestRequest)
1328	c$OMP BARRIER
1329	call WaitRequest(TestRequest)
1330	c$OMP BARRIER
1331	c$OMP MASTER
1332	! CALL writedynav_p(histaveid, itau,vcov ,
1333	! , ucov,teta,pk,phi,q,masse,ps,phis)
1334
1335	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
1336	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1337	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1338	c$OMP END MASTER
1339	ENDIF !ok_dynzon
1340	#endif
1341	ENDIF
1342
1343	c-----------------------------------------------------------------------
1344	c ecriture de la bande histoire:
1345	c ------------------------------
1346
1347	c IF( MOD(itau,iecri ).EQ.0) THEN
1348
1349	IF( MOD(itau,iecri*day_step).EQ.0) THEN
1350	c$OMP BARRIER
1351	c$OMP MASTER
1352	nbetat = nbetatdem
1353	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1354
1355	cym unat=0.
1356
1357	ijb=ij_begin
1358	ije=ij_end
1359
1360	if (pole_nord) then
1361	ijb=ij_begin+iip1
1362	unat(1:iip1,:)=0.
1363	endif
1364
1365	if (pole_sud) then
1366	ije=ij_end-iip1
1367	unat(ij_end-iip1+1:ij_end,:)=0.
1368	endif
1369
1370	do l=1,llm
1371	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1372	enddo
1373
1374	ijb=ij_begin
1375	ije=ij_end
1376	if (pole_sud) ije=ij_end-iip1
1377
1378	do l=1,llm
1379	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1380	enddo
1381
1382	#ifdef CPP_IOIPSL
1383
1384	! CALL writehist_p(histid,histvid, itau,vcov,
1385	! & ucov,teta,phi,q,masse,ps,phis)
1386
1387	#endif
1388	! For some Grads outputs of fields
1389	if (output_grads_dyn) then
1390	! Ehouarn: hope this works the way I think it does:
1391	call Gather_Field(unat,ip1jmp1,llm,0)
1392	call Gather_Field(vnat,ip1jm,llm,0)
1393	call Gather_Field(teta,ip1jmp1,llm,0)
1394	call Gather_Field(ps,ip1jmp1,1,0)
1395	do iq=1,nqtot
1396	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1397	enddo
1398	if (mpi_rank==0) then
1399	#include "write_grads_dyn.h"
1400	endif
1401	endif ! of if (output_grads_dyn)
1402	c$OMP END MASTER
1403	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1404
1405	IF(itau.EQ.itaufin) THEN
1406
1407	c$OMP BARRIER
1408	c$OMP MASTER
1409
1410	if (planet_type.eq."earth") then
1411	! Write an Earth-format restart file
1412	CALL dynredem1_p("restart.nc",0.0,
1413	& vcov,ucov,teta,q,masse,ps)
1414	endif ! of if (planet_type.eq."earth")
1415
1416	! CLOSE(99)
1417	c$OMP END MASTER
1418	ENDIF ! of IF (itau.EQ.itaufin)
1419
1420	c-----------------------------------------------------------------------
1421	c gestion de l'integration temporelle:
1422	c ------------------------------------
1423
1424	IF( MOD(itau,iperiod).EQ.0 ) THEN
1425	GO TO 1
1426	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1427
1428	IF( forward ) THEN
1429	c fin du pas forward et debut du pas backward
1430
1431	forward = .FALSE.
1432	leapf = .FALSE.
1433	GO TO 2
1434
1435	ELSE
1436	c fin du pas backward et debut du premier pas leapfrog
1437
1438	leapf = .TRUE.
1439	dt = 2.*dtvr
1440	GO TO 2
1441	END IF
1442	ELSE
1443
1444	c ...... pas leapfrog .....
1445
1446	leapf = .TRUE.
1447	dt = 2.*dtvr
1448	GO TO 2
1449	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1450	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1451
1452
1453	ELSE ! of IF (.not.purmats)
1454
1455	c ........................................................
1456	c .............. schema matsuno ...............
1457	c ........................................................
1458	IF( forward ) THEN
1459
1460	itau = itau + 1
1461	! iday = day_ini+itau/day_step
1462	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1463	!
1464	! IF(time.GT.1.) THEN
1465	! time = time-1.
1466	! iday = iday+1
1467	! ENDIF
1468
1469	forward = .FALSE.
1470	IF( itau. EQ. itaufinp1 ) then
1471	c$OMP MASTER
1472	call fin_getparam
1473	call finalize_parallel
1474	c$OMP END MASTER
1475	abort_message = 'Simulation finished'
1476	call abort_gcm(modname,abort_message,0)
1477	RETURN
1478	ENDIF
1479	GO TO 2
1480
1481	ELSE ! of IF(forward)
1482
1483	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1484	IF(itau.EQ.itaufin) THEN
1485	iav=1
1486	ELSE
1487	iav=0
1488	ENDIF
1489	#ifdef CPP_IOIPSL
1490	IF (ok_dynzon) THEN
1491	c$OMP BARRIER
1492
1493	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1494	call SendRequest(TestRequest)
1495	c$OMP BARRIER
1496	call WaitRequest(TestRequest)
1497
1498	c$OMP BARRIER
1499	c$OMP MASTER
1500	! CALL writedynav_p(histaveid, itau,vcov ,
1501	! , ucov,teta,pk,phi,q,masse,ps,phis)
1502	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1503	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1504	c$OMP END MASTER
1505	END IF !ok_dynzon
1506	#endif
1507	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1508
1509
1510	c IF(MOD(itau,iecri ).EQ.0) THEN
1511	IF(MOD(itau,iecri*day_step).EQ.0) THEN
1512	c$OMP BARRIER
1513	c$OMP MASTER
1514	nbetat = nbetatdem
1515	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1516
1517	cym unat=0.
1518	ijb=ij_begin
1519	ije=ij_end
1520
1521	if (pole_nord) then
1522	ijb=ij_begin+iip1
1523	unat(1:iip1,:)=0.
1524	endif
1525
1526	if (pole_sud) then
1527	ije=ij_end-iip1
1528	unat(ij_end-iip1+1:ij_end,:)=0.
1529	endif
1530
1531	do l=1,llm
1532	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1533	enddo
1534
1535	ijb=ij_begin
1536	ije=ij_end
1537	if (pole_sud) ije=ij_end-iip1
1538
1539	do l=1,llm
1540	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1541	enddo
1542
1543	#ifdef CPP_IOIPSL
1544
1545	! CALL writehist_p(histid, histvid, itau,vcov ,
1546	! & ucov,teta,phi,q,masse,ps,phis)
1547	#endif
1548	! For some Grads output (but does it work?)
1549	if (output_grads_dyn) then
1550	call Gather_Field(unat,ip1jmp1,llm,0)
1551	call Gather_Field(vnat,ip1jm,llm,0)
1552	call Gather_Field(teta,ip1jmp1,llm,0)
1553	call Gather_Field(ps,ip1jmp1,1,0)
1554	do iq=1,nqtot
1555	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1556	enddo
1557	c
1558	if (mpi_rank==0) then
1559	#include "write_grads_dyn.h"
1560	endif
1561	endif ! of if (output_grads_dyn)
1562
1563	c$OMP END MASTER
1564	ENDIF ! of IF(MOD(itau,iecri*day_step).EQ.0)
1565
1566	IF(itau.EQ.itaufin) THEN
1567	if (planet_type.eq."earth") then
1568	c$OMP MASTER
1569	CALL dynredem1_p("restart.nc",0.0,
1570	. vcov,ucov,teta,q,masse,ps)
1571	c$OMP END MASTER
1572	endif ! of if (planet_type.eq."earth")
1573	ENDIF ! of IF(itau.EQ.itaufin)
1574
1575	forward = .TRUE.
1576	GO TO 1
1577
1578	ENDIF ! of IF (forward)
1579
1580	END IF ! of IF(.not.purmats)
1581	c$OMP MASTER
1582	call fin_getparam
1583	call finalize_parallel
1584	c$OMP END MASTER
1585	RETURN
1586	END

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

Download in other formats:

Original Format