Context Navigation

leapfrog.F @ 4350

Last change on this file since 4350 was 3852, checked in by dcugnet, 3 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

the number of tracers on the first line
one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

Format 2: read the sections from the "tracer.def" file.
format 3: read one section each "tracer_<section name>.def" file.
the first line of a section is "&<section name>
the other lines start with a tracer name followed by <key>=<val> pairs.
the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".

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.7 KB

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

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

Download in other formats:

Original Format