Context Navigation

leapfrog_loc.F

Last change on this file was 1705, checked in by Ehouarn Millour, 12 years ago

Added arch files for ADA (IDRIS IBMx3750) and made the following code modifications:
phylmd/printflag.F : removed "print" of unset variable (radpas0)
dyn3dmem/integrd_loc.F : removed unecessary "include mpif.h"
dyn3dmem/leapfrog_loc.F : removed unecessary "include mpif.h" and allocate saved variables at first call
dyn3dmem/mod_filtreg_p.F : added matmul() alternatives to call to BLAS routine SGEMM (which was incorectly set as DGEMM; which would fail if running with -r4)
filtrez/filtreg.F: changed calls to DGEMM into calls to SGEMM, so that code works with either -r4 or -r8 (the later being used in conjunction with "BLAS SGEMV=DGEMV SGEMM=DGEMM" preprocessing statements)
EM

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

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

Download in other formats:

Original Format