Context Navigation

leapfrog.F @ 1578

Last change on this file since 1578 was 1577, checked in by Laurent Fairhead, 13 years ago

Modifications au code qui permettent de commencer une simulation à n'importe
quelle heure de la journée. On fait toujours un nombre entier de jours de
simulation.
On spécifie cette heure de départ dans la variable starttime du run.def (la
valeur est en jour et elle est à zéro par défaut).
La valeur est sauvegardée dans le fichier restart.nc. Les valeurs lues dans
le fichier start et le run.def sont comparées en début de simulation. La
simulation s'arrête si elles ne sont pas égales sauf si une remise à zéro de
la date a été demandée.
Par ailleurs, la fréquence de lecture des conditions aux limites a été modifiée
pour qu'à chaque changement de jour, celles-ci soient mises à jour (jusqu'à
maintenant elles étaient mises à jour à une fréquence donnée qui, en cas de
départ de simulation à une heure différente de minuit, ne correspondait pas
forcèment à un changement dans la date).
Validation effectuée en traçant le flux solaire descendant au sommet de
l'atmosphère à différentes heures de la journée, après un redémarrage, en
s'assurant que le maximum est bien là où il est sensé être.

Modifications to the code to enable it to be started at any time of the day.
The code still runs for an integer number of days.
The start time is specified using variable starttime in the run.def file (the
value is in days and is zero by default).
The start time is saved in the restart.nc file at the end of the simulation.
The values read in from the start.nc file and the run.def file are compared
at the start of the simulation. If they differ, the simulation is aborted
unless the raz_date variable has been set.
Furthermore, the frequency at which boundary conditions are read in has been
modified so that they are updated everyday at midnight (until now, they were
updated at a certain frequency that, in case of a simulation starting at a time
other than midnight, did not ensure that those conditions would be updated each
day at midnight)
The modifications were validated by plotting the downward solaf flux at TOA at
different times of the day (and after having restarted the simulation) and
ensuring that the maximum of flux was at the right place according to local
time.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 24.8 KB

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

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

Download in other formats:

Original Format