Context Navigation

leapfrog_p.F @ 1793

Last change on this file since 1793 was 1793, checked in by Ehouarn Millour, 11 years ago

Improved sponge layer:

Sponge tendencies are now computed analytically, instead than using a Forward Euler approximation.
Sponge tendencies are added within top_bound, and the sponge is applied after physics tendencies have been taken into account.

These changes imply that GCM results (when using sponge layer) will be differentwrt bench test cases using previous revisions.
EM

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

File size: 55.0 KB

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

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

Download in other formats:

Original Format