Context Navigation

leapfrog.F @ 4054

Last change on this file since 4054 was 3947, checked in by Ehouarn Millour, 3 years ago
Nudging: fix the problem of having one extra step in the output guide_ins.nc file. This is done by reseting ok_guide to .false. once the restart file is generated, and thus in the ultimate forward Matsuno step that is done afterwards, nudging is not called. 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: 29.0 KB

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

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

Download in other formats:

Original Format