Context Navigation

leapfrog_p.F @ 2402

Last change on this file since 2402 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

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

File size: 54.8 KB

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

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

Download in other formats:

Original Format