Context Navigation

leapfrog_loc.F @ 2339

Last change on this file since 2339 was 2286, checked in by Ehouarn Millour, 10 years ago

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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

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

Download in other formats:

Original Format