Context Navigation

leapfrog.F90 @ 5139

Last change on this file since 5139 was 5136, checked in by abarral, 3 months ago
Put comgeom.h, comgeom2.h into modules
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.8 KB

Line
1	! $Id: leapfrog.F90 5136 2024-07-28 14:17:54Z 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 lmdz_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 lmdz_strings, ONLY: msg
26	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
27	USE lmdz_description, ONLY: descript
28	USE lmdz_iniprint, ONLY: lunout, prt_level
29	USE lmdz_ssum_scopy, ONLY: scopy, ssum
30	USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4
31	USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
32	tetagrot, tetatemp, coefdis, vert_prof_dissip
33	USE lmdz_comgeom
34
35	IMPLICIT NONE
36
37	! ...... Version du 10/01/98 ..........
38
39	! avec coordonnees verticales hybrides
40	! avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
41
42	!=======================================================================
43	!
44	! Auteur: P. Le Van /L. Fairhead/F.Hourdin
45	! -------
46	!
47	! Objet:
48	! ------
49	!
50	! GCM LMD nouvelle grille
51	!
52	!=======================================================================
53	!
54	! ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
55	! et possibilite d'appeler une fonction f(y) a derivee tangente
56	! hyperbolique a la place de la fonction a derivee sinusoidale.
57
58	! ... Possibilite de choisir le shema pour l'advection de
59	! q , en modifiant iadv dans traceur.def (10/02) .
60	!
61	! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
62	! Pour Van-Leer iadv=10
63	!
64	!-----------------------------------------------------------------------
65	! Declarations:
66	! -------------
67
68	INCLUDE "dimensions.h"
69	INCLUDE "paramet.h"
70
71	REAL, INTENT(IN) :: time_0 ! not used
72
73	! dynamical variables:
74	REAL, INTENT(INOUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind
75	REAL, INTENT(INOUT) :: vcov(ip1jm, llm) ! meridional covariant wind
76	REAL, INTENT(INOUT) :: teta(ip1jmp1, llm) ! potential temperature
77	REAL, INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
78	REAL, INTENT(INOUT) :: masse(ip1jmp1, llm) ! air mass
79	REAL, INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
80	REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers
81
82	REAL :: p (ip1jmp1, llmp1) ! interlayer pressure
83	REAL :: pks(ip1jmp1) ! exner at the surface
84	REAL :: pk(ip1jmp1, llm) ! exner at mid-layer
85	REAL :: pkf(ip1jmp1, llm) ! filtered exner at mid-layer
86	REAL :: phi(ip1jmp1, llm) ! geopotential
87	REAL :: w(ip1jmp1, llm) ! vertical velocity
88
89	REAL :: zqmin, zqmax
90
91	! variables dynamiques intermediaire pour le transport
92	REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
93
94	! variables dynamiques au pas -1
95	REAL :: vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
96	REAL :: tetam1(ip1jmp1, llm), psm1(ip1jmp1)
97	REAL :: massem1(ip1jmp1, llm)
98
99	! tendances dynamiques
100	REAL :: dv(ip1jm, llm), du(ip1jmp1, llm)
101	REAL :: dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqtot), dp(ip1jmp1)
102
103	! tendances de la dissipation
104	REAL :: dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
105	REAL :: dtetadis(ip1jmp1, llm)
106
107	! tendances physiques
108	REAL :: dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
109	REAL :: dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqtot), dpfi(ip1jmp1)
110
111	! variables pour le fichier histoire
112	REAL :: dtav ! intervalle de temps elementaire
113
114	REAL :: tppn(iim), tpps(iim), tpn, tps
115	!
116	INTEGER :: itau, itaufinp1, iav
117	! INTEGER iday ! jour julien
118	REAL :: time
119
120	! REAL finvmaold(ip1jmp1,llm)
121
122	!ym 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	!IM : pour sortir les param. du modele dans un fis. netcdf 110106
138	save first
139	data first/.TRUE./
140	REAL :: dt_cum
141	CHARACTER(LEN = 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(LEN = 10) :: string10
155
156	REAL :: flxw(ip1jmp1, llm) ! flux de masse verticale
157
158	!+jld variables test conservation energie
159	REAL :: ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
160	! Tendance de la temp. potentiel d (theta)/ d t due a la
161	! tansformation d'energie cinetique en energie thermique
162	! 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(len = 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	!-jld
172
173	CHARACTER(LEN = 80) :: dynhist_file, dynhistave_file
174	CHARACTER(LEN = *), parameter :: modname = "leapfrog"
175	CHARACTER(LEN = 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	INTEGER :: itau_w ! pas de temps ecriture = itap + itau_phy
190
191	IF (nday>=0) THEN
192	itaufin = nday * day_step
193	else
194	itaufin = -nday
195	ENDIF
196	itaufinp1 = itaufin + 1
197	itau = 0
198	physic = .TRUE.
199	IF (iflag_phys==0.OR.iflag_phys==2) physic = .FALSE.
200
201	! iday = day_ini+itau/day_step
202	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
203	! IF(time.GT.1.) THEN
204	! time = time-1.
205	! iday = iday+1
206	! ENDIF
207
208
209	!-----------------------------------------------------------------------
210	! On initialise la pression et la fonction d'Exner :
211	! --------------------------------------------------
212
213	dq(:, :, :) = 0.
214	CALL pression (ip1jmp1, ap, bp, ps, p)
215	IF (pressure_exner) THEN
216	CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
217	else
218	CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
219	ENDIF
220
221	!-----------------------------------------------------------------------
222	! Debut de l'integration temporelle:
223	! ----------------------------------
224
225	1 CONTINUE ! Matsuno Forward step begins here
226
227	! date: (NB: date remains unchanged for Backward step)
228	! -----
229
230	jD_cur = jD_ref + day_ini - day_ref + &
231	(itau + 1) / day_step
232	jH_cur = jH_ref + start_time + &
233	mod(itau + 1, day_step) / float(day_step)
234	jD_cur = jD_cur + int(jH_cur)
235	jH_cur = jH_cur - int(jH_cur)
236
237	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 321')
238
239	IF (ok_guide) THEN
240	CALL guide_main(itau, ucov, vcov, teta, q, masse, ps)
241	ENDIF
242
243
244	!
245	! IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
246	! CALL test_period ( ucov,vcov,teta,q,p,phis )
247	! PRINT *,' ---- Test_period apres continue OK ! -----', itau
248	! ENDIF
249	!
250
251	! Save fields obtained at previous time step as '...m1'
252	CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)
253	CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)
254	CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)
255	CALL SCOPY(ijp1llm, masse, 1, massem1, 1)
256	CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)
257
258	forward = .TRUE.
259	leapf = .FALSE.
260	dt = dtvr
261
262	! ... P.Le Van .26/04/94 ....
263	! Ehouarn: finvmaold is actually not used
264	! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
265	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
266
267	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 400')
268
269	2 CONTINUE ! Matsuno backward or leapfrog step begins here
270
271	!-----------------------------------------------------------------------
272
273	! date: (NB: only leapfrog step requires recomputing date)
274	! -----
275
276	IF (leapf) THEN
277	jD_cur = jD_ref + day_ini - day_ref + &
278	(itau + 1) / day_step
279	jH_cur = jH_ref + start_time + &
280	mod(itau + 1, day_step) / float(day_step)
281	jD_cur = jD_cur + int(jH_cur)
282	jH_cur = jH_cur - int(jH_cur)
283	ENDIF
284
285
286	! gestion des appels de la physique et des dissipations:
287	! ------------------------------------------------------
288	!
289	! ... P.Le Van ( 6/02/95 ) ....
290
291	apphys = .FALSE.
292	statcl = .FALSE.
293	conser = .FALSE.
294	apdiss = .FALSE.
295
296	IF(purmats) THEN
297	! Purely Matsuno time stepping
298	IF(MOD(itau, iconser) ==0.AND. forward) conser = .TRUE.
299	IF(MOD(itau, dissip_period)==0.AND..NOT.forward) &
300	apdiss = .TRUE.
301	IF(MOD(itau, iphysiq)==0.AND..NOT.forward &
302	.AND. physic) apphys = .TRUE.
303	ELSE
304	! Leapfrog/Matsuno time stepping
305	IF(MOD(itau, iconser) == 0) conser = .TRUE.
306	IF(MOD(itau + 1, dissip_period)==0 .AND. .NOT. forward) &
307	apdiss = .TRUE.
308	IF(MOD(itau + 1, iphysiq)==0.AND.physic) apphys = .TRUE.
309	END IF
310
311	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
312	! supress dissipation step
313	IF (llm==1) THEN
314	apdiss = .FALSE.
315	ENDIF
316
317	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 589')
318
319	!-----------------------------------------------------------------------
320	! calcul des tendances dynamiques:
321	! --------------------------------
322
323	! compute geopotential phi()
324	CALL geopot (ip1jmp1, teta, pk, pks, phis, phi)
325
326	time = jD_cur + jH_cur
327	CALL caldyn &
328	(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
329	phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)
330
331
332	!-----------------------------------------------------------------------
333	! calcul des tendances advection des traceurs (dont l'humidite)
334	! -------------------------------------------------------------
335
336	CALL check_isotopes_seq(q, ip1jmp1, &
337	'leapfrog 686: avant caladvtrac')
338
339	IF(forward .OR. leapf) THEN
340	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
341	CALL caladvtrac(q, pbaru, pbarv, &
342	p, masse, dq, teta, &
343	flxw, pk)
344	!WRITE(,) 'caladvtrac 346'
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	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, dtvr * iperiod, dtvr * day_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	! END IF ! 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, dtvr * iperiod, dtvr * day_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	! END IF ! 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