Context Navigation

leapfrog_p.F @ 2237

Last change on this file since 2237 was 2221, checked in by Ehouarn Millour, 10 years ago
Some cleanup: remove (unused) clesph0 from dynamics. EM
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.7 KB

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

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

Download in other formats:

Original Format