Context Navigation

leapfrog.F90 @ 5105

Last change on this file since 5105 was 5105, checked in by abarral, 8 weeks ago
Replace 1DUTILS.h by module lmdz_1dutils.f90 Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files) Convert .F to .f90 Fix gradsdef.h formatting Remove unnecessary "RETURN" at the end of functions/subroutines
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: 25.6 KB

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

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

Download in other formats:

Original Format