Context Navigation

leapfrog.F @ 2019

Last change on this file since 2019 was 1987, checked in by Ehouarn Millour, 12 years ago
Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency). Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid). EM
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.2 KB

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

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

Download in other formats:

Original Format