Context Navigation

leapfrog_p.F @ 1999

Last change on this file since 1999 was 1999, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1920:1997 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.5 KB

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

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

Download in other formats:

Original Format