Context Navigation

leapfrog.F90 @ 5300

Last change on this file since 5300 was 5292, checked in by abarral, 4 days ago
Move academic.h chem.h chem_spla.h to module
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.4 KB

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

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

Download in other formats:

Original Format