Context Navigation

leapfrog_p.F @ 5160

Last change on this file since 5160 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

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

File size: 54.0 KB

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

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

Download in other formats:

Original Format