Context Navigation

leapfrog_loc.F @ 2435

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

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

Download in other formats:

Original Format