Context Navigation

leapfrog_p.F @ 1363

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

Added possibility to run in "Shallow Water" mode. SW mode is automatically set if the number of atmospheric layers is 1.
To use SW mode, the model should be compiled without physics (makelmdz_fcm -p nophys) and/or without calls to the physics (i.e. set iflag_phys=0 in gcm.def).

-Updated some definitions and comments in run.def & gcm.def

-Fixed misplaced #ifdef CPP_EARTH in calfis.F + some write(lunout,*)

-Specific initialisation of fields for SW are done in sw_case_williamson91_6 (called by iniacademic, when read_start=.false.)

Checked (on Vargas & Brodie) that these changes don't alter usual bench GCM outputs.

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

File size: 53.1 KB

Line
1	!
2	! $Id: leapfrog_p.F 1363 2010-04-16 09:50:10Z emillour $
3	!
4	c
5	c
6
7	SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
8	& time_0)
9
10	USE misc_mod
11	USE parallel
12	USE times
13	USE mod_hallo
14	USE Bands
15	USE Write_Field
16	USE Write_Field_p
17	USE vampir
18	USE timer_filtre, ONLY : print_filtre_timer
19	USE infotrac
20	USE guide_p_mod, ONLY : guide_main
21	USE getparam
22	USE control_mod
23
24	IMPLICIT NONE
25
26	c ...... Version du 10/01/98 ..........
27
28	c avec coordonnees verticales hybrides
29	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
30
31	c=======================================================================
32	c
33	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
34	c -------
35	c
36	c Objet:
37	c ------
38	c
39	c GCM LMD nouvelle grille
40	c
41	c=======================================================================
42	c
43	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
44	c et possibilite d'appeler une fonction f(y) a derivee tangente
45	c hyperbolique a la place de la fonction a derivee sinusoidale.
46
47	c ... Possibilite de choisir le shema pour l'advection de
48	c q , en modifiant iadv dans traceur.def (10/02) .
49	c
50	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
51	c Pour Van-Leer iadv=10
52	c
53	c-----------------------------------------------------------------------
54	c Declarations:
55	c -------------
56
57	#include "dimensions.h"
58	#include "paramet.h"
59	#include "comconst.h"
60	#include "comdissnew.h"
61	#include "comvert.h"
62	#include "comgeom.h"
63	#include "logic.h"
64	#include "temps.h"
65	#include "ener.h"
66	#include "description.h"
67	#include "serre.h"
68	!#include "com_io_dyn.h"
69	#include "iniprint.h"
70	#include "academic.h"
71
72	INTEGER longcles
73	PARAMETER ( longcles = 20 )
74	REAL clesphy0( longcles )
75
76	real zqmin,zqmax
77	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	IF (offline) THEN
599	Cmaf stokage du flux de masse pour traceurs OFF-LINE
600
601	#ifdef CPP_IOIPSL
602	CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
603	. dtvr, itau)
604	#endif
605
606
607	ENDIF ! of IF (offline)
608	c
609	ENDIF ! of IF( forward. OR . leapf )
610	cc$OMP END PARALLEL
611
612	c-----------------------------------------------------------------------
613	c integrations dynamique et traceurs:
614	c ----------------------------------
615
616	c$OMP MASTER
617	call VTb(VTintegre)
618	c$OMP END MASTER
619	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
620	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
621	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
622	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
623	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
624	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
625	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
626	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
627	cc$OMP PARALLEL DEFAULT(SHARED)
628	c$OMP BARRIER
629	! CALL FTRACE_REGION_BEGIN("integrd")
630
631	CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
632	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
633	$ finvmaold )
634
635	! CALL FTRACE_REGION_END("integrd")
636	c$OMP BARRIER
637	cc$OMP MASTER
638	c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
639	c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
640	c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
641	c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
642	c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
643	c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
644	c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
645	c call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
646	c
647	c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
648	c do j=1,nqtot
649	c call WriteField_p('q'//trim(int2str(j)),
650	c . reshape(q(:,:,j),(/iip1,jmp1,llm/)))
651	c call WriteField_p('dq'//trim(int2str(j)),
652	c . reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
653	c enddo
654	cc$OMP END MASTER
655
656
657	c$OMP MASTER
658	call VTe(VTintegre)
659	c$OMP END MASTER
660	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
661	c
662	c-----------------------------------------------------------------------
663	c calcul des tendances physiques:
664	c -------------------------------
665	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
666	c
667	IF( purmats ) THEN
668	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
669	ELSE
670	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
671	ENDIF
672
673	cc$OMP END PARALLEL
674
675	c
676	c
677	IF( apphys ) THEN
678	c
679	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
680	c
681	cc$OMP PARALLEL DEFAULT(SHARED)
682	cc$OMP+ PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
683
684	c$OMP MASTER
685	call suspend_timer(timer_caldyn)
686
687	if (prt_level >= 10) then
688	write(lunout,*)
689	& 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
690	endif
691	c$OMP END MASTER
692
693	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
694
695	c$OMP BARRIER
696	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
697	c$OMP BARRIER
698	jD_cur = jD_ref + day_ini - day_ref
699	$ + int (itau * dtvr / daysec)
700	jH_cur = jH_ref + &
701	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
702	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
703
704	c rajout debug
705	c lafin = .true.
706
707
708	c Inbterface avec les routines de phylmd (phymars ... )
709	c -----------------------------------------------------
710
711	c+jld
712
713	c Diagnostique de conservation de l'energie : initialisation
714	IF (ip_ebil_dyn.ge.1 ) THEN
715	ztit='bil dyn'
716	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
717	IF (planet_type.eq."earth") THEN
718	CALL diagedyn(ztit,2,1,1,dtphys
719	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
720	ENDIF
721	ENDIF
722	c-jld
723	c$OMP BARRIER
724	c$OMP MASTER
725	call VTb(VThallo)
726	c$OMP END MASTER
727
728	call SetTag(Request_physic,800)
729
730	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
731	* jj_Nb_physic,2,2,Request_physic)
732
733	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
734	* jj_Nb_physic,2,2,Request_physic)
735
736	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
737	* jj_Nb_physic,2,2,Request_physic)
738
739	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
740	* jj_Nb_physic,1,2,Request_physic)
741
742	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
743	* jj_Nb_physic,2,2,Request_physic)
744
745	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
746	* jj_Nb_physic,2,2,Request_physic)
747
748	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
749	* jj_Nb_physic,2,2,Request_physic)
750
751	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
752	* jj_Nb_physic,2,2,Request_physic)
753
754	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
755	* jj_Nb_physic,2,2,Request_physic)
756
757	c call SetDistrib(jj_nb_vanleer)
758	do j=1,nqtot
759
760	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
761	* jj_Nb_physic,2,2,Request_physic)
762	enddo
763
764	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
765	* jj_Nb_physic,2,2,Request_physic)
766
767	call SendRequest(Request_Physic)
768	c$OMP BARRIER
769	call WaitRequest(Request_Physic)
770
771	c$OMP BARRIER
772	c$OMP MASTER
773	call SetDistrib(jj_nb_Physic)
774	call VTe(VThallo)
775
776	call VTb(VTphysiq)
777	c$OMP END MASTER
778	c$OMP BARRIER
779
780	cc$OMP MASTER
781	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
782	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
783	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
784	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
785	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
786	cc$OMP END MASTER
787	cc$OMP BARRIER
788	! CALL FTRACE_REGION_BEGIN("calfis")
789	CALL calfis_p(lafin ,jD_cur, jH_cur,
790	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
791	$ du,dv,dteta,dq,
792	$ flxw,
793	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
794	! CALL FTRACE_REGION_END("calfis")
795	ijb=ij_begin
796	ije=ij_end
797	if ( .not. pole_nord) then
798	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
799	DO l=1,llm
800	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
801	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
802	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
803	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
804	ENDDO
805	c$OMP END DO NOWAIT
806
807	c$OMP MASTER
808	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
809	c$OMP END MASTER
810	endif ! of if ( .not. pole_nord)
811
812	c$OMP BARRIER
813	c$OMP MASTER
814	call SetDistrib(jj_nb_Physic_bis)
815
816	call VTb(VThallo)
817	c$OMP END MASTER
818	c$OMP BARRIER
819
820	call Register_Hallo(dufi,ip1jmp1,llm,
821	* 1,0,0,1,Request_physic)
822
823	call Register_Hallo(dvfi,ip1jm,llm,
824	* 1,0,0,1,Request_physic)
825
826	call Register_Hallo(dtetafi,ip1jmp1,llm,
827	* 1,0,0,1,Request_physic)
828
829	call Register_Hallo(dpfi,ip1jmp1,1,
830	* 1,0,0,1,Request_physic)
831
832	do j=1,nqtot
833	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
834	* 1,0,0,1,Request_physic)
835	enddo
836
837	call SendRequest(Request_Physic)
838	c$OMP BARRIER
839	call WaitRequest(Request_Physic)
840
841	c$OMP BARRIER
842	c$OMP MASTER
843	call VTe(VThallo)
844
845	call SetDistrib(jj_nb_Physic)
846	c$OMP END MASTER
847	c$OMP BARRIER
848	ijb=ij_begin
849	if (.not. pole_nord) then
850
851	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
852	DO l=1,llm
853	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
854	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
855	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
856	& +dtetafi_tmp(1:iip1,l)
857	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
858	& + dqfi_tmp(1:iip1,l,:)
859	ENDDO
860	c$OMP END DO NOWAIT
861
862	c$OMP MASTER
863	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
864	c$OMP END MASTER
865
866	endif ! of if (.not. pole_nord)
867	c$OMP BARRIER
868	cc$OMP MASTER
869	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
870	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
871	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
872	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
873	cc$OMP END MASTER
874	c
875	c do j=1,nqtot
876	c call WriteField_p('dqfi'//trim(int2str(j)),
877	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
878	c enddo
879
880	c ajout des tendances physiques:
881	c ------------------------------
882	IF (ok_strato) THEN
883	CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
884	ENDIF
885
886	CALL addfi_p( dtphys, leapf, forward ,
887	$ ucov, vcov, teta , q ,ps ,
888	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
889
890	c$OMP BARRIER
891	c$OMP MASTER
892	call VTe(VTphysiq)
893
894	call VTb(VThallo)
895	c$OMP END MASTER
896
897	call SetTag(Request_physic,800)
898	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
899	* jj_Nb_caldyn,Request_physic)
900
901	call Register_SwapField(vcov,vcov,ip1jm,llm,
902	* jj_Nb_caldyn,Request_physic)
903
904	call Register_SwapField(teta,teta,ip1jmp1,llm,
905	* jj_Nb_caldyn,Request_physic)
906
907	call Register_SwapField(masse,masse,ip1jmp1,llm,
908	* jj_Nb_caldyn,Request_physic)
909
910	call Register_SwapField(p,p,ip1jmp1,llmp1,
911	* jj_Nb_caldyn,Request_physic)
912
913	call Register_SwapField(pk,pk,ip1jmp1,llm,
914	* jj_Nb_caldyn,Request_physic)
915
916	call Register_SwapField(phis,phis,ip1jmp1,1,
917	* jj_Nb_caldyn,Request_physic)
918
919	call Register_SwapField(phi,phi,ip1jmp1,llm,
920	* jj_Nb_caldyn,Request_physic)
921
922	call Register_SwapField(w,w,ip1jmp1,llm,
923	* jj_Nb_caldyn,Request_physic)
924
925	do j=1,nqtot
926
927	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
928	* jj_Nb_caldyn,Request_physic)
929
930	enddo
931
932	call SendRequest(Request_Physic)
933	c$OMP BARRIER
934	call WaitRequest(Request_Physic)
935
936	c$OMP BARRIER
937	c$OMP MASTER
938	call VTe(VThallo)
939	call SetDistrib(jj_Nb_caldyn)
940	c$OMP END MASTER
941	c$OMP BARRIER
942	c
943	c Diagnostique de conservation de l'energie : difference
944	IF (ip_ebil_dyn.ge.1 ) THEN
945	ztit='bil phys'
946	CALL diagedyn(ztit,2,1,1,dtphys
947	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
948	ENDIF
949
950	cc$OMP MASTER
951	c if (debug) then
952	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
953	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
954	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
955	c endif
956	cc$OMP END MASTER
957
958
959	c-jld
960	c$OMP MASTER
961	call resume_timer(timer_caldyn)
962	if (FirstPhysic) then
963	ok_start_timer=.TRUE.
964	FirstPhysic=.false.
965	endif
966	c$OMP END MASTER
967	ENDIF ! of IF( apphys )
968
969	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
970	c Calcul academique de la physique = Rappel Newtonien + fritcion
971	c --------------------------------------------------------------
972	cym teta(:,:)=teta(:,:)
973	cym s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
974	ijb=ij_begin
975	ije=ij_end
976	teta(ijb:ije,:)=teta(ijb:ije,:)
977	s -iphysiqdtvr(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
978
979	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
980	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
981	call SendRequest(Request_Physic)
982	c$OMP BARRIER
983	call WaitRequest(Request_Physic)
984	c$OMP BARRIER
985	call friction_p(ucov,vcov,iphysiq*dtvr)
986	ENDIF ! of IF(iflag_phys.EQ.2)
987
988
989	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
990	c$OMP BARRIER
991	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
992	c$OMP BARRIER
993
994	cc$OMP END PARALLEL
995
996	c-----------------------------------------------------------------------
997	c dissipation horizontale et verticale des petites echelles:
998	c ----------------------------------------------------------
999
1000	IF(apdiss) THEN
1001	cc$OMP PARALLEL DEFAULT(SHARED)
1002	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1003	c$OMP MASTER
1004	call suspend_timer(timer_caldyn)
1005
1006	c print*,'Entree dans la dissipation : Iteration No ',true_itau
1007	c calcul de l'energie cinetique avant dissipation
1008	c print *,'Passage dans la dissipation'
1009
1010	call VTb(VThallo)
1011	c$OMP END MASTER
1012
1013	c$OMP BARRIER
1014
1015	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1016	* jj_Nb_dissip,1,1,Request_dissip)
1017
1018	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1019	* jj_Nb_dissip,1,1,Request_dissip)
1020
1021	call Register_SwapField(teta,teta,ip1jmp1,llm,
1022	* jj_Nb_dissip,Request_dissip)
1023
1024	call Register_SwapField(p,p,ip1jmp1,llmp1,
1025	* jj_Nb_dissip,Request_dissip)
1026
1027	call Register_SwapField(pk,pk,ip1jmp1,llm,
1028	* jj_Nb_dissip,Request_dissip)
1029
1030	call SendRequest(Request_dissip)
1031	c$OMP BARRIER
1032	call WaitRequest(Request_dissip)
1033
1034	c$OMP BARRIER
1035	c$OMP MASTER
1036	call SetDistrib(jj_Nb_dissip)
1037	call VTe(VThallo)
1038	call VTb(VTdissipation)
1039	call start_timer(timer_dissip)
1040	c$OMP END MASTER
1041	c$OMP BARRIER
1042
1043	call covcont_p(llm,ucov,vcov,ucont,vcont)
1044	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1045
1046	c dissipation
1047
1048	! CALL FTRACE_REGION_BEGIN("dissip")
1049	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1050	! CALL FTRACE_REGION_END("dissip")
1051
1052	ijb=ij_begin
1053	ije=ij_end
1054	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1055	DO l=1,llm
1056	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1057	ENDDO
1058	c$OMP END DO NOWAIT
1059	if (pole_sud) ije=ije-iip1
1060	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1061	DO l=1,llm
1062	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1063	ENDDO
1064	c$OMP END DO NOWAIT
1065
1066	c teta=teta+dtetadis
1067
1068
1069	c------------------------------------------------------------------------
1070	if (dissip_conservative) then
1071	C On rajoute la tendance due a la transform. Ec -> E therm. cree
1072	C lors de la dissipation
1073	c$OMP BARRIER
1074	c$OMP MASTER
1075	call suspend_timer(timer_dissip)
1076	call VTb(VThallo)
1077	c$OMP END MASTER
1078	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1079	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1080	call SendRequest(Request_Dissip)
1081	c$OMP BARRIER
1082	call WaitRequest(Request_Dissip)
1083	c$OMP MASTER
1084	call VTe(VThallo)
1085	call resume_timer(timer_dissip)
1086	c$OMP END MASTER
1087	c$OMP BARRIER
1088	call covcont_p(llm,ucov,vcov,ucont,vcont)
1089	call enercin_p(vcov,ucov,vcont,ucont,ecin)
1090
1091	ijb=ij_begin
1092	ije=ij_end
1093	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1094	do l=1,llm
1095	do ij=ijb,ije
1096	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1097	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1098	enddo
1099	enddo
1100	c$OMP END DO NOWAIT
1101	endif ! of if (dissip_conservative)
1102
1103	ijb=ij_begin
1104	ije=ij_end
1105	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1106	do l=1,llm
1107	do ij=ijb,ije
1108	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1109	enddo
1110	enddo
1111	c$OMP END DO NOWAIT
1112	c------------------------------------------------------------------------
1113
1114
1115	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1116	c ... Calcul de la valeur moyenne, unique de h aux poles .....
1117	c
1118
1119	ijb=ij_begin
1120	ije=ij_end
1121
1122	if (pole_nord) then
1123	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1124	DO l = 1, llm
1125	DO ij = 1,iim
1126	tppn(ij) = aire( ij ) * teta( ij ,l)
1127	ENDDO
1128	tpn = SSUM(iim,tppn,1)/apoln
1129
1130	DO ij = 1, iip1
1131	teta( ij ,l) = tpn
1132	ENDDO
1133	ENDDO
1134	c$OMP END DO NOWAIT
1135
1136	c$OMP MASTER
1137	DO ij = 1,iim
1138	tppn(ij) = aire( ij ) * ps ( ij )
1139	ENDDO
1140	tpn = SSUM(iim,tppn,1)/apoln
1141
1142	DO ij = 1, iip1
1143	ps( ij ) = tpn
1144	ENDDO
1145	c$OMP END MASTER
1146	endif
1147
1148	if (pole_sud) then
1149	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1150	DO l = 1, llm
1151	DO ij = 1,iim
1152	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1153	ENDDO
1154	tps = SSUM(iim,tpps,1)/apols
1155
1156	DO ij = 1, iip1
1157	teta(ij+ip1jm,l) = tps
1158	ENDDO
1159	ENDDO
1160	c$OMP END DO NOWAIT
1161
1162	c$OMP MASTER
1163	DO ij = 1,iim
1164	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1165	ENDDO
1166	tps = SSUM(iim,tpps,1)/apols
1167
1168	DO ij = 1, iip1
1169	ps(ij+ip1jm) = tps
1170	ENDDO
1171	c$OMP END MASTER
1172	endif
1173
1174
1175	c$OMP BARRIER
1176	c$OMP MASTER
1177	call VTe(VTdissipation)
1178
1179	call stop_timer(timer_dissip)
1180
1181	call VTb(VThallo)
1182	c$OMP END MASTER
1183	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1184	* jj_Nb_caldyn,Request_dissip)
1185
1186	call Register_SwapField(vcov,vcov,ip1jm,llm,
1187	* jj_Nb_caldyn,Request_dissip)
1188
1189	call Register_SwapField(teta,teta,ip1jmp1,llm,
1190	* jj_Nb_caldyn,Request_dissip)
1191
1192	call Register_SwapField(p,p,ip1jmp1,llmp1,
1193	* jj_Nb_caldyn,Request_dissip)
1194
1195	call Register_SwapField(pk,pk,ip1jmp1,llm,
1196	* jj_Nb_caldyn,Request_dissip)
1197
1198	call SendRequest(Request_dissip)
1199	c$OMP BARRIER
1200	call WaitRequest(Request_dissip)
1201
1202	c$OMP BARRIER
1203	c$OMP MASTER
1204	call SetDistrib(jj_Nb_caldyn)
1205	call VTe(VThallo)
1206	call resume_timer(timer_caldyn)
1207	c print *,'fin dissipation'
1208	c$OMP END MASTER
1209	c$OMP BARRIER
1210	END IF ! of IF(apdiss)
1211
1212	cc$OMP END PARALLEL
1213
1214	c ajout debug
1215	c IF( lafin ) then
1216	c abort_message = 'Simulation finished'
1217	c call abort_gcm(modname,abort_message,0)
1218	c ENDIF
1219
1220	c ********************************************************************
1221	c ********************************************************************
1222	c .... fin de l'integration dynamique et physique pour le pas itau ..
1223	c ********************************************************************
1224	c ********************************************************************
1225
1226	c preparation du pas d'integration suivant ......
1227	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1228	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1229	c$OMP MASTER
1230	call stop_timer(timer_caldyn)
1231	c$OMP END MASTER
1232	IF (itau==itaumax) then
1233	c$OMP MASTER
1234	call allgather_timer_average
1235
1236	if (mpi_rank==0) then
1237
1238	print ,'********************************'
1239	print ,'*** TIMER CALDYN ****'
1240	do i=0,mpi_size-1
1241	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1242	& ' : temps moyen :',
1243	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1244	enddo
1245
1246	print ,'********************************'
1247	print ,'*** TIMER VANLEER ****'
1248	do i=0,mpi_size-1
1249	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1250	& ' : temps moyen :',
1251	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1252	enddo
1253
1254	print ,'********************************'
1255	print ,'*** TIMER DISSIP ****'
1256	do i=0,mpi_size-1
1257	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1258	& ' : temps moyen :',
1259	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1260	enddo
1261
1262	print ,'********************************'
1263	print ,'*** TIMER PHYSIC ****'
1264	do i=0,mpi_size-1
1265	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1266	& ' : temps moyen :',
1267	& timer_average(jj_nb_physic(i),timer_physic,i)
1268	enddo
1269
1270	endif
1271
1272	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1273	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1274	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1275	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1276	CALL print_filtre_timer
1277	call fin_getparam
1278	call finalize_parallel
1279	c$OMP END MASTER
1280	c$OMP BARRIER
1281	RETURN
1282	ENDIF
1283
1284	IF ( .NOT.purmats ) THEN
1285	c ........................................................
1286	c .............. schema matsuno + leapfrog ..............
1287	c ........................................................
1288
1289	IF(forward. OR. leapf) THEN
1290	itau= itau + 1
1291	! iday= day_ini+itau/day_step
1292	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1293	! IF(time.GT.1.) THEN
1294	! time = time-1.
1295	! iday = iday+1
1296	! ENDIF
1297	ENDIF
1298
1299
1300	IF( itau. EQ. itaufinp1 ) then
1301
1302	if (flag_verif) then
1303	write(79,*) 'ucov',ucov
1304	write(80,*) 'vcov',vcov
1305	write(81,*) 'teta',teta
1306	write(82,*) 'ps',ps
1307	write(83,*) 'q',q
1308	WRITE(85,*) 'q1 = ',q(:,:,1)
1309	WRITE(86,*) 'q3 = ',q(:,:,3)
1310	endif
1311
1312
1313	c$OMP MASTER
1314	call fin_getparam
1315	call finalize_parallel
1316	c$OMP END MASTER
1317	abort_message = 'Simulation finished'
1318	call abort_gcm(modname,abort_message,0)
1319	RETURN
1320	ENDIF
1321	c-----------------------------------------------------------------------
1322	c ecriture du fichier histoire moyenne:
1323	c -------------------------------------
1324
1325	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1326	c$OMP BARRIER
1327	IF(itau.EQ.itaufin) THEN
1328	iav=1
1329	ELSE
1330	iav=0
1331	ENDIF
1332	#ifdef CPP_IOIPSL
1333	IF (ok_dynzon) THEN
1334	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1335	call SendRequest(TestRequest)
1336	c$OMP BARRIER
1337	call WaitRequest(TestRequest)
1338	c$OMP BARRIER
1339	c$OMP MASTER
1340	! CALL writedynav_p(histaveid, itau,vcov ,
1341	! , ucov,teta,pk,phi,q,masse,ps,phis)
1342
1343	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
1344	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1345	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1346	c$OMP END MASTER
1347	ENDIF !ok_dynzon
1348	#endif
1349	IF (ok_dyn_ave) THEN
1350	!$OMP MASTER
1351	#ifdef CPP_IOIPSL
1352	! Ehouarn: Gather fields and make master send to output
1353	call Gather_Field(vcov,ip1jm,llm,0)
1354	call Gather_Field(ucov,ip1jmp1,llm,0)
1355	call Gather_Field(teta,ip1jmp1,llm,0)
1356	call Gather_Field(pk,ip1jmp1,llm,0)
1357	call Gather_Field(phi,ip1jmp1,llm,0)
1358	do iq=1,nqtot
1359	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1360	enddo
1361	call Gather_Field(masse,ip1jmp1,llm,0)
1362	call Gather_Field(ps,ip1jmp1,1,0)
1363	call Gather_Field(phis,ip1jmp1,1,0)
1364	if (mpi_rank==0) then
1365	CALL writedynav(itau,vcov,
1366	& ucov,teta,pk,phi,q,masse,ps,phis)
1367	endif
1368	#endif
1369	!$OMP END MASTER
1370	ENDIF ! of IF (ok_dyn_ave)
1371	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1372
1373	c-----------------------------------------------------------------------
1374	c ecriture de la bande histoire:
1375	c ------------------------------
1376
1377	IF( MOD(itau,iecri).EQ.0) THEN
1378	! Ehouarn: output only during LF or Backward Matsuno
1379	if (leapf.or.(.not.leapf.and.(.not.forward))) then
1380	c$OMP BARRIER
1381	c$OMP MASTER
1382	nbetat = nbetatdem
1383	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1384
1385	cym unat=0.
1386
1387	ijb=ij_begin
1388	ije=ij_end
1389
1390	if (pole_nord) then
1391	ijb=ij_begin+iip1
1392	unat(1:iip1,:)=0.
1393	endif
1394
1395	if (pole_sud) then
1396	ije=ij_end-iip1
1397	unat(ij_end-iip1+1:ij_end,:)=0.
1398	endif
1399
1400	do l=1,llm
1401	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1402	enddo
1403
1404	ijb=ij_begin
1405	ije=ij_end
1406	if (pole_sud) ije=ij_end-iip1
1407
1408	do l=1,llm
1409	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1410	enddo
1411
1412	#ifdef CPP_IOIPSL
1413	if (ok_dyn_ins) then
1414	! Ehouarn: Gather fields and make master write to output
1415	call Gather_Field(vcov,ip1jm,llm,0)
1416	call Gather_Field(ucov,ip1jmp1,llm,0)
1417	call Gather_Field(teta,ip1jmp1,llm,0)
1418	call Gather_Field(phi,ip1jmp1,llm,0)
1419	do iq=1,nqtot
1420	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1421	enddo
1422	call Gather_Field(masse,ip1jmp1,llm,0)
1423	call Gather_Field(ps,ip1jmp1,1,0)
1424	call Gather_Field(phis,ip1jmp1,1,0)
1425	if (mpi_rank==0) then
1426	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1427	endif
1428	! CALL writehist_p(histid,histvid, itau,vcov,
1429	! & ucov,teta,phi,q,masse,ps,phis)
1430	! or use writefield_p
1431	! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1432	! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1433	! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1434	! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1435	endif ! of if (ok_dyn_ins)
1436	#endif
1437	! For some Grads outputs of fields
1438	if (output_grads_dyn) then
1439	! Ehouarn: hope this works the way I think it does:
1440	call Gather_Field(unat,ip1jmp1,llm,0)
1441	call Gather_Field(vnat,ip1jm,llm,0)
1442	call Gather_Field(teta,ip1jmp1,llm,0)
1443	call Gather_Field(ps,ip1jmp1,1,0)
1444	do iq=1,nqtot
1445	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1446	enddo
1447	if (mpi_rank==0) then
1448	#include "write_grads_dyn.h"
1449	endif
1450	endif ! of if (output_grads_dyn)
1451	c$OMP END MASTER
1452	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1453	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1454
1455	IF(itau.EQ.itaufin) THEN
1456
1457	c$OMP BARRIER
1458	c$OMP MASTER
1459
1460	if (planet_type.eq."earth") then
1461	! Write an Earth-format restart file
1462	CALL dynredem1_p("restart.nc",0.0,
1463	& vcov,ucov,teta,q,masse,ps)
1464	endif ! of if (planet_type.eq."earth")
1465
1466	! CLOSE(99)
1467	c$OMP END MASTER
1468	ENDIF ! of IF (itau.EQ.itaufin)
1469
1470	c-----------------------------------------------------------------------
1471	c gestion de l'integration temporelle:
1472	c ------------------------------------
1473
1474	IF( MOD(itau,iperiod).EQ.0 ) THEN
1475	GO TO 1
1476	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1477
1478	IF( forward ) THEN
1479	c fin du pas forward et debut du pas backward
1480
1481	forward = .FALSE.
1482	leapf = .FALSE.
1483	GO TO 2
1484
1485	ELSE
1486	c fin du pas backward et debut du premier pas leapfrog
1487
1488	leapf = .TRUE.
1489	dt = 2.*dtvr
1490	GO TO 2
1491	END IF
1492	ELSE
1493
1494	c ...... pas leapfrog .....
1495
1496	leapf = .TRUE.
1497	dt = 2.*dtvr
1498	GO TO 2
1499	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1500	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1501
1502
1503	ELSE ! of IF (.not.purmats)
1504
1505	c ........................................................
1506	c .............. schema matsuno ...............
1507	c ........................................................
1508	IF( forward ) THEN
1509
1510	itau = itau + 1
1511	! iday = day_ini+itau/day_step
1512	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1513	!
1514	! IF(time.GT.1.) THEN
1515	! time = time-1.
1516	! iday = iday+1
1517	! ENDIF
1518
1519	forward = .FALSE.
1520	IF( itau. EQ. itaufinp1 ) then
1521	c$OMP MASTER
1522	call fin_getparam
1523	call finalize_parallel
1524	c$OMP END MASTER
1525	abort_message = 'Simulation finished'
1526	call abort_gcm(modname,abort_message,0)
1527	RETURN
1528	ENDIF
1529	GO TO 2
1530
1531	ELSE ! of IF(forward) i.e. backward step
1532
1533	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1534	IF(itau.EQ.itaufin) THEN
1535	iav=1
1536	ELSE
1537	iav=0
1538	ENDIF
1539	#ifdef CPP_IOIPSL
1540	IF (ok_dynzon) THEN
1541	c$OMP BARRIER
1542	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1543	call SendRequest(TestRequest)
1544	c$OMP BARRIER
1545	call WaitRequest(TestRequest)
1546	c$OMP BARRIER
1547	c$OMP MASTER
1548	! CALL writedynav_p(histaveid, itau,vcov ,
1549	! , ucov,teta,pk,phi,q,masse,ps,phis)
1550	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1551	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1552	c$OMP END MASTER
1553	END IF !ok_dynzon
1554	#endif
1555	IF (ok_dyn_ave) THEN
1556	!$OMP MASTER
1557	#ifdef CPP_IOIPSL
1558	! Ehouarn: Gather fields and make master send to output
1559	call Gather_Field(vcov,ip1jm,llm,0)
1560	call Gather_Field(ucov,ip1jmp1,llm,0)
1561	call Gather_Field(teta,ip1jmp1,llm,0)
1562	call Gather_Field(pk,ip1jmp1,llm,0)
1563	call Gather_Field(phi,ip1jmp1,llm,0)
1564	do iq=1,nqtot
1565	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1566	enddo
1567	call Gather_Field(masse,ip1jmp1,llm,0)
1568	call Gather_Field(ps,ip1jmp1,1,0)
1569	call Gather_Field(phis,ip1jmp1,1,0)
1570	if (mpi_rank==0) then
1571	CALL writedynav(itau,vcov,
1572	& ucov,teta,pk,phi,q,masse,ps,phis)
1573	endif
1574	#endif
1575	!$OMP END MASTER
1576	ENDIF ! of IF (ok_dyn_ave)
1577
1578	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1579
1580
1581	IF(MOD(itau,iecri ).EQ.0) THEN
1582	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1583	c$OMP BARRIER
1584	c$OMP MASTER
1585	nbetat = nbetatdem
1586	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1587
1588	cym unat=0.
1589	ijb=ij_begin
1590	ije=ij_end
1591
1592	if (pole_nord) then
1593	ijb=ij_begin+iip1
1594	unat(1:iip1,:)=0.
1595	endif
1596
1597	if (pole_sud) then
1598	ije=ij_end-iip1
1599	unat(ij_end-iip1+1:ij_end,:)=0.
1600	endif
1601
1602	do l=1,llm
1603	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1604	enddo
1605
1606	ijb=ij_begin
1607	ije=ij_end
1608	if (pole_sud) ije=ij_end-iip1
1609
1610	do l=1,llm
1611	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1612	enddo
1613
1614	#ifdef CPP_IOIPSL
1615	if (ok_dyn_ins) then
1616	! Ehouarn: Gather fields and make master send to output
1617	call Gather_Field(vcov,ip1jm,llm,0)
1618	call Gather_Field(ucov,ip1jmp1,llm,0)
1619	call Gather_Field(teta,ip1jmp1,llm,0)
1620	call Gather_Field(phi,ip1jmp1,llm,0)
1621	do iq=1,nqtot
1622	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1623	enddo
1624	call Gather_Field(masse,ip1jmp1,llm,0)
1625	call Gather_Field(ps,ip1jmp1,1,0)
1626	call Gather_Field(phis,ip1jmp1,1,0)
1627	if (mpi_rank==0) then
1628	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1629	endif
1630	! CALL writehist_p(histid, histvid, itau,vcov ,
1631	! & ucov,teta,phi,q,masse,ps,phis)
1632	endif ! of if (ok_dyn_ins)
1633	#endif
1634	! For some Grads output (but does it work?)
1635	if (output_grads_dyn) then
1636	call Gather_Field(unat,ip1jmp1,llm,0)
1637	call Gather_Field(vnat,ip1jm,llm,0)
1638	call Gather_Field(teta,ip1jmp1,llm,0)
1639	call Gather_Field(ps,ip1jmp1,1,0)
1640	do iq=1,nqtot
1641	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1642	enddo
1643	c
1644	if (mpi_rank==0) then
1645	#include "write_grads_dyn.h"
1646	endif
1647	endif ! of if (output_grads_dyn)
1648
1649	c$OMP END MASTER
1650	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1651
1652	IF(itau.EQ.itaufin) THEN
1653	if (planet_type.eq."earth") then
1654	c$OMP MASTER
1655	CALL dynredem1_p("restart.nc",0.0,
1656	. vcov,ucov,teta,q,masse,ps)
1657	c$OMP END MASTER
1658	endif ! of if (planet_type.eq."earth")
1659	ENDIF ! of IF(itau.EQ.itaufin)
1660
1661	forward = .TRUE.
1662	GO TO 1
1663
1664	ENDIF ! of IF (forward)
1665
1666	END IF ! of IF(.not.purmats)
1667	c$OMP MASTER
1668	call fin_getparam
1669	call finalize_parallel
1670	c$OMP END MASTER
1671	RETURN
1672	END

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

Download in other formats:

Original Format