Context Navigation

leapfrog_p.F @ 2170

Last change on this file since 2170 was 2170, checked in by acozic, 10 years ago
Add modification to manage xios with inca model
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.9 KB

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

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

Download in other formats:

Original Format