Context Navigation

leapfrog_loc.F @ 2242

Last change on this file since 2242 was 2221, checked in by Ehouarn Millour, 9 years ago
Some cleanup: remove (unused) clesph0 from dynamics. EM
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
File size: 51.6 KB

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

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

Download in other formats:

Original Format