Context Navigation

leapfrog_p.F @ 1825

Last change on this file since 1825 was 1825, checked in by Ehouarn Millour, 11 years ago
Première étape de l'implémentation de XIOS. Modifications isolées dans des flags CPP_XIOS. Sorties opérationnelles (sauf stations et régionalisation) en modes séquentiel et omp, pas mpi. UG ........................................... First step of the XIOS implementation. Modifications are confined into CPP_XIOS flags. Output is operationnal (except for stations and regionalization) in sequential and omp modes (not mpi). UG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.2 KB

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

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

Download in other formats:

Original Format