Context Navigation

leapfrog_p.F @ 2176

Last change on this file since 2176 was 2171, checked in by acozic, 10 years ago
There are some commits that we must not do just before holiday .... so be back to rev 2168
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.8 KB

Line
1	!
2	! $Id: leapfrog_p.F 2171 2014-12-19 15:21:08Z jescribano $
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	call finalize_parallel
1413	c$OMP END MASTER
1414	abort_message = 'Simulation finished'
1415	call abort_gcm(modname,abort_message,0)
1416	RETURN
1417	ENDIF
1418	c-----------------------------------------------------------------------
1419	c ecriture du fichier histoire moyenne:
1420	c -------------------------------------
1421
1422	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1423	c$OMP BARRIER
1424	IF(itau.EQ.itaufin) THEN
1425	iav=1
1426	ELSE
1427	iav=0
1428	ENDIF
1429	#ifdef CPP_IOIPSL
1430	IF (ok_dynzon) THEN
1431	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1432	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1433	ENDIF !ok_dynzon
1434	#endif
1435	IF (ok_dyn_ave) THEN
1436	!$OMP MASTER
1437	#ifdef CPP_IOIPSL
1438	! Ehouarn: Gather fields and make master send to output
1439	call Gather_Field(vcov,ip1jm,llm,0)
1440	call Gather_Field(ucov,ip1jmp1,llm,0)
1441	call Gather_Field(teta,ip1jmp1,llm,0)
1442	call Gather_Field(pk,ip1jmp1,llm,0)
1443	call Gather_Field(phi,ip1jmp1,llm,0)
1444	do iq=1,nqtot
1445	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1446	enddo
1447	call Gather_Field(masse,ip1jmp1,llm,0)
1448	call Gather_Field(ps,ip1jmp1,1,0)
1449	call Gather_Field(phis,ip1jmp1,1,0)
1450	if (mpi_rank==0) then
1451	CALL writedynav(itau,vcov,
1452	& ucov,teta,pk,phi,q,masse,ps,phis)
1453	endif
1454	#endif
1455	!$OMP END MASTER
1456	ENDIF ! of IF (ok_dyn_ave)
1457	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1458
1459	c-----------------------------------------------------------------------
1460	c ecriture de la bande histoire:
1461	c ------------------------------
1462
1463	IF( MOD(itau,iecri).EQ.0) THEN
1464	! Ehouarn: output only during LF or Backward Matsuno
1465	if (leapf.or.(.not.leapf.and.(.not.forward))) then
1466	c$OMP BARRIER
1467	c$OMP MASTER
1468	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1469
1470	cym unat=0.
1471
1472	ijb=ij_begin
1473	ije=ij_end
1474
1475	if (pole_nord) then
1476	ijb=ij_begin+iip1
1477	unat(1:iip1,:)=0.
1478	endif
1479
1480	if (pole_sud) then
1481	ije=ij_end-iip1
1482	unat(ij_end-iip1+1:ij_end,:)=0.
1483	endif
1484
1485	do l=1,llm
1486	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1487	enddo
1488
1489	ijb=ij_begin
1490	ije=ij_end
1491	if (pole_sud) ije=ij_end-iip1
1492
1493	do l=1,llm
1494	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1495	enddo
1496
1497	#ifdef CPP_IOIPSL
1498	if (ok_dyn_ins) then
1499	! Ehouarn: Gather fields and make master write to output
1500	call Gather_Field(vcov,ip1jm,llm,0)
1501	call Gather_Field(ucov,ip1jmp1,llm,0)
1502	call Gather_Field(teta,ip1jmp1,llm,0)
1503	call Gather_Field(phi,ip1jmp1,llm,0)
1504	do iq=1,nqtot
1505	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1506	enddo
1507	call Gather_Field(masse,ip1jmp1,llm,0)
1508	call Gather_Field(ps,ip1jmp1,1,0)
1509	call Gather_Field(phis,ip1jmp1,1,0)
1510	if (mpi_rank==0) then
1511	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1512	endif
1513	! CALL writehist_p(histid,histvid, itau,vcov,
1514	! & ucov,teta,phi,q,masse,ps,phis)
1515	! or use writefield_p
1516	! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1517	! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1518	! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1519	! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1520	endif ! of if (ok_dyn_ins)
1521	#endif
1522	! For some Grads outputs of fields
1523	if (output_grads_dyn) then
1524	! Ehouarn: hope this works the way I think it does:
1525	call Gather_Field(unat,ip1jmp1,llm,0)
1526	call Gather_Field(vnat,ip1jm,llm,0)
1527	call Gather_Field(teta,ip1jmp1,llm,0)
1528	call Gather_Field(ps,ip1jmp1,1,0)
1529	do iq=1,nqtot
1530	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1531	enddo
1532	if (mpi_rank==0) then
1533	#include "write_grads_dyn.h"
1534	endif
1535	endif ! of if (output_grads_dyn)
1536	c$OMP END MASTER
1537	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1538	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1539
1540	IF(itau.EQ.itaufin) THEN
1541
1542	c$OMP BARRIER
1543	c$OMP MASTER
1544
1545	! if (planet_type.eq."earth") then
1546	! Write an Earth-format restart file
1547	CALL dynredem1_p("restart.nc",0.0,
1548	& vcov,ucov,teta,q,masse,ps)
1549	! endif ! of if (planet_type.eq."earth")
1550
1551	! CLOSE(99)
1552	c$OMP END MASTER
1553	ENDIF ! of IF (itau.EQ.itaufin)
1554
1555	c-----------------------------------------------------------------------
1556	c gestion de l'integration temporelle:
1557	c ------------------------------------
1558
1559	IF( MOD(itau,iperiod).EQ.0 ) THEN
1560	GO TO 1
1561	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1562
1563	IF( forward ) THEN
1564	c fin du pas forward et debut du pas backward
1565
1566	forward = .FALSE.
1567	leapf = .FALSE.
1568	GO TO 2
1569
1570	ELSE
1571	c fin du pas backward et debut du premier pas leapfrog
1572
1573	leapf = .TRUE.
1574	dt = 2.*dtvr
1575	GO TO 2
1576	END IF
1577	ELSE
1578
1579	c ...... pas leapfrog .....
1580
1581	leapf = .TRUE.
1582	dt = 2.*dtvr
1583	GO TO 2
1584	END IF ! of IF (MOD(itau,iperiod).EQ.0)
1585	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1586
1587
1588	ELSE ! of IF (.not.purmats)
1589
1590	c ........................................................
1591	c .............. schema matsuno ...............
1592	c ........................................................
1593	IF( forward ) THEN
1594
1595	itau = itau + 1
1596	! iday = day_ini+itau/day_step
1597	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1598	!
1599	! IF(time.GT.1.) THEN
1600	! time = time-1.
1601	! iday = iday+1
1602	! ENDIF
1603
1604	forward = .FALSE.
1605	IF( itau. EQ. itaufinp1 ) then
1606	c$OMP MASTER
1607	call fin_getparam
1608	call finalize_parallel
1609	c$OMP END MASTER
1610	abort_message = 'Simulation finished'
1611	call abort_gcm(modname,abort_message,0)
1612	RETURN
1613	ENDIF
1614	GO TO 2
1615
1616	ELSE ! of IF(forward) i.e. backward step
1617
1618	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1619	IF(itau.EQ.itaufin) THEN
1620	iav=1
1621	ELSE
1622	iav=0
1623	ENDIF
1624	#ifdef CPP_IOIPSL
1625	IF (ok_dynzon) THEN
1626	CALL bilan_dyn_p(2,dtvriperiod,dtvrday_step*periodav,
1627	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1628	END IF !ok_dynzon
1629	#endif
1630	IF (ok_dyn_ave) THEN
1631	!$OMP MASTER
1632	#ifdef CPP_IOIPSL
1633	! Ehouarn: Gather fields and make master send to output
1634	call Gather_Field(vcov,ip1jm,llm,0)
1635	call Gather_Field(ucov,ip1jmp1,llm,0)
1636	call Gather_Field(teta,ip1jmp1,llm,0)
1637	call Gather_Field(pk,ip1jmp1,llm,0)
1638	call Gather_Field(phi,ip1jmp1,llm,0)
1639	do iq=1,nqtot
1640	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1641	enddo
1642	call Gather_Field(masse,ip1jmp1,llm,0)
1643	call Gather_Field(ps,ip1jmp1,1,0)
1644	call Gather_Field(phis,ip1jmp1,1,0)
1645	if (mpi_rank==0) then
1646	CALL writedynav(itau,vcov,
1647	& ucov,teta,pk,phi,q,masse,ps,phis)
1648	endif
1649	#endif
1650	!$OMP END MASTER
1651	ENDIF ! of IF (ok_dyn_ave)
1652
1653	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1654
1655
1656	IF(MOD(itau,iecri ).EQ.0) THEN
1657	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1658	c$OMP BARRIER
1659	c$OMP MASTER
1660	CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1661
1662	cym unat=0.
1663	ijb=ij_begin
1664	ije=ij_end
1665
1666	if (pole_nord) then
1667	ijb=ij_begin+iip1
1668	unat(1:iip1,:)=0.
1669	endif
1670
1671	if (pole_sud) then
1672	ije=ij_end-iip1
1673	unat(ij_end-iip1+1:ij_end,:)=0.
1674	endif
1675
1676	do l=1,llm
1677	unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1678	enddo
1679
1680	ijb=ij_begin
1681	ije=ij_end
1682	if (pole_sud) ije=ij_end-iip1
1683
1684	do l=1,llm
1685	vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1686	enddo
1687
1688	#ifdef CPP_IOIPSL
1689	if (ok_dyn_ins) then
1690	! Ehouarn: Gather fields and make master send to output
1691	call Gather_Field(vcov,ip1jm,llm,0)
1692	call Gather_Field(ucov,ip1jmp1,llm,0)
1693	call Gather_Field(teta,ip1jmp1,llm,0)
1694	call Gather_Field(phi,ip1jmp1,llm,0)
1695	do iq=1,nqtot
1696	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1697	enddo
1698	call Gather_Field(masse,ip1jmp1,llm,0)
1699	call Gather_Field(ps,ip1jmp1,1,0)
1700	call Gather_Field(phis,ip1jmp1,1,0)
1701	if (mpi_rank==0) then
1702	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1703	endif
1704	! CALL writehist_p(histid, histvid, itau,vcov ,
1705	! & ucov,teta,phi,q,masse,ps,phis)
1706	endif ! of if (ok_dyn_ins)
1707	#endif
1708	! For some Grads output (but does it work?)
1709	if (output_grads_dyn) then
1710	call Gather_Field(unat,ip1jmp1,llm,0)
1711	call Gather_Field(vnat,ip1jm,llm,0)
1712	call Gather_Field(teta,ip1jmp1,llm,0)
1713	call Gather_Field(ps,ip1jmp1,1,0)
1714	do iq=1,nqtot
1715	call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1716	enddo
1717	c
1718	if (mpi_rank==0) then
1719	#include "write_grads_dyn.h"
1720	endif
1721	endif ! of if (output_grads_dyn)
1722
1723	c$OMP END MASTER
1724	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1725
1726	IF(itau.EQ.itaufin) THEN
1727	! if (planet_type.eq."earth") then
1728	c$OMP MASTER
1729	CALL dynredem1_p("restart.nc",0.0,
1730	. vcov,ucov,teta,q,masse,ps)
1731	c$OMP END MASTER
1732	! endif ! of if (planet_type.eq."earth")
1733	ENDIF ! of IF(itau.EQ.itaufin)
1734
1735	forward = .TRUE.
1736	GO TO 1
1737
1738	ENDIF ! of IF (forward)
1739
1740	END IF ! of IF(.not.purmats)
1741	c$OMP MASTER
1742	call fin_getparam
1743	call finalize_parallel
1744	c$OMP END MASTER
1745	RETURN
1746	END

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

Download in other formats:

Original Format