Context Navigation

leapfrog_p.F @ 101

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

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

Download in other formats:

Original Format