Context Navigation

leapfrog_p.F @ 1207

Last change on this file since 1207 was 1207, checked in by Laurent Fairhead, 15 years ago
Probleme sur la forme des commentaires laisses par emacs, openmp pourrait mal les comprendre L'argument heure disparait de la liste des arguments d'appel a calfis
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 48.9 KB

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

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

Download in other formats:

Original Format