Context Navigation

leapfrog.F @ 5409

Last change on this file since 5409 was 2375, checked in by Ehouarn Millour, 9 years ago
Fix in the computation of the date; the convention is that it corresponds to the time at the end of the current dynamics or physics step (except when in backward Matsuno step where it remains unchanged as it was correctly updated during the forward part of the step). Note that the relationship between itau and date is a bit tricky as itau is incremented between Matsuno forward and backward steps (and from a leapfrog step to the next) but unchanged from Matsuno bacward step to leapfrog step. Because of change in date when calling physics, bench results will change with this revision. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.1 KB

Line
1	!
2	! $Id: leapfrog.F 2375 2015-10-18 06:38:58Z evignon $
3	!
4	c
5	c
6	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
7
8
9	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
10	#ifdef CPP_IOIPSL
11	use IOIPSL
12	#endif
13	USE infotrac, ONLY: nqtot,ok_iso_verif
14	USE guide_mod, ONLY : guide_main
15	USE write_field, ONLY: writefield
16	USE control_mod, ONLY: nday, day_step, planet_type, offline,
17	& iconser, iphysiq, iperiod, dissip_period,
18	& iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
19	& periodav, ok_dyn_ave, output_grads_dyn
20	use exner_hyb_m, only: exner_hyb
21	use exner_milieu_m, only: exner_milieu
22
23	IMPLICIT NONE
24
25	c ...... Version du 10/01/98 ..........
26
27	c avec coordonnees verticales hybrides
28	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
29
30	c=======================================================================
31	c
32	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
33	c -------
34	c
35	c Objet:
36	c ------
37	c
38	c GCM LMD nouvelle grille
39	c
40	c=======================================================================
41	c
42	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
43	c et possibilite d'appeler une fonction f(y) a derivee tangente
44	c hyperbolique a la place de la fonction a derivee sinusoidale.
45
46	c ... Possibilite de choisir le shema pour l'advection de
47	c q , en modifiant iadv dans traceur.def (10/02) .
48	c
49	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
50	c Pour Van-Leer iadv=10
51	c
52	c-----------------------------------------------------------------------
53	c Declarations:
54	c -------------
55
56	#include "dimensions.h"
57	#include "paramet.h"
58	#include "comconst.h"
59	#include "comdissnew.h"
60	#include "comvert.h"
61	#include "comgeom.h"
62	#include "logic.h"
63	#include "temps.h"
64	#include "ener.h"
65	#include "description.h"
66	#include "serre.h"
67	!#include "com_io_dyn.h"
68	#include "iniprint.h"
69	#include "academic.h"
70
71	REAL,INTENT(IN) :: time_0 ! not used
72
73	c dynamical variables:
74	REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
75	REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
76	REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
77	REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
78	REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
79	REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
80	REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
81
82	REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
83	REAL pks(ip1jmp1) ! exner at the surface
84	REAL pk(ip1jmp1,llm) ! exner at mid-layer
85	REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
86	REAL phi(ip1jmp1,llm) ! geopotential
87	REAL w(ip1jmp1,llm) ! vertical velocity
88
89	real zqmin,zqmax
90
91	c variables dynamiques intermediaire pour le transport
92	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
93
94	c variables dynamiques au pas -1
95	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
96	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
97	REAL massem1(ip1jmp1,llm)
98
99	c tendances dynamiques
100	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
101	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
102
103	c tendances de la dissipation
104	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
105	REAL dtetadis(ip1jmp1,llm)
106
107	c tendances physiques
108	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
109	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
110
111	c variables pour le fichier histoire
112	REAL dtav ! intervalle de temps elementaire
113
114	REAL tppn(iim),tpps(iim),tpn,tps
115	c
116	INTEGER itau,itaufinp1,iav
117	! INTEGER iday ! jour julien
118	REAL time
119
120	REAL SSUM
121	! REAL finvmaold(ip1jmp1,llm)
122
123	cym LOGICAL lafin
124	LOGICAL :: lafin=.false.
125	INTEGER ij,iq,l
126	INTEGER ik
127
128	real time_step, t_wrt, t_ops
129
130	! REAL rdayvrai,rdaym_ini
131	! jD_cur: jour julien courant
132	! jH_cur: heure julienne courante
133	REAL :: jD_cur, jH_cur
134	INTEGER :: an, mois, jour
135	REAL :: secondes
136
137	LOGICAL first,callinigrads
138	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
139	save first
140	data first/.true./
141	real dt_cum
142	character*10 infile
143	integer zan, tau0, thoriid
144	integer nid_ctesGCM
145	save nid_ctesGCM
146	real degres
147	real rlong(iip1), rlatg(jjp1)
148	real zx_tmp_2d(iip1,jjp1)
149	integer ndex2d(iip1*jjp1)
150	logical ok_sync
151	parameter (ok_sync = .true.)
152	logical physic
153
154	data callinigrads/.true./
155	character*10 string10
156
157	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
158
159	c+jld variables test conservation energie
160	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
161	C Tendance de la temp. potentiel d (theta)/ d t due a la
162	C tansformation d'energie cinetique en energie thermique
163	C cree par la dissipation
164	REAL dtetaecdt(ip1jmp1,llm)
165	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
166	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
167	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
168	CHARACTER*15 ztit
169	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
170	!IM SAVE ip_ebil_dyn
171	!IM DATA ip_ebil_dyn/0/
172	c-jld
173
174	character*80 dynhist_file, dynhistave_file
175	character(len=*),parameter :: modname="leapfrog"
176	character*80 abort_message
177
178	logical dissip_conservative
179	save dissip_conservative
180	data dissip_conservative/.true./
181
182	LOGICAL prem
183	save prem
184	DATA prem/.true./
185	INTEGER testita
186	PARAMETER (testita = 9)
187
188	logical , parameter :: flag_verif = .false.
189
190
191	integer itau_w ! pas de temps ecriture = itap + itau_phy
192
193
194	if (nday>=0) then
195	itaufin = nday*day_step
196	else
197	itaufin = -nday
198	endif
199	itaufinp1 = itaufin +1
200	itau = 0
201	physic=.true.
202	if (iflag_phys==0.or.iflag_phys==2) physic=.false.
203
204	c iday = day_ini+itau/day_step
205	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
206	c IF(time.GT.1.) THEN
207	c time = time-1.
208	c iday = iday+1
209	c ENDIF
210
211
212	c-----------------------------------------------------------------------
213	c On initialise la pression et la fonction d'Exner :
214	c --------------------------------------------------
215
216	dq(:,:,:)=0.
217	CALL pression ( ip1jmp1, ap, bp, ps, p )
218	if (pressure_exner) then
219	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
220	else
221	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
222	endif
223
224	c-----------------------------------------------------------------------
225	c Debut de l'integration temporelle:
226	c ----------------------------------
227
228	1 CONTINUE ! Matsuno Forward step begins here
229
230	c date: (NB: date remains unchanged for Backward step)
231	c -----
232
233	jD_cur = jD_ref + day_ini - day_ref + &
234	& (itau+1)/day_step
235	jH_cur = jH_ref + start_time + &
236	& mod(itau+1,day_step)/float(day_step)
237	jD_cur = jD_cur + int(jH_cur)
238	jH_cur = jH_cur - int(jH_cur)
239
240	if (ok_iso_verif) then
241	call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
242	endif !if (ok_iso_verif) then
243
244	#ifdef CPP_IOIPSL
245	if (ok_guide) then
246	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
247	endif
248	#endif
249
250
251	c
252	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
253	c CALL test_period ( ucov,vcov,teta,q,p,phis )
254	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
255	c ENDIF
256	c
257
258	! Save fields obtained at previous time step as '...m1'
259	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
260	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
261	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
262	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
263	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
264
265	forward = .TRUE.
266	leapf = .FALSE.
267	dt = dtvr
268
269	c ... P.Le Van .26/04/94 ....
270	! Ehouarn: finvmaold is actually not used
271	! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
272	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
273
274	if (ok_iso_verif) then
275	call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
276	endif !if (ok_iso_verif) then
277
278	2 CONTINUE ! Matsuno backward or leapfrog step begins here
279
280	c-----------------------------------------------------------------------
281
282	c date: (NB: only leapfrog step requires recomputing date)
283	c -----
284
285	IF (leapf) THEN
286	jD_cur = jD_ref + day_ini - day_ref +
287	& (itau+1)/day_step
288	jH_cur = jH_ref + start_time +
289	& mod(itau+1,day_step)/float(day_step)
290	jD_cur = jD_cur + int(jH_cur)
291	jH_cur = jH_cur - int(jH_cur)
292	ENDIF
293
294
295	c gestion des appels de la physique et des dissipations:
296	c ------------------------------------------------------
297	c
298	c ... P.Le Van ( 6/02/95 ) ....
299
300	apphys = .FALSE.
301	statcl = .FALSE.
302	conser = .FALSE.
303	apdiss = .FALSE.
304
305	IF( purmats ) THEN
306	! Purely Matsuno time stepping
307	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
308	IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
309	s apdiss = .TRUE.
310	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
311	s .and. physic ) apphys = .TRUE.
312	ELSE
313	! Leapfrog/Matsuno time stepping
314	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
315	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
316	s apdiss = .TRUE.
317	IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE.
318	END IF
319
320	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
321	! supress dissipation step
322	if (llm.eq.1) then
323	apdiss=.false.
324	endif
325
326
327	if (ok_iso_verif) then
328	call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
329	endif !if (ok_iso_verif) then
330
331	c-----------------------------------------------------------------------
332	c calcul des tendances dynamiques:
333	c --------------------------------
334
335	! compute geopotential phi()
336	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
337
338	time = jD_cur + jH_cur
339	CALL caldyn
340	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
341	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
342
343
344	c-----------------------------------------------------------------------
345	c calcul des tendances advection des traceurs (dont l'humidite)
346	c -------------------------------------------------------------
347
348	if (ok_iso_verif) then
349	call check_isotopes_seq(q,ip1jmp1,
350	& 'leapfrog 686: avant caladvtrac')
351	endif !if (ok_iso_verif) then
352
353	IF( forward. OR . leapf ) THEN
354	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
355	CALL caladvtrac(q,pbaru,pbarv,
356	* p, masse, dq, teta,
357	. flxw, pk)
358	!write(,) 'caladvtrac 346'
359
360
361	IF (offline) THEN
362	Cmaf stokage du flux de masse pour traceurs OFF-LINE
363
364	#ifdef CPP_IOIPSL
365	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
366	. dtvr, itau)
367	#endif
368
369
370	ENDIF ! of IF (offline)
371	c
372	ENDIF ! of IF( forward. OR . leapf )
373
374
375	c-----------------------------------------------------------------------
376	c integrations dynamique et traceurs:
377	c ----------------------------------
378
379	if (ok_iso_verif) then
380	write(,) 'leapfrog 720'
381	call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
382	endif !if (ok_iso_verif) then
383
384	CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
385	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
386	! $ finvmaold )
387
388	if (ok_iso_verif) then
389	write(,) 'leapfrog 724'
390	call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
391	endif !if (ok_iso_verif) then
392
393	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
394	c
395	c-----------------------------------------------------------------------
396	c calcul des tendances physiques:
397	c -------------------------------
398	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
399	c
400	IF( purmats ) THEN
401	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
402	ELSE
403	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
404	ENDIF
405	c
406	c
407	IF( apphys ) THEN
408	c
409	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
410	c
411
412	CALL pression ( ip1jmp1, ap, bp, ps, p )
413	if (pressure_exner) then
414	CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
415	else
416	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
417	endif
418
419	! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
420	! avec dyn3dmem
421	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
422
423	! rdaym_ini = itau * dtvr / daysec
424	! rdayvrai = rdaym_ini + day_ini
425	! jD_cur = jD_ref + day_ini - day_ref
426	! $ + int (itau * dtvr / daysec)
427	! jH_cur = jH_ref + &
428	! & (itau * dtvr / daysec - int(itau * dtvr / daysec))
429	jD_cur = jD_ref + day_ini - day_ref + &
430	& (itau+1)/day_step
431
432	IF (planet_type .eq."generic") THEN
433	! AS: we make jD_cur to be pday
434	jD_cur = int(day_ini + itau/day_step)
435	ENDIF
436
437	jH_cur = jH_ref + start_time + &
438	& mod(itau+1,day_step)/float(day_step)
439	jD_cur = jD_cur + int(jH_cur)
440	jH_cur = jH_cur - int(jH_cur)
441	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
442	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
443	! write(lunout,*)'current date = ',an, mois, jour, secondes
444
445	c rajout debug
446	c lafin = .true.
447
448
449	c Inbterface avec les routines de phylmd (phymars ... )
450	c -----------------------------------------------------
451
452	c+jld
453
454	c Diagnostique de conservation de l'énergie : initialisation
455	IF (ip_ebil_dyn.ge.1 ) THEN
456	ztit='bil dyn'
457	! Ehouarn: be careful, diagedyn is Earth-specific!
458	IF (planet_type.eq."earth") THEN
459	CALL diagedyn(ztit,2,1,1,dtphys
460	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
461	ENDIF
462	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
463	c-jld
464	#ifdef CPP_IOIPSL
465	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
466	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
467	c IF (first) THEN
468	c first=.false.
469	c#include "ini_paramLMDZ_dyn.h"
470	c ENDIF
471	c
472	c#include "write_paramLMDZ_dyn.h"
473	c
474	#endif
475	! #endif of #ifdef CPP_IOIPSL
476	#ifdef CPP_PHYS
477	CALL calfis( lafin , jD_cur, jH_cur,
478	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
479	$ du,dv,dteta,dq,
480	$ flxw,dufi,dvfi,dtetafi,dqfi,dpfi )
481	#endif
482	c ajout des tendances physiques:
483	c ------------------------------
484	CALL addfi( dtphys, leapf, forward ,
485	$ ucov, vcov, teta , q ,ps ,
486	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
487	! since addfi updates ps(), also update p(), masse() and pk()
488	CALL pression (ip1jmp1,ap,bp,ps,p)
489	CALL massdair(p,masse)
490	if (pressure_exner) then
491	CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
492	else
493	CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
494	endif
495
496	IF (ok_strato) THEN
497	CALL top_bound( vcov,ucov,teta,masse,dtphys)
498	ENDIF
499
500	c
501	c Diagnostique de conservation de l'énergie : difference
502	IF (ip_ebil_dyn.ge.1 ) THEN
503	ztit='bil phys'
504	IF (planet_type.eq."earth") THEN
505	CALL diagedyn(ztit,2,1,1,dtphys
506	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
507	ENDIF
508	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
509
510	ENDIF ! of IF( apphys )
511
512	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
513	! Academic case : Simple friction and Newtonan relaxation
514	! -------------------------------------------------------
515	DO l=1,llm
516	DO ij=1,ip1jmp1
517	teta(ij,l)=teta(ij,l)-dtvr*
518	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
519	ENDDO
520	ENDDO ! of DO l=1,llm
521
522	if (planet_type.eq."giant") then
523	! add an intrinsic heat flux at the base of the atmosphere
524	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
525	endif
526
527	call friction(ucov,vcov,dtvr)
528
529	! Sponge layer (if any)
530	IF (ok_strato) THEN
531	! dufi(:,:)=0.
532	! dvfi(:,:)=0.
533	! dtetafi(:,:)=0.
534	! dqfi(:,:,:)=0.
535	! dpfi(:)=0.
536	! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
537	CALL top_bound( vcov,ucov,teta,masse,dtvr)
538	! CALL addfi( dtvr, leapf, forward ,
539	! $ ucov, vcov, teta , q ,ps ,
540	! $ dufi, dvfi, dtetafi , dqfi ,dpfi )
541	ENDIF ! of IF (ok_strato)
542	ENDIF ! of IF (iflag_phys.EQ.2)
543
544
545	c-jld
546
547	CALL pression ( ip1jmp1, ap, bp, ps, p )
548	if (pressure_exner) then
549	CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
550	else
551	CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
552	endif
553	CALL massdair(p,masse)
554
555	if (ok_iso_verif) then
556	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
557	endif !if (ok_iso_verif) then
558
559	c-----------------------------------------------------------------------
560	c dissipation horizontale et verticale des petites echelles:
561	c ----------------------------------------------------------
562
563	IF(apdiss) THEN
564
565
566	c calcul de l'energie cinetique avant dissipation
567	call covcont(llm,ucov,vcov,ucont,vcont)
568	call enercin(vcov,ucov,vcont,ucont,ecin0)
569
570	c dissipation
571	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
572	ucov=ucov+dudis
573	vcov=vcov+dvdis
574	c teta=teta+dtetadis
575
576
577	c------------------------------------------------------------------------
578	if (dissip_conservative) then
579	C On rajoute la tendance due a la transform. Ec -> E therm. cree
580	C lors de la dissipation
581	call covcont(llm,ucov,vcov,ucont,vcont)
582	call enercin(vcov,ucov,vcont,ucont,ecin)
583	dtetaecdt= (ecin0-ecin)/ pk
584	c teta=teta+dtetaecdt
585	dtetadis=dtetadis+dtetaecdt
586	endif
587	teta=teta+dtetadis
588	c------------------------------------------------------------------------
589
590
591	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
592	c ... Calcul de la valeur moyenne, unique de h aux poles .....
593	c
594
595	DO l = 1, llm
596	DO ij = 1,iim
597	tppn(ij) = aire( ij ) * teta( ij ,l)
598	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
599	ENDDO
600	tpn = SSUM(iim,tppn,1)/apoln
601	tps = SSUM(iim,tpps,1)/apols
602
603	DO ij = 1, iip1
604	teta( ij ,l) = tpn
605	teta(ij+ip1jm,l) = tps
606	ENDDO
607	ENDDO
608
609	if (1 == 0) then
610	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
611	!!! 2) should probably not be here anyway
612	!!! but are kept for those who would want to revert to previous behaviour
613	DO ij = 1,iim
614	tppn(ij) = aire( ij ) * ps ( ij )
615	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
616	ENDDO
617	tpn = SSUM(iim,tppn,1)/apoln
618	tps = SSUM(iim,tpps,1)/apols
619
620	DO ij = 1, iip1
621	ps( ij ) = tpn
622	ps(ij+ip1jm) = tps
623	ENDDO
624	endif ! of if (1 == 0)
625
626	END IF ! of IF(apdiss)
627
628	c ajout debug
629	c IF( lafin ) then
630	c abort_message = 'Simulation finished'
631	c call abort_gcm(modname,abort_message,0)
632	c ENDIF
633
634	c ********************************************************************
635	c ********************************************************************
636	c .... fin de l'integration dynamique et physique pour le pas itau ..
637	c ********************************************************************
638	c ********************************************************************
639
640	c preparation du pas d'integration suivant ......
641
642	if (ok_iso_verif) then
643	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
644	endif !if (ok_iso_verif) then
645
646	IF ( .NOT.purmats ) THEN
647	c ........................................................
648	c .............. schema matsuno + leapfrog ..............
649	c ........................................................
650
651	IF(forward. OR. leapf) THEN
652	itau= itau + 1
653	c iday= day_ini+itau/day_step
654	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
655	c IF(time.GT.1.) THEN
656	c time = time-1.
657	c iday = iday+1
658	c ENDIF
659	ENDIF
660
661
662	IF( itau. EQ. itaufinp1 ) then
663	if (flag_verif) then
664	write(79,*) 'ucov',ucov
665	write(80,*) 'vcov',vcov
666	write(81,*) 'teta',teta
667	write(82,*) 'ps',ps
668	write(83,*) 'q',q
669	WRITE(85,*) 'q1 = ',q(:,:,1)
670	WRITE(86,*) 'q3 = ',q(:,:,3)
671	endif
672
673	abort_message = 'Simulation finished'
674
675	call abort_gcm(modname,abort_message,0)
676	ENDIF
677	c-----------------------------------------------------------------------
678	c ecriture du fichier histoire moyenne:
679	c -------------------------------------
680
681	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
682	IF(itau.EQ.itaufin) THEN
683	iav=1
684	ELSE
685	iav=0
686	ENDIF
687
688	IF (ok_dynzon) THEN
689	#ifdef CPP_IOIPSL
690	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
691	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
692	#endif
693	END IF
694	IF (ok_dyn_ave) THEN
695	#ifdef CPP_IOIPSL
696	CALL writedynav(itau,vcov,
697	& ucov,teta,pk,phi,q,masse,ps,phis)
698	#endif
699	ENDIF
700
701	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
702
703	if (ok_iso_verif) then
704	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
705	endif !if (ok_iso_verif) then
706
707	c-----------------------------------------------------------------------
708	c ecriture de la bande histoire:
709	c ------------------------------
710
711	IF( MOD(itau,iecri).EQ.0) THEN
712	! Ehouarn: output only during LF or Backward Matsuno
713	if (leapf.or.(.not.leapf.and.(.not.forward))) then
714	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
715	unat=0.
716	do l=1,llm
717	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
718	vnat(:,l)=vcov(:,l)/cv(:)
719	enddo
720	#ifdef CPP_IOIPSL
721	if (ok_dyn_ins) then
722	! write(lunout,*) "leapfrog: call writehist, itau=",itau
723	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
724	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
725	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
726	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
727	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
728	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
729	endif ! of if (ok_dyn_ins)
730	#endif
731	! For some Grads outputs of fields
732	if (output_grads_dyn) then
733	#include "write_grads_dyn.h"
734	endif
735	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
736	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
737
738	IF(itau.EQ.itaufin) THEN
739
740
741	! if (planet_type.eq."earth") then
742	! Write an Earth-format restart file
743	CALL dynredem1("restart.nc",start_time,
744	& vcov,ucov,teta,q,masse,ps)
745	! endif ! of if (planet_type.eq."earth")
746
747	CLOSE(99)
748	!!! Ehouarn: Why not stop here and now?
749	ENDIF ! of IF (itau.EQ.itaufin)
750
751	c-----------------------------------------------------------------------
752	c gestion de l'integration temporelle:
753	c ------------------------------------
754
755	IF( MOD(itau,iperiod).EQ.0 ) THEN
756	GO TO 1
757	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
758
759	IF( forward ) THEN
760	c fin du pas forward et debut du pas backward
761
762	forward = .FALSE.
763	leapf = .FALSE.
764	GO TO 2
765
766	ELSE
767	c fin du pas backward et debut du premier pas leapfrog
768
769	leapf = .TRUE.
770	dt = 2.*dtvr
771	GO TO 2
772	END IF ! of IF (forward)
773	ELSE
774
775	c ...... pas leapfrog .....
776
777	leapf = .TRUE.
778	dt = 2.*dtvr
779	GO TO 2
780	END IF ! of IF (MOD(itau,iperiod).EQ.0)
781	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
782
783	ELSE ! of IF (.not.purmats)
784
785	if (ok_iso_verif) then
786	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
787	endif !if (ok_iso_verif) then
788
789	c ........................................................
790	c .............. schema matsuno ...............
791	c ........................................................
792	IF( forward ) THEN
793
794	itau = itau + 1
795	c iday = day_ini+itau/day_step
796	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
797	c
798	c IF(time.GT.1.) THEN
799	c time = time-1.
800	c iday = iday+1
801	c ENDIF
802
803	forward = .FALSE.
804	IF( itau. EQ. itaufinp1 ) then
805	abort_message = 'Simulation finished'
806	call abort_gcm(modname,abort_message,0)
807	ENDIF
808	GO TO 2
809
810	ELSE ! of IF(forward) i.e. backward step
811
812	if (ok_iso_verif) then
813	call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
814	endif !if (ok_iso_verif) then
815
816	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
817	IF(itau.EQ.itaufin) THEN
818	iav=1
819	ELSE
820	iav=0
821	ENDIF
822
823	IF (ok_dynzon) THEN
824	#ifdef CPP_IOIPSL
825	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
826	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
827	#endif
828	ENDIF
829	IF (ok_dyn_ave) THEN
830	#ifdef CPP_IOIPSL
831	CALL writedynav(itau,vcov,
832	& ucov,teta,pk,phi,q,masse,ps,phis)
833	#endif
834	ENDIF
835
836	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
837
838	IF(MOD(itau,iecri ).EQ.0) THEN
839	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
840	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
841	unat=0.
842	do l=1,llm
843	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
844	vnat(:,l)=vcov(:,l)/cv(:)
845	enddo
846	#ifdef CPP_IOIPSL
847	if (ok_dyn_ins) then
848	! write(lunout,*) "leapfrog: call writehist (b)",
849	! & itau,iecri
850	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
851	endif ! of if (ok_dyn_ins)
852	#endif
853	! For some Grads outputs
854	if (output_grads_dyn) then
855	#include "write_grads_dyn.h"
856	endif
857
858	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
859
860	IF(itau.EQ.itaufin) THEN
861	! if (planet_type.eq."earth") then
862	CALL dynredem1("restart.nc",start_time,
863	& vcov,ucov,teta,q,masse,ps)
864	! endif ! of if (planet_type.eq."earth")
865	ENDIF ! of IF(itau.EQ.itaufin)
866
867	forward = .TRUE.
868	GO TO 1
869
870	ENDIF ! of IF (forward)
871
872	END IF ! of IF(.not.purmats)
873
874	STOP
875	END

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

Download in other formats:

Original Format