Context Navigation

leapfrog.F @ 3084

Last change on this file since 3084 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

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

File size: 23.5 KB

Line
1	!
2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z oboucher $
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
152	data callinigrads/.true./
153	character*10 string10
154
155	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
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=20) :: modname
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	itaufin = nday*day_step
194	itaufinp1 = itaufin +1
195	modname="leapfrog"
196
197
198	itau = 0
199	c iday = day_ini+itau/day_step
200	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
201	c IF(time.GT.1.) THEN
202	c time = time-1.
203	c iday = iday+1
204	c ENDIF
205
206
207	c-----------------------------------------------------------------------
208	c On initialise la pression et la fonction d'Exner :
209	c --------------------------------------------------
210
211	dq(:,:,:)=0.
212	CALL pression ( ip1jmp1, ap, bp, ps, p )
213	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
214
215	c-----------------------------------------------------------------------
216	c Debut de l'integration temporelle:
217	c ----------------------------------
218
219	1 CONTINUE
220
221	jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
222	jH_cur = jH_ref + &
223	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
224
225
226	#ifdef CPP_IOIPSL
227	if (ok_guide) then
228	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
229	endif
230	#endif
231
232
233	c
234	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
235	c CALL test_period ( ucov,vcov,teta,q,p,phis )
236	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
237	c ENDIF
238	c
239
240	! Save fields obtained at previous time step as '...m1'
241	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
242	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
243	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
244	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
245	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
246
247	forward = .TRUE.
248	leapf = .FALSE.
249	dt = dtvr
250
251	c ... P.Le Van .26/04/94 ....
252
253	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
254	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
255
256	2 CONTINUE
257
258	c-----------------------------------------------------------------------
259
260	c date:
261	c -----
262
263
264	c gestion des appels de la physique et des dissipations:
265	c ------------------------------------------------------
266	c
267	c ... P.Le Van ( 6/02/95 ) ....
268
269	apphys = .FALSE.
270	statcl = .FALSE.
271	conser = .FALSE.
272	apdiss = .FALSE.
273
274	IF( purmats ) THEN
275	! Purely Matsuno time stepping
276	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
277	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
278	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
279	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
280	ELSE
281	! Leapfrog/Matsuno time stepping
282	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
283	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
284	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
285	END IF
286
287	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
288	! supress dissipation step
289	if (llm.eq.1) then
290	apdiss=.false.
291	endif
292
293	c-----------------------------------------------------------------------
294	c calcul des tendances dynamiques:
295	c --------------------------------
296
297	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
298
299	time = jD_cur + jH_cur
300	CALL caldyn
301	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
302	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
303
304
305	c-----------------------------------------------------------------------
306	c calcul des tendances advection des traceurs (dont l'humidite)
307	c -------------------------------------------------------------
308
309	IF( forward. OR . leapf ) THEN
310
311	CALL caladvtrac(q,pbaru,pbarv,
312	* p, masse, dq, teta,
313	. flxw, pk)
314
315	IF (offline) THEN
316	Cmaf stokage du flux de masse pour traceurs OFF-LINE
317
318	#ifdef CPP_IOIPSL
319	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
320	. dtvr, itau)
321	#endif
322
323
324	ENDIF ! of IF (offline)
325	c
326	ENDIF ! of IF( forward. OR . leapf )
327
328
329	c-----------------------------------------------------------------------
330	c integrations dynamique et traceurs:
331	c ----------------------------------
332
333
334	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
335	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
336	$ finvmaold )
337
338
339	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
340	c
341	c-----------------------------------------------------------------------
342	c calcul des tendances physiques:
343	c -------------------------------
344	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
345	c
346	IF( purmats ) THEN
347	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
348	ELSE
349	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
350	ENDIF
351	c
352	c
353	IF( apphys ) THEN
354	c
355	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
356	c
357
358	CALL pression ( ip1jmp1, ap, bp, ps, p )
359	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
360
361	! rdaym_ini = itau * dtvr / daysec
362	! rdayvrai = rdaym_ini + day_ini
363	jD_cur = jD_ref + day_ini - day_ref
364	$ + int (itau * dtvr / daysec)
365	jH_cur = jH_ref + &
366	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
367	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
368	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
369	! write(lunout,*)'current date = ',an, mois, jour, secondes
370
371	c rajout debug
372	c lafin = .true.
373
374
375	c Inbterface avec les routines de phylmd (phymars ... )
376	c -----------------------------------------------------
377
378	c+jld
379
380	c Diagnostique de conservation de l'énergie : initialisation
381	IF (ip_ebil_dyn.ge.1 ) THEN
382	ztit='bil dyn'
383	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
384	IF (planet_type.eq."earth") THEN
385	CALL diagedyn(ztit,2,1,1,dtphys
386	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
387	ENDIF
388	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
389	c-jld
390	#ifdef CPP_IOIPSL
391	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
392	IF (first) THEN
393	first=.false.
394	#include "ini_paramLMDZ_dyn.h"
395	ENDIF
396	c
397	#include "write_paramLMDZ_dyn.h"
398	c
399	#endif
400	! #endif of #ifdef CPP_IOIPSL
401	CALL calfis( lafin , jD_cur, jH_cur,
402	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
403	$ du,dv,dteta,dq,
404	$ flxw,
405	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
406
407	IF (ok_strato) THEN
408	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
409	ENDIF
410
411	c ajout des tendances physiques:
412	c ------------------------------
413	CALL addfi( dtphys, leapf, forward ,
414	$ ucov, vcov, teta , q ,ps ,
415	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
416	c
417	c Diagnostique de conservation de l'énergie : difference
418	IF (ip_ebil_dyn.ge.1 ) THEN
419	ztit='bil phys'
420	IF (planet_type.eq."earth") THEN
421	CALL diagedyn(ztit,2,1,1,dtphys
422	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
423	ENDIF
424	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
425
426	ENDIF ! of IF( apphys )
427
428	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
429	! Academic case : Simple friction and Newtonan relaxation
430	! -------------------------------------------------------
431	DO l=1,llm
432	DO ij=1,ip1jmp1
433	teta(ij,l)=teta(ij,l)-dtvr*
434	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
435	ENDDO
436	ENDDO ! of DO l=1,llm
437	call friction(ucov,vcov,dtvr)
438
439	! Sponge layer (if any)
440	IF (ok_strato) THEN
441	dufi(:,:)=0.
442	dvfi(:,:)=0.
443	dtetafi(:,:)=0.
444	dqfi(:,:,:)=0.
445	dpfi(:)=0.
446	CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
447	CALL addfi( dtvr, leapf, forward ,
448	$ ucov, vcov, teta , q ,ps ,
449	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
450	ENDIF ! of IF (ok_strato)
451	ENDIF ! of IF (iflag_phys.EQ.2)
452
453
454	c-jld
455
456	CALL pression ( ip1jmp1, ap, bp, ps, p )
457	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
458
459
460	c-----------------------------------------------------------------------
461	c dissipation horizontale et verticale des petites echelles:
462	c ----------------------------------------------------------
463
464	IF(apdiss) THEN
465
466
467	c calcul de l'energie cinetique avant dissipation
468	call covcont(llm,ucov,vcov,ucont,vcont)
469	call enercin(vcov,ucov,vcont,ucont,ecin0)
470
471	c dissipation
472	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
473	ucov=ucov+dudis
474	vcov=vcov+dvdis
475	c teta=teta+dtetadis
476
477
478	c------------------------------------------------------------------------
479	if (dissip_conservative) then
480	C On rajoute la tendance due a la transform. Ec -> E therm. cree
481	C lors de la dissipation
482	call covcont(llm,ucov,vcov,ucont,vcont)
483	call enercin(vcov,ucov,vcont,ucont,ecin)
484	dtetaecdt= (ecin0-ecin)/ pk
485	c teta=teta+dtetaecdt
486	dtetadis=dtetadis+dtetaecdt
487	endif
488	teta=teta+dtetadis
489	c------------------------------------------------------------------------
490
491
492	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
493	c ... Calcul de la valeur moyenne, unique de h aux poles .....
494	c
495
496	DO l = 1, llm
497	DO ij = 1,iim
498	tppn(ij) = aire( ij ) * teta( ij ,l)
499	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
500	ENDDO
501	tpn = SSUM(iim,tppn,1)/apoln
502	tps = SSUM(iim,tpps,1)/apols
503
504	DO ij = 1, iip1
505	teta( ij ,l) = tpn
506	teta(ij+ip1jm,l) = tps
507	ENDDO
508	ENDDO
509
510	DO ij = 1,iim
511	tppn(ij) = aire( ij ) * ps ( ij )
512	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
513	ENDDO
514	tpn = SSUM(iim,tppn,1)/apoln
515	tps = SSUM(iim,tpps,1)/apols
516
517	DO ij = 1, iip1
518	ps( ij ) = tpn
519	ps(ij+ip1jm) = tps
520	ENDDO
521
522
523	END IF ! of IF(apdiss)
524
525	c ajout debug
526	c IF( lafin ) then
527	c abort_message = 'Simulation finished'
528	c call abort_gcm(modname,abort_message,0)
529	c ENDIF
530
531	c ********************************************************************
532	c ********************************************************************
533	c .... fin de l'integration dynamique et physique pour le pas itau ..
534	c ********************************************************************
535	c ********************************************************************
536
537	c preparation du pas d'integration suivant ......
538
539	IF ( .NOT.purmats ) THEN
540	c ........................................................
541	c .............. schema matsuno + leapfrog ..............
542	c ........................................................
543
544	IF(forward. OR. leapf) THEN
545	itau= itau + 1
546	c iday= day_ini+itau/day_step
547	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
548	c IF(time.GT.1.) THEN
549	c time = time-1.
550	c iday = iday+1
551	c ENDIF
552	ENDIF
553
554
555	IF( itau. EQ. itaufinp1 ) then
556	if (flag_verif) then
557	write(79,*) 'ucov',ucov
558	write(80,*) 'vcov',vcov
559	write(81,*) 'teta',teta
560	write(82,*) 'ps',ps
561	write(83,*) 'q',q
562	WRITE(85,*) 'q1 = ',q(:,:,1)
563	WRITE(86,*) 'q3 = ',q(:,:,3)
564	endif
565
566	abort_message = 'Simulation finished'
567
568	call abort_gcm(modname,abort_message,0)
569	ENDIF
570	c-----------------------------------------------------------------------
571	c ecriture du fichier histoire moyenne:
572	c -------------------------------------
573
574	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
575	IF(itau.EQ.itaufin) THEN
576	iav=1
577	ELSE
578	iav=0
579	ENDIF
580
581	IF (ok_dynzon) THEN
582	#ifdef CPP_IOIPSL
583	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
584	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
585	#endif
586	END IF
587	IF (ok_dyn_ave) THEN
588	#ifdef CPP_IOIPSL
589	CALL writedynav(itau,vcov,
590	& ucov,teta,pk,phi,q,masse,ps,phis)
591	#endif
592	ENDIF
593
594	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
595
596	c-----------------------------------------------------------------------
597	c ecriture de la bande histoire:
598	c ------------------------------
599
600	IF( MOD(itau,iecri).EQ.0) THEN
601	! Ehouarn: output only during LF or Backward Matsuno
602	if (leapf.or.(.not.leapf.and.(.not.forward))) then
603	nbetat = nbetatdem
604	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
605	unat=0.
606	do l=1,llm
607	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
608	vnat(:,l)=vcov(:,l)/cv(:)
609	enddo
610	#ifdef CPP_IOIPSL
611	if (ok_dyn_ins) then
612	! write(lunout,*) "leapfrog: call writehist, itau=",itau
613	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
614	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
615	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
616	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
617	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
618	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
619	endif ! of if (ok_dyn_ins)
620	#endif
621	! For some Grads outputs of fields
622	if (output_grads_dyn) then
623	#include "write_grads_dyn.h"
624	endif
625	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
626	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
627
628	IF(itau.EQ.itaufin) THEN
629
630
631	! if (planet_type.eq."earth") then
632	! Write an Earth-format restart file
633	CALL dynredem1("restart.nc",0.0,
634	& vcov,ucov,teta,q,masse,ps)
635	! endif ! of if (planet_type.eq."earth")
636
637	CLOSE(99)
638	ENDIF ! of IF (itau.EQ.itaufin)
639
640	c-----------------------------------------------------------------------
641	c gestion de l'integration temporelle:
642	c ------------------------------------
643
644	IF( MOD(itau,iperiod).EQ.0 ) THEN
645	GO TO 1
646	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
647
648	IF( forward ) THEN
649	c fin du pas forward et debut du pas backward
650
651	forward = .FALSE.
652	leapf = .FALSE.
653	GO TO 2
654
655	ELSE
656	c fin du pas backward et debut du premier pas leapfrog
657
658	leapf = .TRUE.
659	dt = 2.*dtvr
660	GO TO 2
661	END IF ! of IF (forward)
662	ELSE
663
664	c ...... pas leapfrog .....
665
666	leapf = .TRUE.
667	dt = 2.*dtvr
668	GO TO 2
669	END IF ! of IF (MOD(itau,iperiod).EQ.0)
670	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
671
672	ELSE ! of IF (.not.purmats)
673
674	c ........................................................
675	c .............. schema matsuno ...............
676	c ........................................................
677	IF( forward ) THEN
678
679	itau = itau + 1
680	c iday = day_ini+itau/day_step
681	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
682	c
683	c IF(time.GT.1.) THEN
684	c time = time-1.
685	c iday = iday+1
686	c ENDIF
687
688	forward = .FALSE.
689	IF( itau. EQ. itaufinp1 ) then
690	abort_message = 'Simulation finished'
691	call abort_gcm(modname,abort_message,0)
692	ENDIF
693	GO TO 2
694
695	ELSE ! of IF(forward) i.e. backward step
696
697	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
698	IF(itau.EQ.itaufin) THEN
699	iav=1
700	ELSE
701	iav=0
702	ENDIF
703
704	IF (ok_dynzon) THEN
705	#ifdef CPP_IOIPSL
706	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
707	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
708	#endif
709	ENDIF
710	IF (ok_dyn_ave) THEN
711	#ifdef CPP_IOIPSL
712	CALL writedynav(itau,vcov,
713	& ucov,teta,pk,phi,q,masse,ps,phis)
714	#endif
715	ENDIF
716
717	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
718
719	IF(MOD(itau,iecri ).EQ.0) THEN
720	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
721	nbetat = nbetatdem
722	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
723	unat=0.
724	do l=1,llm
725	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
726	vnat(:,l)=vcov(:,l)/cv(:)
727	enddo
728	#ifdef CPP_IOIPSL
729	if (ok_dyn_ins) then
730	! write(lunout,*) "leapfrog: call writehist (b)",
731	! & itau,iecri
732	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
733	endif ! of if (ok_dyn_ins)
734	#endif
735	! For some Grads outputs
736	if (output_grads_dyn) then
737	#include "write_grads_dyn.h"
738	endif
739
740	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
741
742	IF(itau.EQ.itaufin) THEN
743	! if (planet_type.eq."earth") then
744	CALL dynredem1("restart.nc",0.0,
745	& vcov,ucov,teta,q,masse,ps)
746	! endif ! of if (planet_type.eq."earth")
747	ENDIF ! of IF(itau.EQ.itaufin)
748
749	forward = .TRUE.
750	GO TO 1
751
752	ENDIF ! of IF (forward)
753
754	END IF ! of IF(.not.purmats)
755
756	STOP
757	END

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

Download in other formats:

Original Format