Context Navigation

leapfrog.F90 @ 5127

Last change on this file since 5127 was 5123, checked in by abarral, 3 months ago
Correct various minor mistakes from previous commits
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 5123 2024-07-25 06:45:50Z snguyen $
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
31	IMPLICIT NONE
32
33	! ...... Version du 10/01/98 ..........
34
35	! avec coordonnees verticales hybrides
36	! avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
37
38	!=======================================================================
39	!
40	! Auteur: P. Le Van /L. Fairhead/F.Hourdin
41	! -------
42	!
43	! Objet:
44	! ------
45	!
46	! GCM LMD nouvelle grille
47	!
48	!=======================================================================
49	!
50	! ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
51	! et possibilite d'appeler une fonction f(y) a derivee tangente
52	! hyperbolique a la place de la fonction a derivee sinusoidale.
53
54	! ... Possibilite de choisir le shema pour l'advection de
55	! q , en modifiant iadv dans traceur.def (10/02) .
56	!
57	! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
58	! Pour Van-Leer iadv=10
59	!
60	!-----------------------------------------------------------------------
61	! Declarations:
62	! -------------
63
64	include "dimensions.h"
65	include "paramet.h"
66	include "comdissnew.h"
67	include "comgeom.h"
68	include "academic.h"
69
70	REAL, INTENT(IN) :: time_0 ! not used
71
72	! dynamical variables:
73	REAL, INTENT(INOUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind
74	REAL, INTENT(INOUT) :: vcov(ip1jm, llm) ! meridional covariant wind
75	REAL, INTENT(INOUT) :: teta(ip1jmp1, llm) ! potential temperature
76	REAL, INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
77	REAL, INTENT(INOUT) :: masse(ip1jmp1, llm) ! air mass
78	REAL, INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
79	REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers
80
81	REAL :: p (ip1jmp1, llmp1) ! interlayer pressure
82	REAL :: pks(ip1jmp1) ! exner at the surface
83	REAL :: pk(ip1jmp1, llm) ! exner at mid-layer
84	REAL :: pkf(ip1jmp1, llm) ! filtered exner at mid-layer
85	REAL :: phi(ip1jmp1, llm) ! geopotential
86	REAL :: w(ip1jmp1, llm) ! vertical velocity
87
88	REAL :: zqmin, zqmax
89
90	! variables dynamiques intermediaire pour le transport
91	REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
92
93	! variables dynamiques au pas -1
94	REAL :: vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
95	REAL :: tetam1(ip1jmp1, llm), psm1(ip1jmp1)
96	REAL :: massem1(ip1jmp1, llm)
97
98	! tendances dynamiques
99	REAL :: dv(ip1jm, llm), du(ip1jmp1, llm)
100	REAL :: dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqtot), dp(ip1jmp1)
101
102	! tendances de la dissipation
103	REAL :: dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
104	REAL :: dtetadis(ip1jmp1, llm)
105
106	! tendances physiques
107	REAL :: dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
108	REAL :: dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqtot), dpfi(ip1jmp1)
109
110	! variables pour le fichier histoire
111	REAL :: dtav ! intervalle de temps elementaire
112
113	REAL :: tppn(iim), tpps(iim), tpn, tps
114	!
115	INTEGER :: itau, itaufinp1, iav
116	! INTEGER iday ! jour julien
117	REAL :: time
118
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	ENDIF ! of IF (offline)
352	!
353	ENDIF ! of IF( forward .OR. leapf )
354
355
356	!-----------------------------------------------------------------------
357	! integrations dynamique et traceurs:
358	! ----------------------------------
359
360	CALL msg('720', modname, isoCheck)
361	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 756')
362
363	CALL integrd (nqtot, vcovm1, ucovm1, tetam1, psm1, massem1, &
364	dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis)
365	! $ finvmaold )
366
367	CALL msg('724', modname, isoCheck)
368	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 762')
369
370	! .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
371	!
372	!-----------------------------------------------------------------------
373	! calcul des tendances physiques:
374	! -------------------------------
375	! ######## P.Le Van ( Modif le 6/02/95 ) ###########
376	!
377	IF(purmats) THEN
378	IF(itau==itaufin.AND..NOT.forward) lafin = .TRUE.
379	ELSE
380	IF(itau + 1 == itaufin) lafin = .TRUE.
381	ENDIF
382	!
383	!
384	IF(apphys) THEN
385	!
386	! ....... Ajout P.Le Van ( 17/04/96 ) ...........
387	!
388
389	CALL pression (ip1jmp1, ap, bp, ps, p)
390	IF (pressure_exner) THEN
391	CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
392	else
393	CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
394	endif
395
396	! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
397	! avec dyn3dmem
398	CALL geopot (ip1jmp1, teta, pk, pks, phis, phi)
399
400	! rdaym_ini = itau * dtvr / daysec
401	! rdayvrai = rdaym_ini + day_ini
402	! jD_cur = jD_ref + day_ini - day_ref
403	! $ + int (itau * dtvr / daysec)
404	! jH_cur = jH_ref + &
405	! & (itau * dtvr / daysec - int(itau * dtvr / daysec))
406	jD_cur = jD_ref + day_ini - day_ref + &
407	(itau + 1) / day_step
408
409	IF (planet_type =="generic") THEN
410	! AS: we make jD_cur to be pday
411	jD_cur = int(day_ini + itau / day_step)
412	ENDIF
413
414	jH_cur = jH_ref + start_time + &
415	mod(itau + 1, day_step) / float(day_step)
416	jD_cur = jD_cur + int(jH_cur)
417	jH_cur = jH_cur - int(jH_cur)
418	! WRITE(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
419	! CALL ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
420	! WRITE(lunout,*)'current date = ',an, mois, jour, secondes
421
422	! rajout debug
423	! lafin = .TRUE.
424
425
426	! Inbterface avec les routines de phylmd (phymars ... )
427	! -----------------------------------------------------
428
429	!+jld
430
431	! Diagnostique de conservation de l'energie : initialisation
432	IF (ip_ebil_dyn>=1) THEN
433	ztit = 'bil dyn'
434	! Ehouarn: be careful, diagedyn is Earth-specific!
435	IF (planet_type=="earth") THEN
436	CALL diagedyn(ztit, 2, 1, 1, dtphys &
437	, ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
438	ENDIF
439	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
440	IF (CPPKEY_PHYS) THEN
441	CALL calfis(lafin, jD_cur, jH_cur, &
442	ucov, vcov, teta, q, masse, ps, p, pk, phis, phi, &
443	du, dv, dteta, dq, &
444	flxw, dufi, dvfi, dtetafi, dqfi, dpfi)
445	END IF
446	! ajout des tendances physiques:
447	! ------------------------------
448	CALL addfi(dtphys, leapf, forward, &
449	ucov, vcov, teta, q, ps, &
450	dufi, dvfi, dtetafi, dqfi, dpfi)
451	! since addfi updates ps(), also update p(), masse() and pk()
452	CALL pression (ip1jmp1, ap, bp, ps, p)
453	CALL massdair(p, masse)
454	IF (pressure_exner) THEN
455	CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
456	else
457	CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
458	endif
459
460	IF (ok_strato) THEN
461	CALL top_bound(vcov, ucov, teta, masse, dtphys)
462	ENDIF
463
464	!
465	! Diagnostique de conservation de l'energie : difference
466	IF (ip_ebil_dyn>=1) THEN
467	ztit = 'bil phys'
468	IF (planet_type=="earth") THEN
469	CALL diagedyn(ztit, 2, 1, 1, dtphys &
470	, ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
471	ENDIF
472	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
473
474	ENDIF ! of IF( apphys )
475
476	IF(iflag_phys==2) THEN ! "Newtonian" case
477	! Academic case : Simple friction and Newtonan relaxation
478	! -------------------------------------------------------
479	DO l = 1, llm
480	DO ij = 1, ip1jmp1
481	teta(ij, l) = teta(ij, l) - dtvr * &
482	(teta(ij, l) - tetarappel(ij, l)) * (knewt_g + knewt_t(l) * clat4(ij))
483	ENDDO
484	ENDDO ! of DO l=1,llm
485
486	IF (planet_type=="giant") THEN
487	! add an intrinsic heat flux at the base of the atmosphere
488	teta(:, 1) = teta(:, 1) + dtvr * aire(:) * ihf / cpp / masse(:, 1)
489	endif
490
491	CALL friction(ucov, vcov, dtvr)
492
493	! Sponge layer (if any)
494	IF (ok_strato) THEN
495	! dufi(:,:)=0.
496	! dvfi(:,:)=0.
497	! dtetafi(:,:)=0.
498	! dqfi(:,:,:)=0.
499	! dpfi(:)=0.
500	! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
501	CALL top_bound(vcov, ucov, teta, masse, dtvr)
502	! CALL addfi( dtvr, leapf, forward ,
503	! $ ucov, vcov, teta , q ,ps ,
504	! $ dufi, dvfi, dtetafi , dqfi ,dpfi )
505	ENDIF ! of IF (ok_strato)
506	ENDIF ! of IF (iflag_phys.EQ.2)
507
508
509	!-jld
510
511	CALL pression (ip1jmp1, ap, bp, ps, p)
512	IF (pressure_exner) THEN
513	CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
514	else
515	CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
516	ENDIF
517	CALL massdair(p, masse)
518
519	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1196')
520
521	!-----------------------------------------------------------------------
522	! dissipation horizontale et verticale des petites echelles:
523	! ----------------------------------------------------------
524
525	IF(apdiss) THEN
526
527
528	! calcul de l'energie cinetique avant dissipation
529	CALL covcont(llm, ucov, vcov, ucont, vcont)
530	CALL enercin(vcov, ucov, vcont, ucont, ecin0)
531
532	! dissipation
533	CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)
534	ucov = ucov + dudis
535	vcov = vcov + dvdis
536	! teta=teta+dtetadis
537
538
539	!------------------------------------------------------------------------
540	IF (dissip_conservative) THEN
541	! On rajoute la tendance due a la transform. Ec -> E therm. cree
542	! lors de la dissipation
543	CALL covcont(llm, ucov, vcov, ucont, vcont)
544	CALL enercin(vcov, ucov, vcont, ucont, ecin)
545	dtetaecdt = (ecin0 - ecin) / pk
546	! teta=teta+dtetaecdt
547	dtetadis = dtetadis + dtetaecdt
548	endif
549	teta = teta + dtetadis
550	!------------------------------------------------------------------------
551
552
553	! ....... P. Le Van ( ajout le 17/04/96 ) ...........
554	! ... Calcul de la valeur moyenne, unique de h aux poles .....
555	!
556
557	DO l = 1, llm
558	DO ij = 1, iim
559	tppn(ij) = aire(ij) * teta(ij, l)
560	tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
561	ENDDO
562	tpn = SSUM(iim, tppn, 1) / apoln
563	tps = SSUM(iim, tpps, 1) / apols
564
565	DO ij = 1, iip1
566	teta(ij, l) = tpn
567	teta(ij + ip1jm, l) = tps
568	ENDDO
569	ENDDO
570
571	IF (1 == 0) THEN
572	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
573	!!! 2) should probably not be here anyway
574	!!! but are kept for those who would want to revert to previous behaviour
575	DO ij = 1, iim
576	tppn(ij) = aire(ij) * ps (ij)
577	tpps(ij) = aire(ij + ip1jm) * ps (ij + ip1jm)
578	ENDDO
579	tpn = SSUM(iim, tppn, 1) / apoln
580	tps = SSUM(iim, tpps, 1) / apols
581
582	DO ij = 1, iip1
583	ps(ij) = tpn
584	ps(ij + ip1jm) = tps
585	ENDDO
586	endif ! of if (1 == 0)
587
588	END IF ! of IF(apdiss)
589
590	! ajout debug
591	! IF( lafin ) THEN
592	! abort_message = 'Simulation finished'
593	! CALL abort_gcm(modname,abort_message,0)
594	! ENDIF
595
596	! ********************************************************************
597	! ********************************************************************
598	! .... fin de l'integration dynamique et physique pour le pas itau ..
599	! ********************************************************************
600	! ********************************************************************
601
602	! preparation du pas d'integration suivant ......
603
604	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1509')
605
606	IF (.NOT.purmats) THEN
607	! ........................................................
608	! .............. schema matsuno + leapfrog ..............
609	! ........................................................
610
611	IF(forward .OR. leapf) THEN
612	itau = itau + 1
613	! iday= day_ini+itau/day_step
614	! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
615	! IF(time.GT.1.) THEN
616	! time = time-1.
617	! iday = iday+1
618	! ENDIF
619	ENDIF
620
621	IF(itau == itaufinp1) THEN
622	IF (flag_verif) THEN
623	WRITE(79, *) 'ucov', ucov
624	WRITE(80, *) 'vcov', vcov
625	WRITE(81, *) 'teta', teta
626	WRITE(82, *) 'ps', ps
627	WRITE(83, *) 'q', q
628	WRITE(85, *) 'q1 = ', q(:, :, 1)
629	WRITE(86, *) 'q3 = ', q(:, :, 3)
630	endif
631
632	abort_message = 'Simulation finished'
633
634	CALL abort_gcm(modname, abort_message, 0)
635	ENDIF
636	!-----------------------------------------------------------------------
637	! ecriture du fichier histoire moyenne:
638	! -------------------------------------
639
640	IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
641	IF(itau==itaufin) THEN
642	iav = 1
643	ELSE
644	iav = 0
645	ENDIF
646
647	! Ehouarn: re-compute geopotential for outputs
648	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
649
650	IF (ok_dynzon) THEN
651	CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
652	ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
653	END IF
654	IF (ok_dyn_ave) THEN
655	CALL writedynav(itau, vcov, &
656	ucov, teta, pk, phi, q, masse, ps, phis)
657	ENDIF
658
659	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
660
661	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1584')
662
663	!-----------------------------------------------------------------------
664	! ecriture de la bande histoire:
665	! ------------------------------
666
667	IF(MOD(itau, iecri)==0) THEN
668	! ! Ehouarn: output only during LF or Backward Matsuno
669	IF (leapf.OR.(.NOT.leapf.AND.(.NOT.forward))) THEN
670	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
671	unat = 0.
672	do l = 1, llm
673	unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
674	vnat(:, l) = vcov(:, l) / cv(:)
675	enddo
676	IF (ok_dyn_ins) THEN
677	! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
678	CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
679	! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
680	! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
681	! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
682	! CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
683	! CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
684	endif ! of if (ok_dyn_ins)
685	! For some Grads outputs of fields
686	IF (output_grads_dyn) THEN
687	include "write_grads_dyn.h"
688	endif
689	endif ! of if (leapf.OR.(.NOT.leapf.AND.(.NOT.forward)))
690	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
691
692	IF(itau==itaufin) THEN
693
694
695	! if (planet_type.EQ."earth") THEN
696	! Write an Earth-format restart file
697	CALL dynredem1("restart.nc", start_time, &
698	vcov, ucov, teta, q, masse, ps)
699	! END IF ! of if (planet_type.EQ."earth")
700
701	CLOSE(99)
702	IF (ok_guide) THEN
703	! ! set ok_guide to false to avoid extra output
704	! ! in following forward step
705	ok_guide = .FALSE.
706	endif
707	! !!! Ehouarn: Why not stop here and now?
708	ENDIF ! of IF (itau.EQ.itaufin)
709
710	!-----------------------------------------------------------------------
711	! gestion de l'integration temporelle:
712	! ------------------------------------
713
714	IF(MOD(itau, iperiod)==0) THEN
715	GO TO 1
716	ELSE IF (MOD(itau - 1, iperiod) == 0) THEN
717
718	IF(forward) THEN
719	! fin du pas forward et debut du pas backward
720
721	forward = .FALSE.
722	leapf = .FALSE.
723	GO TO 2
724
725	ELSE
726	! fin du pas backward et debut du premier pas leapfrog
727
728	leapf = .TRUE.
729	dt = 2. * dtvr
730	GO TO 2
731	END IF ! of IF (forward)
732	ELSE
733
734	! ...... pas leapfrog .....
735
736	leapf = .TRUE.
737	dt = 2. * dtvr
738	GO TO 2
739	END IF ! of IF (MOD(itau,iperiod).EQ.0)
740	! ! ELSEIF (MOD(itau-1,iperiod).EQ.0)
741
742	ELSE ! of IF (.NOT.purmats)
743
744	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1664')
745
746	! ........................................................
747	! .............. schema matsuno ...............
748	! ........................................................
749	IF(forward) THEN
750
751	itau = itau + 1
752	! iday = day_ini+itau/day_step
753	! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
754	!
755	! IF(time.GT.1.) THEN
756	! time = time-1.
757	! iday = iday+1
758	! ENDIF
759
760	forward = .FALSE.
761	IF(itau == itaufinp1) THEN
762	abort_message = 'Simulation finished'
763	CALL abort_gcm(modname, abort_message, 0)
764	ENDIF
765	GO TO 2
766
767	ELSE ! of IF(forward) i.e. backward step
768
769	CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1698')
770
771	IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
772	IF(itau==itaufin) THEN
773	iav = 1
774	ELSE
775	iav = 0
776	ENDIF
777
778	! ! Ehouarn: re-compute geopotential for outputs
779	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
780
781	IF (ok_dynzon) THEN
782	CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
783	ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
784	ENDIF
785	IF (ok_dyn_ave) THEN
786	CALL writedynav(itau, vcov, &
787	ucov, teta, pk, phi, q, masse, ps, phis)
788	ENDIF
789
790	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
791
792	IF(MOD(itau, iecri)==0) THEN
793	! IF(MOD(itau,iecri*day_step).EQ.0) THEN
794	CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
795	unat = 0.
796	do l = 1, llm
797	unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
798	vnat(:, l) = vcov(:, l) / cv(:)
799	enddo
800	IF (ok_dyn_ins) THEN
801	! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
802	! & itau,iecri
803	CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
804	endif ! of if (ok_dyn_ins)
805	! For some Grads outputs
806	IF (output_grads_dyn) THEN
807	include "write_grads_dyn.h"
808	endif
809
810	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
811
812	IF(itau==itaufin) THEN
813	! if (planet_type.EQ."earth") THEN
814	CALL dynredem1("restart.nc", start_time, &
815	vcov, ucov, teta, q, masse, ps)
816	! END IF ! of if (planet_type.EQ."earth")
817	IF (ok_guide) THEN
818	! ! set ok_guide to false to avoid extra output
819	! ! in following forward step
820	ok_guide = .FALSE.
821	endif
822	ENDIF ! of IF(itau.EQ.itaufin)
823
824	forward = .TRUE.
825	GO TO 1
826
827	ENDIF ! of IF (forward)
828
829	END IF ! of IF(.NOT.purmats)
830
831	END SUBROUTINE leapfrog

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

Download in other formats:

Original Format