Context Navigation

leapfrog_p.F @ 1229

Last change on this file since 1229 was 1222, checked in by Ehouarn Millour, 15 years ago

Changes and cleanups to enable compiling without physics
and without ioipsl.

IOIPSL related cleanups:

bibio/writehist.F encapsulate the routine (which needs IOIPSL to function)

with #ifdef IOIPSL flag.

dyn3d/abort_gcm.F, dyn3dpar/abort_gcm.F and dyn3dpar/getparam.F90: use ioipsl_getincom module when not compiling with IOIPSL library, in order to always be able to use getin() routine.
removed unused "use IOIPSL" in dyn3dpar/guide_p_mod.F90
calendar related issue: Initialize day_ref and annee_ref in iniacademic.F (i.e. when they are not read from start.nc file)

Earth-specific programs/routines/modules:
create_etat0.F, fluxstokenc.F, limit_netcdf.F, startvar.F
(versions in dyn3d and dyn3dpar)
These routines and modules, which by design and porpose are made to function with
Earth physics are encapsulated with #CPP_EARTH cpp flag.

Earth-specific instructions:

calls to qminimum (specific treatment of first 2 tracers, i.e. water) in dyn3d/caladvtrac.F, dyn3d/integrd.F, dyn3dpar/caladvtrac_p.F, dyn3dpar/integrd_p.F only if (planet_type == 'earth')

Interaction with parallel physics:

routine dyn3dpar/parallel.F90 uses "surface_data" module (which is in the physics ...) to know value of "type_ocean" . Encapsulated that with #ifdef CPP_EARTH and set to a default type_ocean="dummy" otherwise.
So far, only Earth physics are parallelized, so all the interaction between parallel dynamics and parallel physics are encapsulated with #ifdef CCP_EARTH (this way we can run parallel without any physics). The (dyn3dpar) routines which contains such interaction are: bands.F90, gr_dyn_fi_p.F, gr_fi_dyn_p.F, mod_interface_dyn_phys.F90 This should later (when improving dyn/phys interface) be encapsulated with a more general and appropriate #ifdef CPP_PHYS cpp flag.

