Context Navigation

leapfrog_loc.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.0 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 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	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
673
674	CALL caladvtrac_loc(q,pbaru,pbarv,
675	* p, masse, dq, teta,
676	. flxw,pk, iapptrac)
677
678	! do j=1,nqtot
679	! call WriteField_u('qadv'//trim(int2str(j)),q(:,:,j))
680	! enddo
681
682	c
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
1143	cc$OMP END PARALLEL
1144
1145	c-----------------------------------------------------------------------
1146	c dissipation horizontale et verticale des petites echelles:
1147	c ----------------------------------------------------------
1148
1149	IF(apdiss) THEN
1150
1151	CALL call_dissip(ucov,vcov,teta,p,pk,ps)
1152	!cc$OMP PARALLEL DEFAULT(SHARED)
1153	!cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1154	!c$OMP MASTER
1155	! call suspend_timer(timer_caldyn)
1156	!
1157	!c print*,'Entree dans la dissipation : Iteration No ',true_itau
1158	!c calcul de l'energie cinetique avant dissipation
1159	!c print *,'Passage dans la dissipation'
1160
1161	! call VTb(VThallo)
1162	!c$OMP END MASTER
1163
1164	!c$OMP BARRIER
1165
1166	! call Register_SwapField_u(ucov,ucov,distrib_dissip,
1167	! * Request_dissip,up=1,down=1)
1168
1169	! call Register_SwapField_v(vcov,vcov,distrib_dissip,
1170	! * Request_dissip,up=1,down=1)
1171
1172	! call Register_SwapField_u(teta,teta,distrib_dissip,
1173	! * Request_dissip)
1174
1175	! call Register_SwapField_u(p,p,distrib_dissip,
1176	! * Request_dissip)
1177
1178	! call Register_SwapField_u(pk,pk,distrib_dissip,
1179	! * Request_dissip)
1180
1181	! call SendRequest(Request_dissip)
1182	!c$OMP BARRIER
1183	! call WaitRequest(Request_dissip)
1184
1185	!c$OMP BARRIER
1186	!c$OMP MASTER
1187	! call set_distrib(distrib_dissip)
1188	! call VTe(VThallo)
1189	! call VTb(VTdissipation)
1190	! call start_timer(timer_dissip)
1191	!c$OMP END MASTER
1192	!c$OMP BARRIER
1193
1194	! call covcont_loc(llm,ucov,vcov,ucont,vcont)
1195	! call enercin_loc(vcov,ucov,vcont,ucont,ecin0)
1196
1197	!c dissipation
1198
1199	!! CALL FTRACE_REGION_BEGIN("dissip")
1200	! CALL dissip_loc(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1201
1202	!#ifdef DEBUG_IO
1203	! call WriteField_u('dudis',dudis)
1204	! call WriteField_v('dvdis',dvdis)
1205	! call WriteField_u('dtetadis',dtetadis)
1206	!#endif
1207	!
1208	!! CALL FTRACE_REGION_END("dissip")
1209	!
1210	! ijb=ij_begin
1211	! ije=ij_end
1212	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1213	! DO l=1,llm
1214	! ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1215	! ENDDO
1216	!c$OMP END DO NOWAIT
1217	! if (pole_sud) ije=ije-iip1
1218	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1219	! DO l=1,llm
1220	! vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1221	! ENDDO
1222	!c$OMP END DO NOWAIT
1223
1224	!c teta=teta+dtetadis
1225
1226
1227	!c------------------------------------------------------------------------
1228	! if (dissip_conservative) then
1229	!C On rajoute la tendance due a la transform. Ec -> E therm. cree
1230	!C lors de la dissipation
1231	!c$OMP BARRIER
1232	!c$OMP MASTER
1233	! call suspend_timer(timer_dissip)
1234	! call VTb(VThallo)
1235	!c$OMP END MASTER
1236	! call Register_Hallo_u(ucov,llm,1,1,1,1,Request_Dissip)
1237	! call Register_Hallo_v(vcov,llm,1,1,1,1,Request_Dissip)
1238	! call SendRequest(Request_Dissip)
1239	!c$OMP BARRIER
1240	! call WaitRequest(Request_Dissip)
1241	!c$OMP MASTER
1242	! call VTe(VThallo)
1243	! call resume_timer(timer_dissip)
1244	!c$OMP END MASTER
1245	!c$OMP BARRIER
1246	! call covcont_loc(llm,ucov,vcov,ucont,vcont)
1247	! call enercin_loc(vcov,ucov,vcont,ucont,ecin)
1248	!
1249	! ijb=ij_begin
1250	! ije=ij_end
1251	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1252	! do l=1,llm
1253	! do ij=ijb,ije
1254	! dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1255	! dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1256	! enddo
1257	! enddo
1258	!c$OMP END DO NOWAIT
1259	! endif
1260
1261	! ijb=ij_begin
1262	! ije=ij_end
1263	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1264	! do l=1,llm
1265	! do ij=ijb,ije
1266	! teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1267	! enddo
1268	! enddo
1269	!c$OMP END DO NOWAIT
1270	!c------------------------------------------------------------------------
1271
1272
1273	!c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1274	!c ... Calcul de la valeur moyenne, unique de h aux poles .....
1275	!c
1276
1277	! ijb=ij_begin
1278	! ije=ij_end
1279	!
1280	! if (pole_nord) then
1281	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1282	! DO l = 1, llm
1283	! DO ij = 1,iim
1284	! tppn(ij) = aire( ij ) * teta( ij ,l)
1285	! ENDDO
1286	! tpn = SSUM(iim,tppn,1)/apoln
1287
1288	! DO ij = 1, iip1
1289	! teta( ij ,l) = tpn
1290	! ENDDO
1291	! ENDDO
1292	!c$OMP END DO NOWAIT
1293
1294	!c$OMP MASTER
1295	! DO ij = 1,iim
1296	! tppn(ij) = aire( ij ) * ps ( ij )
1297	! ENDDO
1298	! tpn = SSUM(iim,tppn,1)/apoln
1299	!
1300	! DO ij = 1, iip1
1301	! ps( ij ) = tpn
1302	! ENDDO
1303	!c$OMP END MASTER
1304	! endif
1305	!
1306	! if (pole_sud) then
1307	!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1308	! DO l = 1, llm
1309	! DO ij = 1,iim
1310	! tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1311	! ENDDO
1312	! tps = SSUM(iim,tpps,1)/apols
1313
1314	! DO ij = 1, iip1
1315	! teta(ij+ip1jm,l) = tps
1316	! ENDDO
1317	! ENDDO
1318	!c$OMP END DO NOWAIT
1319
1320	!c$OMP MASTER
1321	! DO ij = 1,iim
1322	! tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1323	! ENDDO
1324	! tps = SSUM(iim,tpps,1)/apols
1325	!
1326	! DO ij = 1, iip1
1327	! ps(ij+ip1jm) = tps
1328	! ENDDO
1329	!c$OMP END MASTER
1330	! endif
1331
1332
1333	!c$OMP BARRIER
1334	!c$OMP MASTER
1335	! call VTe(VTdissipation)
1336
1337	! call stop_timer(timer_dissip)
1338	!
1339	! call VTb(VThallo)
1340	!c$OMP END MASTER
1341	! call Register_SwapField_u(ucov,ucov,distrib_caldyn,
1342	! * Request_dissip)
1343
1344	! call Register_SwapField_v(vcov,vcov,distrib_caldyn,
1345	! * Request_dissip)
1346
1347	! call Register_SwapField_u(teta,teta,distrib_caldyn,
1348	! * Request_dissip)
1349
1350	! call Register_SwapField_u(p,p,distrib_caldyn,
1351	! * Request_dissip)
1352
1353	! call Register_SwapField_u(pk,pk,distrib_caldyn,
1354	! * Request_dissip)
1355
1356	! call SendRequest(Request_dissip)
1357	!c$OMP BARRIER
1358	! call WaitRequest(Request_dissip)
1359
1360	!c$OMP BARRIER
1361	!c$OMP MASTER
1362	! call set_distrib(distrib_caldyn)
1363	! call VTe(VThallo)
1364	! call resume_timer(timer_caldyn)
1365	!c print *,'fin dissipation'
1366	!c$OMP END MASTER
1367	!c$OMP BARRIER
1368	END IF ! of IF(apdiss)
1369
1370	cc$OMP END PARALLEL
1371
1372	c ajout debug
1373	c IF( lafin ) then
1374	c abort_message = 'Simulation finished'
1375	c call abort_gcm(modname,abort_message,0)
1376	c ENDIF
1377
1378	c ********************************************************************
1379	c ********************************************************************
1380	c .... fin de l'integration dynamique et physique pour le pas itau ..
1381	c ********************************************************************
1382	c ********************************************************************
1383
1384	c preparation du pas d'integration suivant ......
1385	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1386	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1387	c$OMP MASTER
1388	call stop_timer(timer_caldyn)
1389	c$OMP END MASTER
1390	IF (itau==itaumax) then
1391	c$OMP MASTER
1392	call allgather_timer_average
1393	call barrier
1394	if (mpi_rank==0) then
1395
1396	print ,'********************************'
1397	print ,'*** TIMER CALDYN ****'
1398	do i=0,mpi_size-1
1399	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1400	& ' : temps moyen :',
1401	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1402	enddo
1403
1404	print ,'********************************'
1405	print ,'*** TIMER VANLEER ****'
1406	do i=0,mpi_size-1
1407	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1408	& ' : temps moyen :',
1409	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1410	enddo
1411
1412	print ,'********************************'
1413	print ,'*** TIMER DISSIP ****'
1414	do i=0,mpi_size-1
1415	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1416	& ' : temps moyen :',
1417	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1418	enddo
1419
1420	print ,'********************************'
1421	print ,'*** TIMER PHYSIC ****'
1422	do i=0,mpi_size-1
1423	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1424	& ' : temps moyen :',
1425	& timer_average(jj_nb_physic(i),timer_physic,i)
1426	enddo
1427
1428	endif
1429	CALL barrier
1430	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1431	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1432	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1433	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1434	CALL print_filtre_timer
1435	c$OMP END MASTER
1436	CALL dynredem1_loc("restart.nc",0.0,
1437	. vcov,ucov,teta,q,masse,ps)
1438	c$OMP MASTER
1439	call fin_getparam
1440	call finalize_parallel
1441	c$OMP END MASTER
1442	c$OMP BARRIER
1443	RETURN
1444	ENDIF
1445
1446	IF ( .NOT.purmats ) THEN
1447	c ........................................................
1448	c .............. schema matsuno + leapfrog ..............
1449	c ........................................................
1450
1451	IF(forward. OR. leapf) THEN
1452	itau= itau + 1
1453	! iday= day_ini+itau/day_step
1454	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1455	! IF(time.GT.1.) THEN
1456	! time = time-1.
1457	! iday = iday+1
1458	! ENDIF
1459	ENDIF
1460
1461
1462	IF( itau. EQ. itaufinp1 ) then
1463
1464	if (flag_verif) then
1465	write(79,*) 'ucov',ucov
1466	write(80,*) 'vcov',vcov
1467	write(81,*) 'teta',teta
1468	write(82,*) 'ps',ps
1469	write(83,*) 'q',q
1470	WRITE(85,*) 'q1 = ',q(:,:,1)
1471	WRITE(86,*) 'q3 = ',q(:,:,3)
1472	endif
1473
1474
1475	c$OMP MASTER
1476	call fin_getparam
1477	call finalize_parallel
1478	c$OMP END MASTER
1479	abort_message = 'Simulation finished'
1480	call abort_gcm(modname,abort_message,0)
1481	RETURN
1482	ENDIF
1483	c-----------------------------------------------------------------------
1484	c ecriture du fichier histoire moyenne:
1485	c -------------------------------------
1486
1487	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1488	c$OMP BARRIER
1489	IF(itau.EQ.itaufin) THEN
1490	iav=1
1491	ELSE
1492	iav=0
1493	ENDIF
1494
1495	#ifdef CPP_IOIPSL
1496	IF (ok_dynzon) THEN
1497
1498	CALL bilan_dyn_loc(2,dtvriperiod,dtvrday_step*periodav,
1499	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1500
1501	ENDIF !ok_dynzon
1502
1503	IF (ok_dyn_ave) THEN
1504	CALL writedynav_loc(itau,vcov,
1505	& ucov,teta,pk,phi,q,masse,ps,phis)
1506	ENDIF
1507	#endif
1508	ENDIF
1509
1510	c-----------------------------------------------------------------------
1511	c ecriture de la bande histoire:
1512	c ------------------------------
1513
1514	IF( MOD(itau,iecri).EQ.0) THEN
1515	! Ehouarn: output only during LF or Backward Matsuno
1516	if (leapf.or.(.not.leapf.and.(.not.forward))) then
1517
1518	c$OMP BARRIER
1519	c$OMP MASTER
1520	CALL geopot_loc(ip1jmp1,teta,pk,pks,phis,phi)
1521	c$OMP END MASTER
1522	c$OMP BARRIER
1523
1524	#ifdef CPP_IOIPSL
1525	if (ok_dyn_ins) then
1526	CALL writehist_loc(itau,vcov,ucov,teta,phi,q,
1527	& masse,ps,phis)
1528	endif
1529	#endif
1530	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1531	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1532
1533	IF(itau.EQ.itaufin) THEN
1534
1535	c$OMP BARRIER
1536
1537	! if (planet_type.eq."earth") then
1538	! Write an Earth-format restart file
1539	CALL dynredem1_loc("restart.nc",0.0,
1540	& vcov,ucov,teta,q,masse,ps)
1541	! endif ! of if (planet_type.eq."earth")
1542
1543	! CLOSE(99)
1544	ENDIF ! of IF (itau.EQ.itaufin)
1545
1546	c-----------------------------------------------------------------------
1547	c gestion de l'integration temporelle:
1548	c ------------------------------------
1549
1550	IF( MOD(itau,iperiod).EQ.0 ) THEN
1551	GO TO 1
1552	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1553
1554	IF( forward ) THEN
1555	c fin du pas forward et debut du pas backward
1556
1557	forward = .FALSE.
1558	leapf = .FALSE.
1559	GO TO 2
1560
1561	ELSE
1562	c fin du pas backward et debut du premier pas leapfrog
1563
1564	leapf = .TRUE.
1565	dt = 2.*dtvr
1566	GO TO 2
1567	END IF
1568	ELSE
1569
1570	c ...... pas leapfrog .....
1571
1572	leapf = .TRUE.
1573	dt = 2.*dtvr
1574	GO TO 2
1575	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1576	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1577
1578
1579	ELSE ! of IF (.not.purmats)
1580
1581	c ........................................................
1582	c .............. schema matsuno ...............
1583	c ........................................................
1584	IF( forward ) THEN
1585
1586	itau = itau + 1
1587	! iday = day_ini+itau/day_step
1588	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1589	!
1590	! IF(time.GT.1.) THEN
1591	! time = time-1.
1592	! iday = iday+1
1593	! ENDIF
1594
1595	forward = .FALSE.
1596	IF( itau. EQ. itaufinp1 ) then
1597	c$OMP MASTER
1598	call fin_getparam
1599	call finalize_parallel
1600	c$OMP END MASTER
1601	abort_message = 'Simulation finished'
1602	call abort_gcm(modname,abort_message,0)
1603	RETURN
1604	ENDIF
1605	GO TO 2
1606
1607	ELSE ! of IF(forward) i.e. backward step
1608
1609	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1610	IF(itau.EQ.itaufin) THEN
1611	iav=1
1612	ELSE
1613	iav=0
1614	ENDIF
1615
1616	#ifdef CPP_IOIPSL
1617	IF (ok_dynzon) THEN
1618	CALL bilan_dyn_loc(2,dtvriperiod,dtvrday_step*periodav,
1619	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1620	ENDIF
1621
1622	IF (ok_dyn_ave) THEN
1623	CALL writedynav_loc(itau,vcov,
1624	& ucov,teta,pk,phi,q,masse,ps,phis)
1625	ENDIF
1626	#endif
1627	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1628
1629
1630	IF(MOD(itau,iecri ).EQ.0) THEN
1631
1632	c$OMP BARRIER
1633	c$OMP MASTER
1634	CALL geopot_loc(ip1jmp1,teta,pk,pks,phis,phi)
1635	c$OMP END MASTER
1636	c$OMP BARRIER
1637
1638
1639	#ifdef CPP_IOIPSL
1640	if (ok_dyn_ins) then
1641	CALL writehist_loc(itau,vcov,ucov,teta,phi,q,
1642	& masse,ps,phis)
1643	endif ! of if (ok_dyn_ins)
1644	#endif
1645	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1646
1647
1648	IF(itau.EQ.itaufin) THEN
1649	! if (planet_type.eq."earth") then
1650	CALL dynredem1_loc("restart.nc",0.0,
1651	. vcov,ucov,teta,q,masse,ps)
1652	! endif ! of if (planet_type.eq."earth")
1653	ENDIF ! of IF(itau.EQ.itaufin)
1654
1655	forward = .TRUE.
1656	GO TO 1
1657
1658	ENDIF ! of IF (forward)
1659
1660	END IF ! of IF(.not.purmats)
1661	c$OMP MASTER
1662	call fin_getparam
1663	call finalize_parallel
1664	c$OMP END MASTER
1665	abort_message = 'Simulation finished'
1666	call abort_gcm(modname,abort_message,0)
1667	RETURN
1668	END

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

Download in other formats:

Original Format