Context Navigation

lmdz_leapfrog.f90 @ 5214

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

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

Download in other formats:

Original Format