I checked that these changes do not alter results (on the simple
32x24x11 bench) on Ciclad (seq & mpi), Brodie (seq, mpi & omp) and
Vargas (seq, mpi & omp).

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 1222 2009-08-07 11:48:33Z idelkadi $
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 - day_ref + 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 - day_ref
687	$ + int (itau * dtvr / daysec)
688	jH_cur = jH_ref + &
689	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
690	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
691
692	c rajout debug
693	c lafin = .true.
694
695
696	c Inbterface avec les routines de phylmd (phymars ... )
697	c -----------------------------------------------------
698
699	c+jld
700
701	c Diagnostique de conservation de l'energie : initialisation
702	IF (ip_ebil_dyn.ge.1 ) THEN
703	ztit='bil dyn'
704	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
705	IF (planet_type.eq."earth") THEN
706	CALL diagedyn(ztit,2,1,1,dtphys
707	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
708	ENDIF
709	ENDIF
710	c-jld
711	c$OMP BARRIER
712	c$OMP MASTER
713	call VTb(VThallo)
714	c$OMP END MASTER
715
716	call SetTag(Request_physic,800)
717
718	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
719	* jj_Nb_physic,2,2,Request_physic)
720
721	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
722	* jj_Nb_physic,2,2,Request_physic)
723
724	call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
725	* jj_Nb_physic,2,2,Request_physic)
726
727	call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
728	* jj_Nb_physic,1,2,Request_physic)
729
730	call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
731	* jj_Nb_physic,2,2,Request_physic)
732
733	call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
734	* jj_Nb_physic,2,2,Request_physic)
735
736	call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
737	* jj_Nb_physic,2,2,Request_physic)
738
739	call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
740	* jj_Nb_physic,2,2,Request_physic)
741
742	call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
743	* jj_Nb_physic,2,2,Request_physic)
744
745	c call SetDistrib(jj_nb_vanleer)
746	do j=1,nqtot
747
748	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
749	* jj_Nb_physic,2,2,Request_physic)
750	enddo
751
752	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
753	* jj_Nb_physic,2,2,Request_physic)
754
755	call SendRequest(Request_Physic)
756	c$OMP BARRIER
757	call WaitRequest(Request_Physic)
758
759	c$OMP BARRIER
760	c$OMP MASTER
761	call SetDistrib(jj_nb_Physic)
762	call VTe(VThallo)
763
764	call VTb(VTphysiq)
765	c$OMP END MASTER
766	c$OMP BARRIER
767
768	cc$OMP MASTER
769	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
770	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
771	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
772	c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
773	c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
774	cc$OMP END MASTER
775	cc$OMP BARRIER
776	! CALL FTRACE_REGION_BEGIN("calfis")
777	CALL calfis_p(lafin ,jD_cur, jH_cur,
778	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
779	$ du,dv,dteta,dq,
780	$ flxw,
781	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
782	! CALL FTRACE_REGION_END("calfis")
783	ijb=ij_begin
784	ije=ij_end
785	if ( .not. pole_nord) then
786	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
787	DO l=1,llm
788	dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
789	dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
790	dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
791	dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
792	ENDDO
793	c$OMP END DO NOWAIT
794
795	c$OMP MASTER
796	dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
797	c$OMP END MASTER
798	endif ! of if ( .not. pole_nord)
799
800	c$OMP BARRIER
801	c$OMP MASTER
802	call SetDistrib(jj_nb_Physic_bis)
803
804	call VTb(VThallo)
805	c$OMP END MASTER
806	c$OMP BARRIER
807
808	call Register_Hallo(dufi,ip1jmp1,llm,
809	* 1,0,0,1,Request_physic)
810
811	call Register_Hallo(dvfi,ip1jm,llm,
812	* 1,0,0,1,Request_physic)
813
814	call Register_Hallo(dtetafi,ip1jmp1,llm,
815	* 1,0,0,1,Request_physic)
816
817	call Register_Hallo(dpfi,ip1jmp1,1,
818	* 1,0,0,1,Request_physic)
819
820	do j=1,nqtot
821	call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
822	* 1,0,0,1,Request_physic)
823	enddo
824
825	call SendRequest(Request_Physic)
826	c$OMP BARRIER
827	call WaitRequest(Request_Physic)
828
829	c$OMP BARRIER
830	c$OMP MASTER
831	call VTe(VThallo)
832
833	call SetDistrib(jj_nb_Physic)
834	c$OMP END MASTER
835	c$OMP BARRIER
836	ijb=ij_begin
837	if (.not. pole_nord) then
838
839	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
840	DO l=1,llm
841	dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
842	dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
843	dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
844	& +dtetafi_tmp(1:iip1,l)
845	dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
846	& + dqfi_tmp(1:iip1,l,:)
847	ENDDO
848	c$OMP END DO NOWAIT
849
850	c$OMP MASTER
851	dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
852	c$OMP END MASTER
853
854	endif ! of if (.not. pole_nord)
855	c$OMP BARRIER
856	cc$OMP MASTER
857	c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
858	c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
859	c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
860	c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
861	cc$OMP END MASTER
862	c
863	c do j=1,nqtot
864	c call WriteField_p('dqfi'//trim(int2str(j)),
865	c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
866	c enddo
867
868	c ajout des tendances physiques:
869	c ------------------------------
870	IF (ok_strato) THEN
871	CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
872	ENDIF
873
874	CALL addfi_p( dtphys, leapf, forward ,
875	$ ucov, vcov, teta , q ,ps ,
876	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
877
878	c$OMP BARRIER
879	c$OMP MASTER
880	call VTe(VTphysiq)
881
882	call VTb(VThallo)
883	c$OMP END MASTER
884
885	call SetTag(Request_physic,800)
886	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
887	* jj_Nb_caldyn,Request_physic)
888
889	call Register_SwapField(vcov,vcov,ip1jm,llm,
890	* jj_Nb_caldyn,Request_physic)
891
892	call Register_SwapField(teta,teta,ip1jmp1,llm,
893	* jj_Nb_caldyn,Request_physic)
894
895	call Register_SwapField(masse,masse,ip1jmp1,llm,
896	* jj_Nb_caldyn,Request_physic)
897
898	call Register_SwapField(p,p,ip1jmp1,llmp1,
899	* jj_Nb_caldyn,Request_physic)
900
901	call Register_SwapField(pk,pk,ip1jmp1,llm,
902	* jj_Nb_caldyn,Request_physic)
903
904	call Register_SwapField(phis,phis,ip1jmp1,1,
905	* jj_Nb_caldyn,Request_physic)
906
907	call Register_SwapField(phi,phi,ip1jmp1,llm,
908	* jj_Nb_caldyn,Request_physic)
909
910	call Register_SwapField(w,w,ip1jmp1,llm,
911	* jj_Nb_caldyn,Request_physic)
912
913	do j=1,nqtot
914
915	call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
916	* jj_Nb_caldyn,Request_physic)
917
918	enddo
919
920	call SendRequest(Request_Physic)
921	c$OMP BARRIER
922	call WaitRequest(Request_Physic)
923
924	c$OMP BARRIER
925	c$OMP MASTER
926	call VTe(VThallo)
927	call SetDistrib(jj_Nb_caldyn)
928	c$OMP END MASTER
929	c$OMP BARRIER
930	c
931	c Diagnostique de conservation de l'energie : difference
932	IF (ip_ebil_dyn.ge.1 ) THEN
933	ztit='bil phys'
934	CALL diagedyn(ztit,2,1,1,dtphys
935	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
936	ENDIF
937
938	cc$OMP MASTER
939	c if (debug) then
940	c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
941	c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
942	c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
943	c endif
944	cc$OMP END MASTER
945
946
947	c-jld
948	c$OMP MASTER
949	call resume_timer(timer_caldyn)
950	if (FirstPhysic) then
951	ok_start_timer=.TRUE.
952	FirstPhysic=.false.
953	endif
954	c$OMP END MASTER
955	ENDIF ! of IF( apphys )
956
957	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
958	c Calcul academique de la physique = Rappel Newtonien + fritcion
959	c --------------------------------------------------------------
960	cym teta(:,:)=teta(:,:)
961	cym s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
962	ijb=ij_begin
963	ije=ij_end
964	teta(ijb:ije,:)=teta(ijb:ije,:)
965	s -iphysiqdtvr(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
966
967	call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
968	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
969	call SendRequest(Request_Physic)
970	call WaitRequest(Request_Physic)
971
972	call friction_p(ucov,vcov,iphysiq*dtvr)
973	ENDIF ! of IF(iflag_phys.EQ.2)
974
975
976	CALL pression_p ( ip1jmp1, ap, bp, ps, p )
977	c$OMP BARRIER
978	CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
979	c$OMP BARRIER
980
981	cc$OMP END PARALLEL
982
983	c-----------------------------------------------------------------------
984	c dissipation horizontale et verticale des petites echelles:
985	c ----------------------------------------------------------
986
987	IF(apdiss) THEN
988	cc$OMP PARALLEL DEFAULT(SHARED)
989	cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
990	c$OMP MASTER
991	call suspend_timer(timer_caldyn)
992
993	c print*,'Entree dans la dissipation : Iteration No ',true_itau
994	c calcul de l'energie cinetique avant dissipation
995	c print *,'Passage dans la dissipation'
996
997	call VTb(VThallo)
998	c$OMP END MASTER
999
1000	c$OMP BARRIER
1001
1002	call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
1003	* jj_Nb_dissip,1,1,Request_dissip)
1004
1005	call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
1006	* jj_Nb_dissip,1,1,Request_dissip)
1007
1008	call Register_SwapField(teta,teta,ip1jmp1,llm,
1009	* jj_Nb_dissip,Request_dissip)
1010
1011	call Register_SwapField(p,p,ip1jmp1,llmp1,
1012	* jj_Nb_dissip,Request_dissip)
1013
1014	call Register_SwapField(pk,pk,ip1jmp1,llm,
1015	* jj_Nb_dissip,Request_dissip)
1016
1017	call SendRequest(Request_dissip)
1018	c$OMP BARRIER
1019	call WaitRequest(Request_dissip)
1020
1021	c$OMP BARRIER
1022	c$OMP MASTER
1023	call SetDistrib(jj_Nb_dissip)
1024	call VTe(VThallo)
1025	call VTb(VTdissipation)
1026	call start_timer(timer_dissip)
1027	c$OMP END MASTER
1028	c$OMP BARRIER
1029
1030	call covcont_p(llm,ucov,vcov,ucont,vcont)
1031	call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1032
1033	c dissipation
1034
1035	! CALL FTRACE_REGION_BEGIN("dissip")
1036	CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1037	! CALL FTRACE_REGION_END("dissip")
1038
1039	ijb=ij_begin
1040	ije=ij_end
1041	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1042	DO l=1,llm
1043	ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1044	ENDDO
1045	c$OMP END DO NOWAIT
1046	if (pole_sud) ije=ije-iip1
1047	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1048	DO l=1,llm
1049	vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1050	ENDDO
1051	c$OMP END DO NOWAIT
1052
1053	c teta=teta+dtetadis
1054
1055
1056	c------------------------------------------------------------------------
1057	if (dissip_conservative) then
1058	C On rajoute la tendance due a la transform. Ec -> E therm. cree
1059	C lors de la dissipation
1060	c$OMP BARRIER
1061	c$OMP MASTER
1062	call suspend_timer(timer_dissip)
1063	call VTb(VThallo)
1064	c$OMP END MASTER
1065	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1066	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1067	call SendRequest(Request_Dissip)
1068	c$OMP BARRIER
1069	call WaitRequest(Request_Dissip)
1070	c$OMP MASTER
1071	call VTe(VThallo)
1072	call resume_timer(timer_dissip)
1073	c$OMP END MASTER
1074	c$OMP BARRIER
1075	call covcont_p(llm,ucov,vcov,ucont,vcont)
1076	call enercin_p(vcov,ucov,vcont,ucont,ecin)
1077
1078	ijb=ij_begin
1079	ije=ij_end
1080	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1081	do l=1,llm
1082	do ij=ijb,ije
1083	dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1084	dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1085	enddo
1086	enddo
1087	c$OMP END DO NOWAIT
1088	endif
1089
1090	ijb=ij_begin
1091	ije=ij_end
1092	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1093	do l=1,llm
1094	do ij=ijb,ije
1095	teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1096	enddo
1097	enddo
1098	c$OMP END DO NOWAIT
1099	c------------------------------------------------------------------------
1100
1101
1102	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1103	c ... Calcul de la valeur moyenne, unique de h aux poles .....
1104	c
1105
1106	ijb=ij_begin
1107	ije=ij_end
1108
1109	if (pole_nord) then
1110	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1111	DO l = 1, llm
1112	DO ij = 1,iim
1113	tppn(ij) = aire( ij ) * teta( ij ,l)
1114	ENDDO
1115	tpn = SSUM(iim,tppn,1)/apoln
1116
1117	DO ij = 1, iip1
1118	teta( ij ,l) = tpn
1119	ENDDO
1120	ENDDO
1121	c$OMP END DO NOWAIT
1122
1123	c$OMP MASTER
1124	DO ij = 1,iim
1125	tppn(ij) = aire( ij ) * ps ( ij )
1126	ENDDO
1127	tpn = SSUM(iim,tppn,1)/apoln
1128
1129	DO ij = 1, iip1
1130	ps( ij ) = tpn
1131	ENDDO
1132	c$OMP END MASTER
1133	endif
1134
1135	if (pole_sud) then
1136	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1137	DO l = 1, llm
1138	DO ij = 1,iim
1139	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1140	ENDDO
1141	tps = SSUM(iim,tpps,1)/apols
1142
1143	DO ij = 1, iip1
1144	teta(ij+ip1jm,l) = tps
1145	ENDDO
1146	ENDDO
1147	c$OMP END DO NOWAIT
1148
1149	c$OMP MASTER
1150	DO ij = 1,iim
1151	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
1152	ENDDO
1153	tps = SSUM(iim,tpps,1)/apols
1154
1155	DO ij = 1, iip1
1156	ps(ij+ip1jm) = tps
1157	ENDDO
1158	c$OMP END MASTER
1159	endif
1160
1161
1162	c$OMP BARRIER
1163	c$OMP MASTER
1164	call VTe(VTdissipation)
1165
1166	call stop_timer(timer_dissip)
1167
1168	call VTb(VThallo)
1169	c$OMP END MASTER
1170	call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1171	* jj_Nb_caldyn,Request_dissip)
1172
1173	call Register_SwapField(vcov,vcov,ip1jm,llm,
1174	* jj_Nb_caldyn,Request_dissip)
1175
1176	call Register_SwapField(teta,teta,ip1jmp1,llm,
1177	* jj_Nb_caldyn,Request_dissip)
1178
1179	call Register_SwapField(p,p,ip1jmp1,llmp1,
1180	* jj_Nb_caldyn,Request_dissip)
1181
1182	call Register_SwapField(pk,pk,ip1jmp1,llm,
1183	* jj_Nb_caldyn,Request_dissip)
1184
1185	call SendRequest(Request_dissip)
1186	c$OMP BARRIER
1187	call WaitRequest(Request_dissip)
1188
1189	c$OMP BARRIER
1190	c$OMP MASTER
1191	call SetDistrib(jj_Nb_caldyn)
1192	call VTe(VThallo)
1193	call resume_timer(timer_caldyn)
1194	c print *,'fin dissipation'
1195	c$OMP END MASTER
1196	c$OMP BARRIER
1197	END IF
1198
1199	cc$OMP END PARALLEL
1200
1201	c ajout debug
1202	c IF( lafin ) then
1203	c abort_message = 'Simulation finished'
1204	c call abort_gcm(modname,abort_message,0)
1205	c ENDIF
1206
1207	c ********************************************************************
1208	c ********************************************************************
1209	c .... fin de l'integration dynamique et physique pour le pas itau ..
1210	c ********************************************************************
1211	c ********************************************************************
1212
1213	c preparation du pas d'integration suivant ......
1214	cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1215	cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1216	c$OMP MASTER
1217	call stop_timer(timer_caldyn)
1218	c$OMP END MASTER
1219	IF (itau==itaumax) then
1220	c$OMP MASTER
1221	call allgather_timer_average
1222
1223	if (mpi_rank==0) then
1224
1225	print ,'********************************'
1226	print ,'*** TIMER CALDYN ****'
1227	do i=0,mpi_size-1
1228	print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1229	& ' : temps moyen :',
1230	& timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1231	enddo
1232
1233	print ,'********************************'
1234	print ,'*** TIMER VANLEER ****'
1235	do i=0,mpi_size-1
1236	print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1237	& ' : temps moyen :',
1238	& timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1239	enddo
1240
1241	print ,'********************************'
1242	print ,'*** TIMER DISSIP ****'
1243	do i=0,mpi_size-1
1244	print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1245	& ' : temps moyen :',
1246	& timer_average(jj_nb_dissip(i),timer_dissip,i)
1247	enddo
1248
1249	print ,'********************************'
1250	print ,'*** TIMER PHYSIC ****'
1251	do i=0,mpi_size-1
1252	print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1253	& ' : temps moyen :',
1254	& timer_average(jj_nb_physic(i),timer_physic,i)
1255	enddo
1256
1257	endif
1258
1259	print ,'Taille du Buffer MPI (REAL8)',MaxBufferSize
1260	print ,'Taille du Buffer MPI utilise (REAL8)',MaxBufferSize_Used
1261	print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1262	print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1263	CALL print_filtre_timer
1264	call fin_getparam
1265	call finalize_parallel
1266	c$OMP END MASTER
1267	c$OMP BARRIER
1268	RETURN
1269	ENDIF
1270
1271	IF ( .NOT.purmats ) THEN
1272	c ........................................................
1273	c .............. schema matsuno + leapfrog ..............
1274	c ........................................................
1275
1276	IF(forward. OR. leapf) THEN
1277	itau= itau + 1
1278	! iday= day_ini+itau/day_step
1279	! time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
1280	! IF(time.GT.1.) THEN
1281	! time = time-1.
1282	! iday = iday+1
1283	! ENDIF
1284	ENDIF
1285
1286
1287	IF( itau. EQ. itaufinp1 ) then
1288
1289	c$OMP MASTER
1290	call fin_getparam
1291	call finalize_parallel
1292	c$OMP END MASTER
1293	abort_message = 'Simulation finished'
1294	call abort_gcm(modname,abort_message,0)
1295	RETURN
1296	ENDIF
1297	c-----------------------------------------------------------------------
1298	c ecriture du fichier histoire moyenne:
1299	c -------------------------------------
1300
1301	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1302	c$OMP BARRIER
1303	IF(itau.EQ.itaufin) THEN
1304	iav=1
1305	ELSE
1306	iav=0
1307	ENDIF
1308	#ifdef CPP_IOIPSL
1309	IF (ok_dynzon) THEN
1310	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1311	call SendRequest(TestRequest)
1312	c$OMP BARRIER
1313	call WaitRequest(TestRequest)
1314	c$OMP BARRIER
1315	c$OMP MASTER
1316	CALL writedynav_p(histaveid, itau,vcov ,
1317	, ucov,teta,pk,phi,q,masse,ps,phis)
1318
1319	c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
1320	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1321	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1322	c$OMP END MASTER
1323	ENDIF !ok_dynzon
1324	#endif
1325	ENDIF
1326
1327	c-----------------------------------------------------------------------
1328	c ecriture de la bande histoire:
1329	c ------------------------------
1330
1331	c IF( MOD(itau,iecri ).EQ.0) THEN
1332
1333	IF( MOD(itau,iecri*day_step).EQ.0) THEN
1334	c$OMP BARRIER
1335	c$OMP MASTER
1336	nbetat = nbetatdem
1337	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1338
1339	cym unat=0.
1340
1341	ijb=ij_begin
1342	ije=ij_end
1343
1344	if (pole_nord) then
1345	ijb=ij_begin+iip1
1346	unat(1:iip1,:)=0.
1347	endif
1348
1349	if (pole_sud) then
1350	ije=ij_end-iip1
1351	unat(ij_end-iip1+1:ij_end,:)=0.
1352	endif
1353
1354	do l=1,llm
1355	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1356	enddo
1357
1358	ijb=ij_begin
1359	ije=ij_end
1360	if (pole_sud) ije=ij_end-iip1
1361
1362	do l=1,llm
1363	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1364	enddo
1365
1366	#ifdef CPP_IOIPSL
1367
1368	CALL writehist_p(histid,histvid, itau,vcov,
1369	& ucov,teta,phi,q,masse,ps,phis)
1370
1371	#endif
1372	! For some Grads outputs of fields
1373	if (output_grads_dyn) then
1374	! Ehouarn: hope this works the way I think it does:
1375	call Gather_Field(unat,ip1jmp1,llm,0)
1376	call Gather_Field(vnat,ip1jm,llm,0)
1377	call Gather_Field(teta,ip1jmp1,llm,0)
1378	call Gather_Field(ps,ip1jmp1,1,0)
1379	do iq=1,nqtot
1380	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1381	enddo
1382	if (mpi_rank==0) then
1383	#include "write_grads_dyn.h"
1384	endif
1385	endif ! of if (output_grads_dyn)
1386	c$OMP END MASTER
1387	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1388
1389	IF(itau.EQ.itaufin) THEN
1390
1391	c$OMP BARRIER
1392	c$OMP MASTER
1393
1394	if (planet_type.eq."earth") then
1395	! Write an Earth-format restart file
1396	CALL dynredem1_p("restart.nc",0.0,
1397	& vcov,ucov,teta,q,masse,ps)
1398	endif ! of if (planet_type.eq."earth")
1399
1400	! CLOSE(99)
1401	c$OMP END MASTER
1402	ENDIF ! of IF (itau.EQ.itaufin)
1403
1404	c-----------------------------------------------------------------------
1405	c gestion de l'integration temporelle:
1406	c ------------------------------------
1407
1408	IF( MOD(itau,iperiod).EQ.0 ) THEN
1409	GO TO 1
1410	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1411
1412	IF( forward ) THEN
1413	c fin du pas forward et debut du pas backward
1414
1415	forward = .FALSE.
1416	leapf = .FALSE.
1417	GO TO 2
1418
1419	ELSE
1420	c fin du pas backward et debut du premier pas leapfrog
1421
1422	leapf = .TRUE.
1423	dt = 2.*dtvr
1424	GO TO 2
1425	END IF
1426	ELSE
1427
1428	c ...... pas leapfrog .....
1429
1430	leapf = .TRUE.
1431	dt = 2.*dtvr
1432	GO TO 2
1433	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1434	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1435
1436
1437	ELSE ! of IF (.not.purmats)
1438
1439	c ........................................................
1440	c .............. schema matsuno ...............
1441	c ........................................................
1442	IF( forward ) THEN
1443
1444	itau = itau + 1
1445	! iday = day_ini+itau/day_step
1446	! time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
1447	!
1448	! IF(time.GT.1.) THEN
1449	! time = time-1.
1450	! iday = iday+1
1451	! ENDIF
1452
1453	forward = .FALSE.
1454	IF( itau. EQ. itaufinp1 ) then
1455	c$OMP MASTER
1456	call fin_getparam
1457	call finalize_parallel
1458	c$OMP END MASTER
1459	abort_message = 'Simulation finished'
1460	call abort_gcm(modname,abort_message,0)
1461	RETURN
1462	ENDIF
1463	GO TO 2
1464
1465	ELSE ! of IF(forward)
1466
1467	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1468	IF(itau.EQ.itaufin) THEN
1469	iav=1
1470	ELSE
1471	iav=0
1472	ENDIF
1473	#ifdef CPP_IOIPSL
1474	IF (ok_dynzon) THEN
1475	c$OMP BARRIER
1476
1477	call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1478	call SendRequest(TestRequest)
1479	c$OMP BARRIER
1480	call WaitRequest(TestRequest)
1481
1482	c$OMP BARRIER
1483	c$OMP MASTER
1484	CALL writedynav_p(histaveid, itau,vcov ,
1485	, ucov,teta,pk,phi,q,masse,ps,phis)
1486	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1487	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1488	c$OMP END MASTER
1489	END IF !ok_dynzon
1490	#endif
1491	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1492
1493
1494	c IF(MOD(itau,iecri ).EQ.0) THEN
1495	IF(MOD(itau,iecri*day_step).EQ.0) THEN
1496	c$OMP BARRIER
1497	c$OMP MASTER
1498	nbetat = nbetatdem
1499	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1500
1501	cym unat=0.
1502	ijb=ij_begin
1503	ije=ij_end
1504
1505	if (pole_nord) then
1506	ijb=ij_begin+iip1
1507	unat(1:iip1,:)=0.
1508	endif
1509
1510	if (pole_sud) then
1511	ije=ij_end-iip1
1512	unat(ij_end-iip1+1:ij_end,:)=0.
1513	endif
1514
1515	do l=1,llm
1516	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1517	enddo
1518
1519	ijb=ij_begin
1520	ije=ij_end
1521	if (pole_sud) ije=ij_end-iip1
1522
1523	do l=1,llm
1524	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1525	enddo
1526
1527	#ifdef CPP_IOIPSL
1528
1529	CALL writehist_p(histid, histvid, itau,vcov ,
1530	& ucov,teta,phi,q,masse,ps,phis)
1531	#endif
1532	! For some Grads output (but does it work?)
1533	if (output_grads_dyn) then
1534	call Gather_Field(unat,ip1jmp1,llm,0)
1535	call Gather_Field(vnat,ip1jm,llm,0)
1536	call Gather_Field(teta,ip1jmp1,llm,0)
1537	call Gather_Field(ps,ip1jmp1,1,0)
1538	do iq=1,nqtot
1539	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1540	enddo
1541	c
1542	if (mpi_rank==0) then
1543	#include "write_grads_dyn.h"
1544	endif
1545	endif ! of if (output_grads_dyn)
1546
1547	c$OMP END MASTER
1548	ENDIF ! of IF(MOD(itau,iecri*day_step).EQ.0)
1549
1550	IF(itau.EQ.itaufin) THEN
1551	if (planet_type.eq."earth") then
1552	c$OMP MASTER
1553	CALL dynredem1_p("restart.nc",0.0,
1554	. vcov,ucov,teta,q,masse,ps)
1555	c$OMP END MASTER
1556	endif ! of if (planet_type.eq."earth")
1557	ENDIF ! of IF(itau.EQ.itaufin)
1558
1559	forward = .TRUE.
1560	GO TO 1
1561
1562	ENDIF ! of IF (forward)
1563
1564	END IF ! of IF(.not.purmats)
1565	c$OMP MASTER
1566	call fin_getparam
1567	call finalize_parallel
1568	c$OMP END MASTER
1569	RETURN
1570	END

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

Download in other formats:

Original Format