Context Navigation

leapfrog_p.F @ 1453

Last change on this file since 1453 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 63.6 KB

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

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

Download in other formats:

Original Format