Context Navigation

leapfrog_loc.F @ 1997

Last change on this file since 1997 was 1987, checked in by Ehouarn Millour, 11 years ago
Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency). Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid). EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 51.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