Context Navigation

leapfrog.F @ 2037

Last change on this file since 2037 was 2021, checked in by lguez, 10 years ago

Removed unused variables pks, pk, pkf from main program unit gcm.

Encapsulated procedures exner_hyb, exner_hyb_p, exner_hyb_loc,
exner_milieu, exner_milieu_p and exner_milieu_loc into
modules. (Compulsory to allow optional arguments.)

In the procedures exner_hyb, exner_hyb_p, exner_hyb_loc, donwgraded
arguments alpha and beta to local variables. In exner_milieu,
exner_milieu_p and exner_milieu_loc, removed beta altogether. In the
six procedures exner_*, made pkf an optional argument. Made some
cosmetic modifications in order to keep the six procedures exner_* as
close as possible.

In the six procedures exner_*, removed the averaging of pks at the
poles: this is not useful because ps is already the same at all
longitudes at the poles. This modification changes the results of the
program. Motivation: had to do this for exner_hyb because we call it
from test_disvert with a few surface pressure values.

In all the procedures calling exner_*, removed the variables alpha and
beta. Also removed variables alpha and beta from module leapfrog_mod
and from module call_calfis_mod.

Removed actual argument pkf in call to exner_hyb* and exner_milieu*
from guide_interp, guide_main, iniacademic and iniacademic_loc (pkf
was not used in those procedures).

Argument workvar of startget_dyn is used only if varname is tpot or

When varname is tpot or q, the actual argument associated to

workvar in etat0_netcdf is not y. So y in etat0_netcdf is a
place-holder, never used. So we remove optional argument y in the
calls to exner_hyb and exner_milieu from etat0_netcdf.

Created procedure test_disvert, called only by etat0_netcdf. This
procedure tests the order of pressure values at half-levels and full
levels.

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: 26.1 KB

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

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

Download in other formats:

Original Format