Context Navigation

leapfrog_loc.F @ 1999

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

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

Download in other formats:

Original Format