Context Navigation

leapfrog_p.F @ 1616

Last change on this file since 1616 was 1616, checked in by Ehouarn Millour, 13 years ago
Some cleanup around what is done during the integration step of dynamical tendencies and namely removed computation of (unused) finvmaold, thereby saving us the expense of a call to the (costly) filter at every dynamical time step. Checked (on Vargas, in seq, omp, mpi and mixed mode) that this doesn't change the GCM results, as expected. EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.7 KB

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

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

Download in other formats:

Original Format