Context Navigation

leapfrog_p.F @ 1671

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

fixed "aquaplanet case" so that initializations (creation of files startphy.nc and limit.nc) now also works in parallel (mpi,omp,mixed).
call to "iniaqua" is now done from within "iniphysiq" ; also added some tests to check consistency between essential variables shared by dynamics and physics (planetary radius, gravity, Cp, ...)
simillarily adapted "phydev" routines, and added necessary routines to also be able to generate startphy.nc/restartphy.nc files there. Also removed common file "comcstphy.h" and replaced it with a module "comcstphy.F90"

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 55.6 KB

Line
1	!
2	! $Id: leapfrog_p.F 1671 2012-10-24 07:10:10Z 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 :: 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	jH_cur = jH_ref + start_time + &
713	& mod(itau,day_step)/float(day_step)
714	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
715	if (jH_cur > 1.0 ) then
716	jD_cur = jD_cur +1.
717	jH_cur = jH_cur -1.
718	endif
719
720	c rajout debug
721	c lafin = .true.
722
723
724	c Inbterface avec les routines de phylmd (phymars ... )
725	c -----------------------------------------------------
726
727	c+jld
728
729	c Diagnostique de conservation de l'energie : initialisation
730	IF (ip_ebil_dyn.ge.1 ) THEN
731	ztit='bil dyn'
732	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
733	IF (planet_type.eq."earth") THEN
734	#ifdef CPP_EARTH
735	CALL diagedyn(ztit,2,1,1,dtphys
736	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
737	#endif
738	ENDIF
739	ENDIF
740	c-jld
741	c$OMP BARRIER
742	c$OMP MASTER
743	call VTb(VThallo)
744	c$OMP END MASTER
745
746	call SetTag(Request_physic,800)
747
748	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
749	* jj_Nb_physic,2,2,Request_physic)
750
751	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
752	* jj_Nb_physic,2,2,Request_physic)
753
754	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
755	* jj_Nb_physic,2,2,Request_physic)
756
757	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
758	* jj_Nb_physic,1,2,Request_physic)
759
760	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
761	* jj_Nb_physic,2,2,Request_physic)
762
763	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
764	* jj_Nb_physic,2,2,Request_physic)
765
766	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
767	* jj_Nb_physic,2,2,Request_physic)
768
769	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
770	* jj_Nb_physic,2,2,Request_physic)
771
772	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
773	* jj_Nb_physic,2,2,Request_physic)
774
775	c call SetDistrib(jj_nb_vanleer)
776	do j=1,nqtot
777
778	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
779	* jj_Nb_physic,2,2,Request_physic)
780	enddo
781
782	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
783	* jj_Nb_physic,2,2,Request_physic)
784
785	call SendRequest(Request_Physic)
786	c$OMP BARRIER
787	call WaitRequest(Request_Physic)
788
789	c$OMP BARRIER
790	c$OMP MASTER
791	call SetDistrib(jj_nb_Physic)
792	call VTe(VThallo)
793
794	call VTb(VTphysiq)
795	c$OMP END MASTER
796	c$OMP BARRIER
797
798	cc$OMP MASTER
799	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
800	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
801	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
802	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
803	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
804	cc$OMP END MASTER
805	cc$OMP BARRIER
806	! CALL FTRACE_REGION_BEGIN("calfis")
807	CALL calfis_p(lafin ,jD_cur, jH_cur,
808	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
809	$ du,dv,dteta,dq,
810	$ flxw,
811	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
812	! CALL FTRACE_REGION_END("calfis")
813	ijb=ij_begin
814	ije=ij_end
815	if ( .not. pole_nord) then
816	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
817	DO l=1,llm
818	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
819	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
820	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
821	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
822	ENDDO
823	c$OMP END DO NOWAIT
824
825	c$OMP MASTER
826	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
827	c$OMP END MASTER
828	endif ! of if ( .not. pole_nord)
829
830	c$OMP BARRIER
831	c$OMP MASTER
832	call SetDistrib(jj_nb_Physic_bis)
833
834	call VTb(VThallo)
835	c$OMP END MASTER
836	c$OMP BARRIER
837
838	call Register_Hallo(dufi,ip1jmp1,llm,
839	* 1,0,0,1,Request_physic)
840
841	call Register_Hallo(dvfi,ip1jm,llm,
842	* 1,0,0,1,Request_physic)
843
844	call Register_Hallo(dtetafi,ip1jmp1,llm,
845	* 1,0,0,1,Request_physic)
846
847	call Register_Hallo(dpfi,ip1jmp1,1,
848	* 1,0,0,1,Request_physic)
849
850	do j=1,nqtot
851	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
852	* 1,0,0,1,Request_physic)
853	enddo
854
855	call SendRequest(Request_Physic)
856	c$OMP BARRIER
857	call WaitRequest(Request_Physic)
858
859	c$OMP BARRIER
860	c$OMP MASTER
861	call VTe(VThallo)
862
863	call SetDistrib(jj_nb_Physic)
864	c$OMP END MASTER
865	c$OMP BARRIER
866	ijb=ij_begin
867	if (.not. pole_nord) then
868
869	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
870	DO l=1,llm
871	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
872	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
873	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
874	& +dtetafi_tmp(1:iip1,l)
875	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
876	& + dqfi_tmp(1:iip1,l,:)
877	ENDDO
878	c$OMP END DO NOWAIT
879
880	c$OMP MASTER
881	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
882	c$OMP END MASTER
883
884	endif ! of if (.not. pole_nord)
885	c$OMP BARRIER
886	cc$OMP MASTER
887	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
888	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
889	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
890	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
891	cc$OMP END MASTER
892	c
893	c do j=1,nqtot
894	c call WriteField_p('dqfi'//trim(int2str(j)),
895	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
896	c enddo
897
898	c ajout des tendances physiques:
899	c ------------------------------
900	IF (ok_strato) THEN
901	CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
902	ENDIF
903
904	CALL addfi_p( dtphys, leapf, forward ,
905	$ ucov, vcov, teta , q ,ps ,
906	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
907
908	c$OMP BARRIER
909	c$OMP MASTER
910	call VTe(VTphysiq)
911
912	call VTb(VThallo)
913	c$OMP END MASTER
914
915	call SetTag(Request_physic,800)
916	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
917	* jj_Nb_caldyn,Request_physic)
918
919	call Register_SwapField(vcov,vcov,ip1jm,llm,
920	* jj_Nb_caldyn,Request_physic)
921
922	call Register_SwapField(teta,teta,ip1jmp1,llm,
923	* jj_Nb_caldyn,Request_physic)
924
925	call Register_SwapField(masse,masse,ip1jmp1,llm,
926	* jj_Nb_caldyn,Request_physic)
927
928	call Register_SwapField(p,p,ip1jmp1,llmp1,
929	* jj_Nb_caldyn,Request_physic)
930
931	call Register_SwapField(pk,pk,ip1jmp1,llm,
932	* jj_Nb_caldyn,Request_physic)
933
934	call Register_SwapField(phis,phis,ip1jmp1,1,
935	* jj_Nb_caldyn,Request_physic)
936
937	call Register_SwapField(phi,phi,ip1jmp1,llm,
938	* jj_Nb_caldyn,Request_physic)
939
940	call Register_SwapField(w,w,ip1jmp1,llm,
941	* jj_Nb_caldyn,Request_physic)
942
943	do j=1,nqtot
944
945	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
946	* jj_Nb_caldyn,Request_physic)
947
948	enddo
949
950	call SendRequest(Request_Physic)
951	c$OMP BARRIER
952	call WaitRequest(Request_Physic)
953
954	c$OMP BARRIER
955	c$OMP MASTER
956	call VTe(VThallo)
957	call SetDistrib(jj_Nb_caldyn)
958	c$OMP END MASTER
959	c$OMP BARRIER
960	c
961	c Diagnostique de conservation de l'energie : difference
962	IF (ip_ebil_dyn.ge.1 ) THEN
963	ztit='bil phys'
964	CALL diagedyn(ztit,2,1,1,dtphys
965	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
966	ENDIF
967
968	cc$OMP MASTER
969	c if (debug) then
970	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
971	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
972	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
973	c endif
974	cc$OMP END MASTER
975
976
977	c-jld
978	c$OMP MASTER
979	call resume_timer(timer_caldyn)
980	if (FirstPhysic) then
981	ok_start_timer=.TRUE.
982	FirstPhysic=.false.
983	endif
984	c$OMP END MASTER
985	ENDIF ! of IF( apphys )
986
987	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
988	! Academic case : Simple friction and Newtonan relaxation
989	! -------------------------------------------------------
990	ijb=ij_begin
991	ije=ij_end
992	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
993	do l=1,llm
994	teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
995	& (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
996	& (knewt_g+knewt_t(l)*clat4(ijb:ije))
997	enddo ! of do l=1,llm
998	!$OMP END DO
999
1000	!$OMP MASTER
1001	if (planet_type.eq."giant") then
1002	! add an intrinsic heat flux at the base of the atmosphere
1003	teta(ijb:ije,1) = teta(ijb:ije,1)
1004	& + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
1005	endif
1006	!$OMP END MASTER
1007	!$OMP BARRIER
1008
1009	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
1010	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
1011	call SendRequest(Request_Physic)
1012	c$OMP BARRIER
1013	call WaitRequest(Request_Physic)
1014	c$OMP BARRIER
1015	call friction_p(ucov,vcov,dtvr)
1016	!$OMP BARRIER
1017
1018	! Sponge layer (if any)
1019	IF (ok_strato) THEN
1020	! set dufi,dvfi,... to zero
1021	ijb=ij_begin
1022	ije=ij_end
1023	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1024	do l=1,llm
1025	dufi(ijb:ije,l)=0
1026	dtetafi(ijb:ije,l)=0
1027	dqfi(ijb:ije,l,1:nqtot)=0
1028	enddo
1029	!$OMP END DO
1030	!$OMP MASTER
1031	dpfi(ijb:ije)=0
1032	!$OMP END MASTER
1033	ijb=ij_begin
1034	ije=ij_end
1035	if (pole_sud) ije=ije-iip1
1036	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1037	do l=1,llm
1038	dvfi(ijb:ije,l)=0
1039	enddo
1040	!$OMP END DO
1041
1042	CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
1043	CALL addfi_p( dtvr, leapf, forward ,
1044	$ ucov, vcov, teta , q ,ps ,
1045	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
1046	!$OMP BARRIER
1047	ENDIF ! of IF (ok_strato)
1048	ENDIF ! of IF(iflag_phys.EQ.2)
1049
1050
1051	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
1052	c$OMP BARRIER
1053	if (pressure_exner) then
1054	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
1055	else
1056	CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
1057	endif
1058	c$OMP BARRIER
1059
1060	cc$OMP END PARALLEL
1061
1062	c-----------------------------------------------------------------------
1063	c dissipation horizontale et verticale des petites echelles:
1064	c ----------------------------------------------------------
1065
1066	IF(apdiss) THEN
1067	cc$OMP PARALLEL DEFAULT(SHARED)
1068	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1069	c$OMP MASTER
1070	call suspend_timer(timer_caldyn)
1071
1072	c print*,'Entree dans la dissipation : Iteration No ',true_itau
1073	c calcul de l'energie cinetique avant dissipation
1074	c print *,'Passage dans la dissipation'
1075
1076	call VTb(VThallo)
1077	c$OMP END MASTER
1078
1079	c$OMP BARRIER
1080
1081	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1082	* jj_Nb_dissip,1,1,Request_dissip)
1083
1084	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1085	* jj_Nb_dissip,1,1,Request_dissip)
1086
1087	call Register_SwapField(teta,teta,ip1jmp1,llm,
1088	* jj_Nb_dissip,Request_dissip)
1089
1090	call Register_SwapField(p,p,ip1jmp1,llmp1,
1091	* jj_Nb_dissip,Request_dissip)
1092
1093	call Register_SwapField(pk,pk,ip1jmp1,llm,
1094	* jj_Nb_dissip,Request_dissip)
1095
1096	call SendRequest(Request_dissip)
1097	c$OMP BARRIER
1098	call WaitRequest(Request_dissip)
1099
1100	c$OMP BARRIER
1101	c$OMP MASTER
1102	call SetDistrib(jj_Nb_dissip)
1103	call VTe(VThallo)
1104	call VTb(VTdissipation)
1105	call start_timer(timer_dissip)
1106	c$OMP END MASTER
1107	c$OMP BARRIER
1108
1109	call covcont_p(llm,ucov,vcov,ucont,vcont)
1110	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1111
1112	c dissipation
1113
1114	! CALL FTRACE_REGION_BEGIN("dissip")
1115	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1116	! CALL FTRACE_REGION_END("dissip")
1117
1118	ijb=ij_begin
1119	ije=ij_end
1120	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1121	DO l=1,llm
1122	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1123	ENDDO
1124	c$OMP END DO NOWAIT
1125	if (pole_sud) ije=ije-iip1
1126	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1127	DO l=1,llm
1128	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1129	ENDDO
1130	c$OMP END DO NOWAIT
1131
1132	c teta=teta+dtetadis
1133
1134
1135	c------------------------------------------------------------------------
1136	if (dissip_conservative) then
1137	C On rajoute la tendance due a la transform. Ec -> E therm. cree
1138	C lors de la dissipation
1139	c$OMP BARRIER
1140	c$OMP MASTER
1141	call suspend_timer(timer_dissip)
1142	call VTb(VThallo)
1143	c$OMP END MASTER
1144	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1145	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1146	call SendRequest(Request_Dissip)
1147	c$OMP BARRIER
1148	call WaitRequest(Request_Dissip)
1149	c$OMP MASTER
1150	call VTe(VThallo)
1151	call resume_timer(timer_dissip)
1152	c$OMP END MASTER
1153	c$OMP BARRIER
1154	call covcont_p(llm,ucov,vcov,ucont,vcont)
1155	call enercin_p(vcov,ucov,vcont,ucont,ecin)
1156
1157	ijb=ij_begin
1158	ije=ij_end
1159	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1160	do l=1,llm
1161	do ij=ijb,ije
1162	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1163	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1164	enddo
1165	enddo
1166	c$OMP END DO NOWAIT
1167	endif ! of if (dissip_conservative)
1168
1169	ijb=ij_begin
1170	ije=ij_end
1171	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1172	do l=1,llm
1173	do ij=ijb,ije
1174	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1175	enddo
1176	enddo
1177	c$OMP END DO NOWAIT
1178	c------------------------------------------------------------------------
1179
1180
1181	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1182	c ... Calcul de la valeur moyenne, unique de h aux poles .....
1183	c
1184
1185	ijb=ij_begin
1186	ije=ij_end
1187
1188	if (pole_nord) then
1189	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1190	DO l = 1, llm
1191	DO ij = 1,iim
1192	tppn(ij) = aire( ij ) * teta( ij ,l)
1193	ENDDO
1194	tpn = SSUM(iim,tppn,1)/apoln
1195
1196	DO ij = 1, iip1
1197	teta( ij ,l) = tpn
1198	ENDDO
1199	ENDDO
1200	c$OMP END DO NOWAIT
1201
1202	if (1 == 0) then
1203	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1204	!!! 2) should probably not be here anyway
1205	!!! but are kept for those who would want to revert to previous behaviour
1206	c$OMP MASTER
1207	DO ij = 1,iim
1208	tppn(ij) = aire( ij ) * ps ( ij )
1209	ENDDO
1210	tpn = SSUM(iim,tppn,1)/apoln
1211
1212	DO ij = 1, iip1
1213	ps( ij ) = tpn
1214	ENDDO
1215	c$OMP END MASTER
1216	endif ! of if (1 == 0)
1217	endif ! of of (pole_nord)
1218
1219	if (pole_sud) then
1220	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1221	DO l = 1, llm
1222	DO ij = 1,iim
1223	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1224	ENDDO
1225	tps = SSUM(iim,tpps,1)/apols
1226
1227	DO ij = 1, iip1
1228	teta(ij+ip1jm,l) = tps
1229	ENDDO
1230	ENDDO
1231	c$OMP END DO NOWAIT
1232
1233	if (1 == 0) then
1234	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1235	!!! 2) should probably not be here anyway
1236	!!! but are kept for those who would want to revert to previous behaviour
1237	c$OMP MASTER
1238	DO ij = 1,iim
1239	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1240	ENDDO
1241	tps = SSUM(iim,tpps,1)/apols
1242
1243	DO ij = 1, iip1
1244	ps(ij+ip1jm) = tps
1245	ENDDO
1246	c$OMP END MASTER
1247	endif ! of if (1 == 0)
1248	endif ! of if (pole_sud)
1249
1250
1251	c$OMP BARRIER
1252	c$OMP MASTER
1253	call VTe(VTdissipation)
1254
1255	call stop_timer(timer_dissip)
1256
1257	call VTb(VThallo)
1258	c$OMP END MASTER
1259	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1260	* jj_Nb_caldyn,Request_dissip)
1261
1262	call Register_SwapField(vcov,vcov,ip1jm,llm,
1263	* jj_Nb_caldyn,Request_dissip)
1264
1265	call Register_SwapField(teta,teta,ip1jmp1,llm,
1266	* jj_Nb_caldyn,Request_dissip)
1267
1268	call Register_SwapField(p,p,ip1jmp1,llmp1,
1269	* jj_Nb_caldyn,Request_dissip)
1270
1271	call Register_SwapField(pk,pk,ip1jmp1,llm,
1272	* jj_Nb_caldyn,Request_dissip)
1273
1274	call SendRequest(Request_dissip)
1275	c$OMP BARRIER
1276	call WaitRequest(Request_dissip)
1277
1278	c$OMP BARRIER
1279	c$OMP MASTER
1280	call SetDistrib(jj_Nb_caldyn)
1281	call VTe(VThallo)
1282	call resume_timer(timer_caldyn)
1283	c print *,'fin dissipation'
1284	c$OMP END MASTER
1285	c$OMP BARRIER
1286	END IF ! of IF(apdiss)
1287
1288	cc$OMP END PARALLEL
1289
1290	c ajout debug
1291	c IF( lafin ) then
1292	c abort_message = 'Simulation finished'
1293	c call abort_gcm(modname,abort_message,0)
1294	c ENDIF
1295
1296	c ********************************************************************
1297	c ********************************************************************
1298	c .... fin de l'integration dynamique et physique pour le pas itau ..
1299	c ********************************************************************
1300	c ********************************************************************
1301
1302	c preparation du pas d'integration suivant ......
1303	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1304	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1305	c$OMP MASTER
1306	call stop_timer(timer_caldyn)
1307	c$OMP END MASTER
1308	IF (itau==itaumax) then
1309	c$OMP MASTER
1310	call allgather_timer_average
1311
1312	if (mpi_rank==0) then
1313
1314	print ,'********************************'
1315	print ,'*** TIMER CALDYN ****'
1316	do i=0,mpi_size-1
1317	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1318	& ' : temps moyen :',
1319	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1320	enddo
1321
1322	print ,'********************************'
1323	print ,'*** TIMER VANLEER ****'
1324	do i=0,mpi_size-1
1325	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1326	& ' : temps moyen :',
1327	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1328	enddo
1329
1330	print ,'********************************'
1331	print ,'*** TIMER DISSIP ****'
1332	do i=0,mpi_size-1
1333	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1334	& ' : temps moyen :',
1335	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1336	enddo
1337
1338	print ,'********************************'
1339	print ,'*** TIMER PHYSIC ****'
1340	do i=0,mpi_size-1
1341	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1342	& ' : temps moyen :',
1343	& timer_average(jj_nb_physic(i),timer_physic,i)
1344	enddo
1345
1346	endif
1347
1348	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1349	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1350	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1351	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1352	CALL print_filtre_timer
1353	call fin_getparam
1354	call finalize_parallel
1355	c$OMP END MASTER
1356	c$OMP BARRIER
1357	RETURN
1358	ENDIF
1359
1360	IF ( .NOT.purmats ) THEN
1361	c ........................................................
1362	c .............. schema matsuno + leapfrog ..............
1363	c ........................................................
1364
1365	IF(forward. OR. leapf) THEN
1366	itau= itau + 1
1367	! iday= day_ini+itau/day_step
1368	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1369	! IF(time.GT.1.) THEN
1370	! time = time-1.
1371	! iday = iday+1
1372	! ENDIF
1373	ENDIF
1374
1375
1376	IF( itau. EQ. itaufinp1 ) then
1377
1378	if (flag_verif) then
1379	write(79,*) 'ucov',ucov
1380	write(80,*) 'vcov',vcov
1381	write(81,*) 'teta',teta
1382	write(82,*) 'ps',ps
1383	write(83,*) 'q',q
1384	WRITE(85,*) 'q1 = ',q(:,:,1)
1385	WRITE(86,*) 'q3 = ',q(:,:,3)
1386	endif
1387
1388
1389	c$OMP MASTER
1390	call fin_getparam
1391	call finalize_parallel
1392	c$OMP END MASTER
1393	abort_message = 'Simulation finished'
1394	call abort_gcm(modname,abort_message,0)
1395	RETURN
1396	ENDIF
1397	c-----------------------------------------------------------------------
1398	c ecriture du fichier histoire moyenne:
1399	c -------------------------------------
1400
1401	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1402	c$OMP BARRIER
1403	IF(itau.EQ.itaufin) THEN
1404	iav=1
1405	ELSE
1406	iav=0
1407	ENDIF
1408	#ifdef CPP_IOIPSL
1409	IF (ok_dynzon) THEN
1410	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1411	call SendRequest(TestRequest)
1412	c$OMP BARRIER
1413	call WaitRequest(TestRequest)
1414	c$OMP BARRIER
1415	c$OMP MASTER
1416	! CALL writedynav_p(histaveid, itau,vcov ,
1417	! , ucov,teta,pk,phi,q,masse,ps,phis)
1418
1419	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
1420	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1421	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1422	c$OMP END MASTER
1423	ENDIF !ok_dynzon
1424	#endif
1425	IF (ok_dyn_ave) THEN
1426	!$OMP MASTER
1427	#ifdef CPP_IOIPSL
1428	! Ehouarn: Gather fields and make master send to output
1429	call Gather_Field(vcov,ip1jm,llm,0)
1430	call Gather_Field(ucov,ip1jmp1,llm,0)
1431	call Gather_Field(teta,ip1jmp1,llm,0)
1432	call Gather_Field(pk,ip1jmp1,llm,0)
1433	call Gather_Field(phi,ip1jmp1,llm,0)
1434	do iq=1,nqtot
1435	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1436	enddo
1437	call Gather_Field(masse,ip1jmp1,llm,0)
1438	call Gather_Field(ps,ip1jmp1,1,0)
1439	call Gather_Field(phis,ip1jmp1,1,0)
1440	if (mpi_rank==0) then
1441	CALL writedynav(itau,vcov,
1442	& ucov,teta,pk,phi,q,masse,ps,phis)
1443	endif
1444	#endif
1445	!$OMP END MASTER
1446	ENDIF ! of IF (ok_dyn_ave)
1447	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1448
1449	c-----------------------------------------------------------------------
1450	c ecriture de la bande histoire:
1451	c ------------------------------
1452
1453	IF( MOD(itau,iecri).EQ.0) THEN
1454	! Ehouarn: output only during LF or Backward Matsuno
1455	if (leapf.or.(.not.leapf.and.(.not.forward))) then
1456	c$OMP BARRIER
1457	c$OMP MASTER
1458	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1459
1460	cym unat=0.
1461
1462	ijb=ij_begin
1463	ije=ij_end
1464
1465	if (pole_nord) then
1466	ijb=ij_begin+iip1
1467	unat(1:iip1,:)=0.
1468	endif
1469
1470	if (pole_sud) then
1471	ije=ij_end-iip1
1472	unat(ij_end-iip1+1:ij_end,:)=0.
1473	endif
1474
1475	do l=1,llm
1476	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1477	enddo
1478
1479	ijb=ij_begin
1480	ije=ij_end
1481	if (pole_sud) ije=ij_end-iip1
1482
1483	do l=1,llm
1484	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1485	enddo
1486
1487	#ifdef CPP_IOIPSL
1488	if (ok_dyn_ins) then
1489	! Ehouarn: Gather fields and make master write to output
1490	call Gather_Field(vcov,ip1jm,llm,0)
1491	call Gather_Field(ucov,ip1jmp1,llm,0)
1492	call Gather_Field(teta,ip1jmp1,llm,0)
1493	call Gather_Field(phi,ip1jmp1,llm,0)
1494	do iq=1,nqtot
1495	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1496	enddo
1497	call Gather_Field(masse,ip1jmp1,llm,0)
1498	call Gather_Field(ps,ip1jmp1,1,0)
1499	call Gather_Field(phis,ip1jmp1,1,0)
1500	if (mpi_rank==0) then
1501	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1502	endif
1503	! CALL writehist_p(histid,histvid, itau,vcov,
1504	! & ucov,teta,phi,q,masse,ps,phis)
1505	! or use writefield_p
1506	! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1507	! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1508	! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1509	! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1510	endif ! of if (ok_dyn_ins)
1511	#endif
1512	! For some Grads outputs of fields
1513	if (output_grads_dyn) then
1514	! Ehouarn: hope this works the way I think it does:
1515	call Gather_Field(unat,ip1jmp1,llm,0)
1516	call Gather_Field(vnat,ip1jm,llm,0)
1517	call Gather_Field(teta,ip1jmp1,llm,0)
1518	call Gather_Field(ps,ip1jmp1,1,0)
1519	do iq=1,nqtot
1520	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1521	enddo
1522	if (mpi_rank==0) then
1523	#include "write_grads_dyn.h"
1524	endif
1525	endif ! of if (output_grads_dyn)
1526	c$OMP END MASTER
1527	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1528	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1529
1530	IF(itau.EQ.itaufin) THEN
1531
1532	c$OMP BARRIER
1533	c$OMP MASTER
1534
1535	! if (planet_type.eq."earth") then
1536	! Write an Earth-format restart file
1537	CALL dynredem1_p("restart.nc",0.0,
1538	& vcov,ucov,teta,q,masse,ps)
1539	! endif ! of if (planet_type.eq."earth")
1540
1541	! CLOSE(99)
1542	c$OMP END MASTER
1543	ENDIF ! of IF (itau.EQ.itaufin)
1544
1545	c-----------------------------------------------------------------------
1546	c gestion de l'integration temporelle:
1547	c ------------------------------------
1548
1549	IF( MOD(itau,iperiod).EQ.0 ) THEN
1550	GO TO 1
1551	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1552
1553	IF( forward ) THEN
1554	c fin du pas forward et debut du pas backward
1555
1556	forward = .FALSE.
1557	leapf = .FALSE.
1558	GO TO 2
1559
1560	ELSE
1561	c fin du pas backward et debut du premier pas leapfrog
1562
1563	leapf = .TRUE.
1564	dt = 2.*dtvr
1565	GO TO 2
1566	END IF
1567	ELSE
1568
1569	c ...... pas leapfrog .....
1570
1571	leapf = .TRUE.
1572	dt = 2.*dtvr
1573	GO TO 2
1574	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1575	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1576
1577
1578	ELSE ! of IF (.not.purmats)
1579
1580	c ........................................................
1581	c .............. schema matsuno ...............
1582	c ........................................................
1583	IF( forward ) THEN
1584
1585	itau = itau + 1
1586	! iday = day_ini+itau/day_step
1587	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1588	!
1589	! IF(time.GT.1.) THEN
1590	! time = time-1.
1591	! iday = iday+1
1592	! ENDIF
1593
1594	forward = .FALSE.
1595	IF( itau. EQ. itaufinp1 ) then
1596	c$OMP MASTER
1597	call fin_getparam
1598	call finalize_parallel
1599	c$OMP END MASTER
1600	abort_message = 'Simulation finished'
1601	call abort_gcm(modname,abort_message,0)
1602	RETURN
1603	ENDIF
1604	GO TO 2
1605
1606	ELSE ! of IF(forward) i.e. backward step
1607
1608	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1609	IF(itau.EQ.itaufin) THEN
1610	iav=1
1611	ELSE
1612	iav=0
1613	ENDIF
1614	#ifdef CPP_IOIPSL
1615	IF (ok_dynzon) THEN
1616	c$OMP BARRIER
1617	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1618	call SendRequest(TestRequest)
1619	c$OMP BARRIER
1620	call WaitRequest(TestRequest)
1621	c$OMP BARRIER
1622	c$OMP MASTER
1623	! CALL writedynav_p(histaveid, itau,vcov ,
1624	! , ucov,teta,pk,phi,q,masse,ps,phis)
1625	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1626	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1627	c$OMP END MASTER
1628	END IF !ok_dynzon
1629	#endif
1630	IF (ok_dyn_ave) THEN
1631	!$OMP MASTER
1632	#ifdef CPP_IOIPSL
1633	! Ehouarn: Gather fields and make master send to output
1634	call Gather_Field(vcov,ip1jm,llm,0)
1635	call Gather_Field(ucov,ip1jmp1,llm,0)
1636	call Gather_Field(teta,ip1jmp1,llm,0)
1637	call Gather_Field(pk,ip1jmp1,llm,0)
1638	call Gather_Field(phi,ip1jmp1,llm,0)
1639	do iq=1,nqtot
1640	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1641	enddo
1642	call Gather_Field(masse,ip1jmp1,llm,0)
1643	call Gather_Field(ps,ip1jmp1,1,0)
1644	call Gather_Field(phis,ip1jmp1,1,0)
1645	if (mpi_rank==0) then
1646	CALL writedynav(itau,vcov,
1647	& ucov,teta,pk,phi,q,masse,ps,phis)
1648	endif
1649	#endif
1650	!$OMP END MASTER
1651	ENDIF ! of IF (ok_dyn_ave)
1652
1653	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1654
1655
1656	IF(MOD(itau,iecri ).EQ.0) THEN
1657	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1658	c$OMP BARRIER
1659	c$OMP MASTER
1660	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1661
1662	cym unat=0.
1663	ijb=ij_begin
1664	ije=ij_end
1665
1666	if (pole_nord) then
1667	ijb=ij_begin+iip1
1668	unat(1:iip1,:)=0.
1669	endif
1670
1671	if (pole_sud) then
1672	ije=ij_end-iip1
1673	unat(ij_end-iip1+1:ij_end,:)=0.
1674	endif
1675
1676	do l=1,llm
1677	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1678	enddo
1679
1680	ijb=ij_begin
1681	ije=ij_end
1682	if (pole_sud) ije=ij_end-iip1
1683
1684	do l=1,llm
1685	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1686	enddo
1687
1688	#ifdef CPP_IOIPSL
1689	if (ok_dyn_ins) then
1690	! Ehouarn: Gather fields and make master send to output
1691	call Gather_Field(vcov,ip1jm,llm,0)
1692	call Gather_Field(ucov,ip1jmp1,llm,0)
1693	call Gather_Field(teta,ip1jmp1,llm,0)
1694	call Gather_Field(phi,ip1jmp1,llm,0)
1695	do iq=1,nqtot
1696	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1697	enddo
1698	call Gather_Field(masse,ip1jmp1,llm,0)
1699	call Gather_Field(ps,ip1jmp1,1,0)
1700	call Gather_Field(phis,ip1jmp1,1,0)
1701	if (mpi_rank==0) then
1702	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1703	endif
1704	! CALL writehist_p(histid, histvid, itau,vcov ,
1705	! & ucov,teta,phi,q,masse,ps,phis)
1706	endif ! of if (ok_dyn_ins)
1707	#endif
1708	! For some Grads output (but does it work?)
1709	if (output_grads_dyn) then
1710	call Gather_Field(unat,ip1jmp1,llm,0)
1711	call Gather_Field(vnat,ip1jm,llm,0)
1712	call Gather_Field(teta,ip1jmp1,llm,0)
1713	call Gather_Field(ps,ip1jmp1,1,0)
1714	do iq=1,nqtot
1715	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1716	enddo
1717	c
1718	if (mpi_rank==0) then
1719	#include "write_grads_dyn.h"
1720	endif
1721	endif ! of if (output_grads_dyn)
1722
1723	c$OMP END MASTER
1724	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1725
1726	IF(itau.EQ.itaufin) THEN
1727	! if (planet_type.eq."earth") then
1728	c$OMP MASTER
1729	CALL dynredem1_p("restart.nc",0.0,
1730	. vcov,ucov,teta,q,masse,ps)
1731	c$OMP END MASTER
1732	! endif ! of if (planet_type.eq."earth")
1733	ENDIF ! of IF(itau.EQ.itaufin)
1734
1735	forward = .TRUE.
1736	GO TO 1
1737
1738	ENDIF ! of IF (forward)
1739
1740	END IF ! of IF(.not.purmats)
1741	c$OMP MASTER
1742	call fin_getparam
1743	call finalize_parallel
1744	c$OMP END MASTER
1745	RETURN
1746	END

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

Download in other formats:

Original Format