Context Navigation

leapfrog_p.F @ 4400

Last change on this file since 4400 was 1676, checked in by aslmd, 12 years ago

Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation

Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon

Author : AS equipe planeto

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

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

Download in other formats:

Original Format