Context Navigation

leapfrog.F90 @ 5258

Last change on this file since 5258 was 5250, checked in by abarral, 5 weeks ago
Wrap uses of cpp key CPP_PHYS
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: 26.7 KB

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

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

Download in other formats:

Original Format