Context Navigation

leapfrog_loc.F @ 5360

Last change on this file since 5360 was 2375, checked in by Ehouarn Millour, 9 years ago
Fix in the computation of the date; the convention is that it corresponds to the time at the end of the current dynamics or physics step (except when in backward Matsuno step where it remains unchanged as it was correctly updated during the forward part of the step). Note that the relationship between itau and date is a bit tricky as itau is incremented between Matsuno forward and backward steps (and from a leapfrog step to the next) but unchanged from Matsuno bacward step to leapfrog step. Because of change in date when calling physics, bench results will change with this revision. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 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