Context Navigation

leapfrog.F90 @ 5182

Last change on this file since 5182 was 5182, checked in by abarral, 13 days ago
(WIP) Replace REPROBUS CPP KEY by logical properly name 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 5182 2024-09-10 14:25:29Z 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 lmdz_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	USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
36	USE lmdz_paramet
37	IMPLICIT NONE
38
39	! ...... Version du 10/01/98 ..........
40
41	! avec coordonnees verticales hybrides
42	! avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
43
44	!=======================================================================
45
46	! Auteur: P. Le Van /L. Fairhead/F.Hourdin
47	! -------
48
49	! Objet:
50	! ------
51
52	! GCM LMD nouvelle grille
53
54	!=======================================================================
55
56	! ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
57	! et possibilite d'appeler une fonction f(y) a derivee tangente
58	! hyperbolique a la place de la fonction a derivee sinusoidale.
59
60	! ... Possibilite de choisir le shema pour l'advection de
61	! q , en modifiant iadv dans traceur.def (10/02) .
62
63	! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
64	! Pour Van-Leer iadv=10
65
66	!-----------------------------------------------------------------------
67	! Declarations:
68	! -------------
69
70
71
72
73	REAL, INTENT(IN) :: time_0 ! not used
74
75	! dynamical variables:
76	REAL, INTENT(INOUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind
77	REAL, INTENT(INOUT) :: vcov(ip1jm, llm) ! meridional covariant wind
78	REAL, INTENT(INOUT) :: teta(ip1jmp1, llm) ! potential temperature
79	REAL, INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
80	REAL, INTENT(INOUT) :: masse(ip1jmp1, llm) ! air mass
81	REAL, INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
82	REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers
83
84	REAL :: p (ip1jmp1, llmp1) ! interlayer pressure
85	REAL :: pks(ip1jmp1) ! exner at the surface
86	REAL :: pk(ip1jmp1, llm) ! exner at mid-layer
87	REAL :: pkf(ip1jmp1, llm) ! filtered exner at mid-layer
88	REAL :: phi(ip1jmp1, llm) ! geopotential
89	REAL :: w(ip1jmp1, llm) ! vertical velocity
90
91	REAL :: zqmin, zqmax
92
93	! variables dynamiques intermediaire pour le transport
94	REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
95
96	! variables dynamiques au pas -1
97	REAL :: vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
98	REAL :: tetam1(ip1jmp1, llm), psm1(ip1jmp1)
99	REAL :: massem1(ip1jmp1, llm)
100
101	! tendances dynamiques
102	REAL :: dv(ip1jm, llm), du(ip1jmp1, llm)
103	REAL :: dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqtot), dp(ip1jmp1)
104
105	! tendances de la dissipation
106	REAL :: dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
107	REAL :: dtetadis(ip1jmp1, llm)
108
109	! tendances physiques
110	REAL :: dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
111	REAL :: dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqtot), dpfi(ip1jmp1)
112
113	! variables pour le fichier histoire
114	REAL :: dtav ! intervalle de temps elementaire
115
116	REAL :: tppn(iim), tpps(iim), tpn, tps
117
118	INTEGER :: itau, itaufinp1, iav
119	! INTEGER iday ! jour julien
120	REAL :: time
121
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	INTEGER :: itau_w ! pas de temps ecriture = itap + itau_phy
192
193	IF (nday>=0) THEN
194	itaufin = nday * day_step
195	else
196	itaufin = -nday
197	ENDIF
198	itaufinp1 = itaufin + 1
199	itau = 0
200	physic = .TRUE.
201	IF (iflag_phys==0.OR.iflag_phys==2) physic = .FALSE.
202
203	! iday = day_ini+itau/day_step
204	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
205	! IF(time.GT.1.) THEN
206	! time = time-1.
207	! iday = iday+1
208	! ENDIF
209
210
211	!-----------------------------------------------------------------------
212	! On initialise la pression et la fonction d'Exner :
213	! --------------------------------------------------
214
215	dq(:, :, :) = 0.
216	CALL pression (ip1jmp1, ap, bp, ps, p)
217	IF (pressure_exner) THEN
218	CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
219	else
220	CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
221	ENDIF
222
223	!-----------------------------------------------------------------------
224	! Debut de l'integration temporelle:
225	! ----------------------------------
226
227	1 CONTINUE ! Matsuno Forward step begins here
228
229	! date: (NB: date remains unchanged for Backward step)
230	! -----
231
232	jD_cur = jD_ref + day_ini - day_ref + &
233	(itau + 1) / day_step
234	jH_cur = jH_ref + start_time + &
235	mod(itau + 1, day_step) / float(day_step)
236	jD_cur = jD_cur + int(jH_cur)
237	jH_cur = jH_cur - int(jH_cur)
238
239	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 321')
240
241	IF (ok_guide) THEN
242	CALL guide_main(itau, ucov, vcov, teta, q, masse, ps)
243	ENDIF
244
245
246
247	! IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
248	! CALL test_period ( ucov,vcov,teta,q,p,phis )
249	! PRINT *,' ---- Test_period apres continue OK ! -----', itau
250	! ENDIF
251	!
252
253	! Save fields obtained at previous time step as '...m1'
254	CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)
255	CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)
256	CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)
257	CALL SCOPY(ijp1llm, masse, 1, massem1, 1)
258	CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)
259
260	forward = .TRUE.
261	leapf = .FALSE.
262	dt = dtvr
263
264	! ... P.Le Van .26/04/94 ....
265	! Ehouarn: finvmaold is actually not used
266	! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
267	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
268
269	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 400')
270
271	2 CONTINUE ! Matsuno backward or leapfrog step begins here
272
273	!-----------------------------------------------------------------------
274
275	! date: (NB: only leapfrog step requires recomputing date)
276	! -----
277
278	IF (leapf) THEN
279	jD_cur = jD_ref + day_ini - day_ref + &
280	(itau + 1) / day_step
281	jH_cur = jH_ref + start_time + &
282	mod(itau + 1, day_step) / float(day_step)
283	jD_cur = jD_cur + int(jH_cur)
284	jH_cur = jH_cur - int(jH_cur)
285	ENDIF
286
287
288	! gestion des appels de la physique et des dissipations:
289	! ------------------------------------------------------
290
291	! ... P.Le Van ( 6/02/95 ) ....
292
293	apphys = .FALSE.
294	statcl = .FALSE.
295	conser = .FALSE.
296	apdiss = .FALSE.
297
298	IF(purmats) THEN
299	! Purely Matsuno time stepping
300	IF(MOD(itau, iconser) ==0.AND. forward) conser = .TRUE.
301	IF(MOD(itau, dissip_period)==0.AND..NOT.forward) &
302	apdiss = .TRUE.
303	IF(MOD(itau, iphysiq)==0.AND..NOT.forward &
304	.AND. physic) apphys = .TRUE.
305	ELSE
306	! Leapfrog/Matsuno time stepping
307	IF(MOD(itau, iconser) == 0) conser = .TRUE.
308	IF(MOD(itau + 1, dissip_period)==0 .AND. .NOT. forward) &
309	apdiss = .TRUE.
310	IF(MOD(itau + 1, iphysiq)==0.AND.physic) apphys = .TRUE.
311	END IF
312
313	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
314	! supress dissipation step
315	IF (llm==1) THEN
316	apdiss = .FALSE.
317	ENDIF
318
319	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 589')
320
321	!-----------------------------------------------------------------------
322	! calcul des tendances dynamiques:
323	! --------------------------------
324
325	! compute geopotential phi()
326	CALL geopot (ip1jmp1, teta, pk, pks, phis, phi)
327
328	time = jD_cur + jH_cur
329	CALL caldyn &
330	(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
331	phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)
332
333
334	!-----------------------------------------------------------------------
335	! calcul des tendances advection des traceurs (dont l'humidite)
336	! -------------------------------------------------------------
337
338	CALL check_isotopes_seq(q, ip1jmp1, &
339	'leapfrog 686: avant caladvtrac')
340
341	IF(forward .OR. leapf) THEN
342	! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
343	CALL caladvtrac(q, pbaru, pbarv, &
344	p, masse, dq, teta, &
345	flxw, pk)
346	!WRITE(,) 'caladvtrac 346'
347
348	IF (offline) THEN
349	!maf stokage du flux de masse pour traceurs OFF-LINE
350
351	CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
352	dtvr, itau)
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==itaufin.AND..NOT.forward) lafin = .TRUE.
382	ELSE
383	IF(itau + 1 == 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 =="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>=1) THEN
436	ztit = 'bil dyn'
437	! Ehouarn: be careful, diagedyn is Earth-specific!
438	IF (planet_type=="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	IF (CPPKEY_PHYS) THEN
444	CALL calfis(lafin, jD_cur, jH_cur, &
445	ucov, vcov, teta, q, masse, ps, p, pk, phis, phi, &
446	du, dv, dteta, dq, &
447	flxw, dufi, dvfi, dtetafi, dqfi, dpfi)
448	END IF
449	! ajout des tendances physiques:
450	! ------------------------------
451	CALL addfi(dtphys, leapf, forward, &
452	ucov, vcov, teta, q, ps, &
453	dufi, dvfi, dtetafi, dqfi, dpfi)
454	! since addfi updates ps(), also update p(), masse() and pk()
455	CALL pression (ip1jmp1, ap, bp, ps, p)
456	CALL massdair(p, masse)
457	IF (pressure_exner) THEN
458	CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
459	else
460	CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
461	endif
462
463	IF (ok_strato) THEN
464	CALL top_bound(vcov, ucov, teta, masse, dtphys)
465	ENDIF
466
467
468	! Diagnostique de conservation de l'energie : difference
469	IF (ip_ebil_dyn>=1) THEN
470	ztit = 'bil phys'
471	IF (planet_type=="earth") THEN
472	CALL diagedyn(ztit, 2, 1, 1, dtphys &
473	, ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
474	ENDIF
475	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
476
477	ENDIF ! of IF( apphys )
478
479	IF(iflag_phys==2) THEN ! "Newtonian" case
480	! Academic case : Simple friction and Newtonan relaxation
481	! -------------------------------------------------------
482	DO l = 1, llm
483	DO ij = 1, ip1jmp1
484	teta(ij, l) = teta(ij, l) - dtvr * &
485	(teta(ij, l) - tetarappel(ij, l)) * (knewt_g + knewt_t(l) * clat4(ij))
486	ENDDO
487	ENDDO ! of DO l=1,llm
488
489	IF (planet_type=="giant") THEN
490	! add an intrinsic heat flux at the base of the atmosphere
491	teta(:, 1) = teta(:, 1) + dtvr * aire(:) * ihf / cpp / masse(:, 1)
492	endif
493
494	CALL friction(ucov, vcov, dtvr)
495
496	! Sponge layer (if any)
497	IF (ok_strato) THEN
498	! dufi(:,:)=0.
499	! dvfi(:,:)=0.
500	! dtetafi(:,:)=0.
501	! dqfi(:,:,:)=0.
502	! dpfi(:)=0.
503	! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
504	CALL top_bound(vcov, ucov, teta, masse, dtvr)
505	! CALL addfi( dtvr, leapf, forward ,
506	! $ ucov, vcov, teta , q ,ps ,
507	! $ dufi, dvfi, dtetafi , dqfi ,dpfi )
508	ENDIF ! of IF (ok_strato)
509	ENDIF ! of IF (iflag_phys.EQ.2)
510
511
512	!-jld
513
514	CALL pression (ip1jmp1, ap, bp, ps, p)
515	IF (pressure_exner) THEN
516	CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
517	else
518	CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
519	ENDIF
520	CALL massdair(p, masse)
521
522	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1196')
523
524	!-----------------------------------------------------------------------
525	! dissipation horizontale et verticale des petites echelles:
526	! ----------------------------------------------------------
527
528	IF(apdiss) THEN
529
530
531	! calcul de l'energie cinetique avant dissipation
532	CALL covcont(llm, ucov, vcov, ucont, vcont)
533	CALL enercin(vcov, ucov, vcont, ucont, ecin0)
534
535	! dissipation
536	CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)
537	ucov = ucov + dudis
538	vcov = vcov + dvdis
539	! teta=teta+dtetadis
540
541
542	!------------------------------------------------------------------------
543	IF (dissip_conservative) THEN
544	! On rajoute la tendance due a la transform. Ec -> E therm. cree
545	! lors de la dissipation
546	CALL covcont(llm, ucov, vcov, ucont, vcont)
547	CALL enercin(vcov, ucov, vcont, ucont, ecin)
548	dtetaecdt = (ecin0 - ecin) / pk
549	! teta=teta+dtetaecdt
550	dtetadis = dtetadis + dtetaecdt
551	endif
552	teta = teta + dtetadis
553	!------------------------------------------------------------------------
554
555
556	! ....... P. Le Van ( ajout le 17/04/96 ) ...........
557	! ... Calcul de la valeur moyenne, unique de h aux poles .....
558	!
559
560	DO l = 1, llm
561	DO ij = 1, iim
562	tppn(ij) = aire(ij) * teta(ij, l)
563	tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
564	ENDDO
565	tpn = SSUM(iim, tppn, 1) / apoln
566	tps = SSUM(iim, tpps, 1) / apols
567
568	DO ij = 1, iip1
569	teta(ij, l) = tpn
570	teta(ij + ip1jm, l) = tps
571	ENDDO
572	ENDDO
573
574	IF (1 == 0) THEN
575	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
576	!!! 2) should probably not be here anyway
577	!!! but are kept for those who would want to revert to previous behaviour
578	DO ij = 1, iim
579	tppn(ij) = aire(ij) * ps (ij)
580	tpps(ij) = aire(ij + ip1jm) * ps (ij + ip1jm)
581	ENDDO
582	tpn = SSUM(iim, tppn, 1) / apoln
583	tps = SSUM(iim, tpps, 1) / apols
584
585	DO ij = 1, iip1
586	ps(ij) = tpn
587	ps(ij + ip1jm) = tps
588	ENDDO
589	endif ! of if (1 == 0)
590
591	END IF ! of IF(apdiss)
592
593	! ajout debug
594	! IF( lafin ) THEN
595	! abort_message = 'Simulation finished'
596	! CALL abort_gcm(modname,abort_message,0)
597	! ENDIF
598
599	! ********************************************************************
600	! ********************************************************************
601	! .... fin de l'integration dynamique et physique pour le pas itau ..
602	! ********************************************************************
603	! ********************************************************************
604
605	! preparation du pas d'integration suivant ......
606
607	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1509')
608
609	IF (.NOT.purmats) THEN
610	! ........................................................
611	! .............. schema matsuno + leapfrog ..............
612	! ........................................................
613
614	IF(forward .OR. leapf) THEN
615	itau = itau + 1
616	! iday= day_ini+itau/day_step
617	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
618	! IF(time.GT.1.) THEN
619	! time = time-1.
620	! iday = iday+1
621	! ENDIF
622	ENDIF
623
624	IF(itau == itaufinp1) THEN
625	IF (flag_verif) THEN
626	WRITE(79, *) 'ucov', ucov
627	WRITE(80, *) 'vcov', vcov
628	WRITE(81, *) 'teta', teta
629	WRITE(82, *) 'ps', ps
630	WRITE(83, *) 'q', q
631	WRITE(85, *) 'q1 = ', q(:, :, 1)
632	WRITE(86, *) 'q3 = ', q(:, :, 3)
633	endif
634
635	abort_message = 'Simulation finished'
636
637	CALL abort_gcm(modname, abort_message, 0)
638	ENDIF
639	!-----------------------------------------------------------------------
640	! ecriture du fichier histoire moyenne:
641	! -------------------------------------
642
643	IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
644	IF(itau==itaufin) THEN
645	iav = 1
646	ELSE
647	iav = 0
648	ENDIF
649
650	! Ehouarn: re-compute geopotential for outputs
651	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
652
653	IF (ok_dynzon) THEN
654	CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
655	ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
656	END IF
657	IF (ok_dyn_ave) THEN
658	CALL writedynav(itau, vcov, &
659	ucov, teta, pk, phi, q, masse, ps, phis)
660	ENDIF
661
662	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
663
664	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1584')
665
666	!-----------------------------------------------------------------------
667	! ecriture de la bande histoire:
668	! ------------------------------
669
670	IF(MOD(itau, iecri)==0) THEN
671	! ! Ehouarn: output only during LF or Backward Matsuno
672	IF (leapf.OR.(.NOT.leapf.AND.(.NOT.forward))) THEN
673	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
674	unat = 0.
675	DO l = 1, llm
676	unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
677	vnat(:, l) = vcov(:, l) / cv(:)
678	enddo
679	IF (ok_dyn_ins) THEN
680	! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
681	CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
682	! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
683	! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
684	! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
685	! CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
686	! CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
687	endif ! of if (ok_dyn_ins)
688	! For some Grads outputs of fields
689	IF (output_grads_dyn) THEN
690	INCLUDE "write_grads_dyn.h"
691	endif
692	endif ! of if (leapf.OR.(.NOT.leapf.AND.(.NOT.forward)))
693	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
694
695	IF(itau==itaufin) THEN
696
697
698	! if (planet_type.EQ."earth") THEN
699	! Write an Earth-format restart file
700	CALL dynredem1("restart.nc", start_time, &
701	vcov, ucov, teta, q, masse, ps)
702	! END IF ! of if (planet_type.EQ."earth")
703
704	CLOSE(99)
705	IF (ok_guide) THEN
706	! ! set ok_guide to false to avoid extra output
707	! ! in following forward step
708	ok_guide = .FALSE.
709	endif
710	! !!! Ehouarn: Why not stop here and now?
711	ENDIF ! of IF (itau.EQ.itaufin)
712
713	!-----------------------------------------------------------------------
714	! gestion de l'integration temporelle:
715	! ------------------------------------
716
717	IF(MOD(itau, iperiod)==0) THEN
718	GO TO 1
719	ELSE IF (MOD(itau - 1, iperiod) == 0) THEN
720
721	IF(forward) THEN
722	! fin du pas forward et debut du pas backward
723
724	forward = .FALSE.
725	leapf = .FALSE.
726	GO TO 2
727
728	ELSE
729	! fin du pas backward et debut du premier pas leapfrog
730
731	leapf = .TRUE.
732	dt = 2. * dtvr
733	GO TO 2
734	END IF ! of IF (forward)
735	ELSE
736
737	! ...... pas leapfrog .....
738
739	leapf = .TRUE.
740	dt = 2. * dtvr
741	GO TO 2
742	END IF ! of IF (MOD(itau,iperiod).EQ.0)
743	! ! ELSEIF (MOD(itau-1,iperiod).EQ.0)
744
745	ELSE ! of IF (.NOT.purmats)
746
747	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1664')
748
749	! ........................................................
750	! .............. schema matsuno ...............
751	! ........................................................
752	IF(forward) THEN
753
754	itau = itau + 1
755	! iday = day_ini+itau/day_step
756	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
757
758	! IF(time.GT.1.) THEN
759	! time = time-1.
760	! iday = iday+1
761	! ENDIF
762
763	forward = .FALSE.
764	IF(itau == itaufinp1) THEN
765	abort_message = 'Simulation finished'
766	CALL abort_gcm(modname, abort_message, 0)
767	ENDIF
768	GO TO 2
769
770	ELSE ! of IF(forward) i.e. backward step
771
772	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1698')
773
774	IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
775	IF(itau==itaufin) THEN
776	iav = 1
777	ELSE
778	iav = 0
779	ENDIF
780
781	! ! Ehouarn: re-compute geopotential for outputs
782	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
783
784	IF (ok_dynzon) THEN
785	CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
786	ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
787	ENDIF
788	IF (ok_dyn_ave) THEN
789	CALL writedynav(itau, vcov, &
790	ucov, teta, pk, phi, q, masse, ps, phis)
791	ENDIF
792
793	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
794
795	IF(MOD(itau, iecri)==0) THEN
796	! IF(MOD(itau,iecri*day_step).EQ.0) THEN
797	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
798	unat = 0.
799	DO l = 1, llm
800	unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
801	vnat(:, l) = vcov(:, l) / cv(:)
802	enddo
803	IF (ok_dyn_ins) THEN
804	! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
805	! & itau,iecri
806	CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
807	endif ! of if (ok_dyn_ins)
808	! For some Grads outputs
809	IF (output_grads_dyn) THEN
810	INCLUDE "write_grads_dyn.h"
811	endif
812
813	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
814
815	IF(itau==itaufin) THEN
816	! if (planet_type.EQ."earth") THEN
817	CALL dynredem1("restart.nc", start_time, &
818	vcov, ucov, teta, q, masse, ps)
819	! END IF ! of if (planet_type.EQ."earth")
820	IF (ok_guide) THEN
821	! ! set ok_guide to false to avoid extra output
822	! ! in following forward step
823	ok_guide = .FALSE.
824	endif
825	ENDIF ! of IF(itau.EQ.itaufin)
826
827	forward = .TRUE.
828	GO TO 1
829
830	ENDIF ! of IF (forward)
831
832	END IF ! of IF(.NOT.purmats)
833
834	END SUBROUTINE leapfrog

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

Download in other formats:

Original Format