Context Navigation

lmdz_leapfrog.f90 @ 5186

Last change on this file since 5186 was 5186, checked in by abarral, 2 months ago
Encapsulate files in 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: 27.5 KB

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

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

Download in other formats:

Original Format