Context Navigation

leapfrog.F @ 1143

Last change on this file since 1143 was 1060, checked in by Laurent Fairhead, 16 years ago
Correction bug cas physique Newtonnienne (iflag_phys=2) pour ne plus passer par la physique EM/FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 19.9 KB

Line
1	!
2	c
3	c
4	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
5	& time_0)
6
7
8	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
9	USE IOIPSL
10
11	IMPLICIT NONE
12
13	c ...... Version du 10/01/98 ..........
14
15	c avec coordonnees verticales hybrides
16	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
17
18	c=======================================================================
19	c
20	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
21	c -------
22	c
23	c Objet:
24	c ------
25	c
26	c GCM LMD nouvelle grille
27	c
28	c=======================================================================
29	c
30	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
31	c et possibilite d'appeler une fonction f(y) a derivee tangente
32	c hyperbolique a la place de la fonction a derivee sinusoidale.
33
34	c ... Possibilite de choisir le shema pour l'advection de
35	c q , en modifiant iadv dans traceur.def (10/02) .
36	c
37	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
38	c Pour Van-Leer iadv=10
39	c
40	c-----------------------------------------------------------------------
41	c Declarations:
42	c -------------
43
44	#include "dimensions.h"
45	#include "paramet.h"
46	#include "comconst.h"
47	#include "comdissnew.h"
48	#include "comvert.h"
49	#include "comgeom.h"
50	#include "logic.h"
51	#include "temps.h"
52	#include "control.h"
53	#include "ener.h"
54	#include "description.h"
55	#include "serre.h"
56	#include "com_io_dyn.h"
57	#include "iniprint.h"
58	#include "advtrac.h"
59	c#include "tracstoke.h"
60
61	#include "academic.h"
62
63	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
64	! #include "clesphys.h"
65
66	integer nq
67
68	INTEGER longcles
69	PARAMETER ( longcles = 20 )
70	REAL clesphy0( longcles )
71
72	real zqmin,zqmax
73	INTEGER nbetatmoy, nbetatdem,nbetat
74
75	c variables dynamiques
76	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
77	REAL teta(ip1jmp1,llm) ! temperature potentielle
78	REAL q(ip1jmp1,llm,nqmx) ! champs advectes
79	REAL ps(ip1jmp1) ! pression au sol
80	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
81	REAL pks(ip1jmp1) ! exner au sol
82	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
83	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
84	REAL masse(ip1jmp1,llm) ! masse d'air
85	REAL phis(ip1jmp1) ! geopotentiel au sol
86	REAL phi(ip1jmp1,llm) ! geopotentiel
87	REAL w(ip1jmp1,llm) ! vitesse verticale
88
89	c variables dynamiques intermediaire pour le transport
90	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
91
92	c variables dynamiques au pas -1
93	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
94	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
95	REAL massem1(ip1jmp1,llm)
96
97	c tendances dynamiques
98	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
99	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
100
101	c tendances de la dissipation
102	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
103	REAL dtetadis(ip1jmp1,llm)
104
105	c tendances physiques
106	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
107	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
108
109	c variables pour le fichier histoire
110	REAL dtav ! intervalle de temps elementaire
111
112	REAL tppn(iim),tpps(iim),tpn,tps
113	c
114	INTEGER itau,itaufinp1,iav
115	INTEGER*4 iday ! jour julien
116	REAL time ! Heure de la journee en fraction d'1 jour
117
118	REAL SSUM
119	REAL time_0 , finvmaold(ip1jmp1,llm)
120
121	cym LOGICAL lafin
122	LOGICAL :: lafin=.false.
123	INTEGER ij,iq,l
124	INTEGER ik
125
126	real time_step, t_wrt, t_ops
127
128	REAL rdayvrai,rdaym_ini
129	LOGICAL first,callinigrads
130	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
131	save first
132	data first/.true./
133	real dt_cum, zjulian
134	character*10 infile
135	integer zan, tau0, thoriid
136	integer nid_ctesGCM
137	save nid_ctesGCM
138	real degres
139	real rlong(iip1), rlatg(jjp1)
140	real zx_tmp_2d(iip1,jjp1)
141	integer ndex2d(iip1*jjp1)
142	logical ok_sync
143	parameter (ok_sync = .true.)
144
145	data callinigrads/.true./
146	character*10 string10
147
148	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
149	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
150
151	c+jld variables test conservation energie
152	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
153	C Tendance de la temp. potentiel d (theta)/ d t due a la
154	C tansformation d'energie cinetique en energie thermique
155	C cree par la dissipation
156	REAL dtetaecdt(ip1jmp1,llm)
157	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
158	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
159	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
160	CHARACTER*15 ztit
161	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
162	!IM SAVE ip_ebil_dyn
163	!IM DATA ip_ebil_dyn/0/
164	c-jld
165
166	character*80 dynhist_file, dynhistave_file
167	character*20 modname
168	character*80 abort_message
169
170	C Calendrier
171	LOGICAL true_calendar
172	PARAMETER (true_calendar = .false.)
173
174	logical dissip_conservative
175	save dissip_conservative
176	data dissip_conservative/.true./
177
178	LOGICAL prem
179	save prem
180	DATA prem/.true./
181	INTEGER testita
182	PARAMETER (testita = 9)
183
184	logical , parameter :: flag_verif = .false.
185
186
187	integer itau_w ! pas de temps ecriture = itap + itau_phy
188
189
190	itaufin = nday*day_step
191	itaufinp1 = itaufin +1
192
193
194	itau = 0
195	iday = day_ini+itau/day_step
196	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
197	IF(time.GT.1.) THEN
198	time = time-1.
199	iday = iday+1
200	ENDIF
201
202
203	c-----------------------------------------------------------------------
204	c On initialise la pression et la fonction d'Exner :
205	c --------------------------------------------------
206
207	dq=0.
208	CALL pression ( ip1jmp1, ap, bp, ps, p )
209	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
210
211	c-----------------------------------------------------------------------
212	c Debut de l'integration temporelle:
213	c ----------------------------------
214
215	1 CONTINUE
216
217
218	#ifdef CPP_IOIPSL
219	if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
220	call guide(itau,ucov,vcov,teta,q,masse,ps)
221	else
222	IF(prt_level>9)WRITE(,)'attention on ne guide pas les ',
223	. '6 dernieres heures'
224	endif
225	#endif
226	c
227	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
228	c CALL test_period ( ucov,vcov,teta,q,p,phis )
229	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
230	c ENDIF
231	c
232	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
233	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
234	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
235	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
236	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
237
238	forward = .TRUE.
239	leapf = .FALSE.
240	dt = dtvr
241
242	c ... P.Le Van .26/04/94 ....
243
244	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
245	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
246
247	call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
248
249	2 CONTINUE
250
251	c-----------------------------------------------------------------------
252
253	c date:
254	c -----
255
256
257	c gestion des appels de la physique et des dissipations:
258	c ------------------------------------------------------
259	c
260	c ... P.Le Van ( 6/02/95 ) ....
261
262	apphys = .FALSE.
263	statcl = .FALSE.
264	conser = .FALSE.
265	apdiss = .FALSE.
266
267	IF( purmats ) THEN
268	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
269	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
270	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
271	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
272	ELSE
273	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
274	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
275	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
276	END IF
277
278	c-----------------------------------------------------------------------
279	c calcul des tendances dynamiques:
280	c --------------------------------
281
282	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
283
284	CALL caldyn
285	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
286	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
287
288	c-----------------------------------------------------------------------
289	c calcul des tendances advection des traceurs (dont l'humidite)
290	c -------------------------------------------------------------
291
292	IF( forward. OR . leapf ) THEN
293
294	CALL caladvtrac(q,pbaru,pbarv,
295	* p, masse, dq, teta,
296	. flxw, pk)
297
298	IF (offline) THEN
299	Cmaf stokage du flux de masse pour traceurs OFF-LINE
300
301	#ifdef CPP_IOIPSL
302	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
303	. dtvr, itau)
304	#endif
305
306
307	ENDIF
308	c
309	ENDIF
310
311
312	c-----------------------------------------------------------------------
313	c integrations dynamique et traceurs:
314	c ----------------------------------
315
316
317	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
318	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
319	$ finvmaold )
320
321
322	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
323	c
324	c-----------------------------------------------------------------------
325	c calcul des tendances physiques:
326	c -------------------------------
327	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
328	c
329	IF( purmats ) THEN
330	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
331	ELSE
332	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
333	ENDIF
334	c
335	c
336	IF( apphys ) THEN
337	c
338	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
339	c
340
341	CALL pression ( ip1jmp1, ap, bp, ps, p )
342	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
343
344	rdaym_ini = itau * dtvr / daysec
345	rdayvrai = rdaym_ini + day_ini
346
347
348	c rajout debug
349	c lafin = .true.
350
351
352	c Inbterface avec les routines de phylmd (phymars ... )
353	c -----------------------------------------------------
354
355	#ifdef CPP_PHYS
356	c+jld
357
358	c Diagnostique de conservation de l'énergie : initialisation
359	IF (ip_ebil_dyn.ge.1 ) THEN
360	ztit='bil dyn'
361	CALL diagedyn(ztit,2,1,1,dtphys
362	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
363	ENDIF
364	c-jld
365	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
366	IF (first) THEN
367	first=.false.
368	#include "ini_paramLMDZ_dyn.h"
369	ENDIF
370	c
371	#include "write_paramLMDZ_dyn.h"
372	c
373
374	CALL calfis( nq, lafin ,rdayvrai,time ,
375	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
376	$ du,dv,dteta,dq,
377	$ flxw,
378	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
379
380	IF (ok_strato) THEN
381	CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi)
382	ENDIF
383
384	c ajout des tendances physiques:
385	c ------------------------------
386	CALL addfi( nqmx, dtphys, leapf, forward ,
387	$ ucov, vcov, teta , q ,ps ,
388	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
389	c
390	c Diagnostique de conservation de l'énergie : difference
391	IF (ip_ebil_dyn.ge.1 ) THEN
392	ztit='bil phys'
393	CALL diagedyn(ztit,2,1,1,dtphys
394	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
395	ENDIF
396	#endif
397	ENDIF ! of IF( apphys )
398
399	IF(iflag_phys.EQ.2) THEN ! "Newtonian physics" case
400	c Calcul academique de la physique = Rappel Newtonien + friction
401	c --------------------------------------------------------------
402	teta(:,:)=teta(:,:)
403	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
404	call friction(ucov,vcov,iphysiq*dtvr)
405	ENDIF
406
407
408	c-jld
409
410	CALL pression ( ip1jmp1, ap, bp, ps, p )
411	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
412
413
414	c-----------------------------------------------------------------------
415	c dissipation horizontale et verticale des petites echelles:
416	c ----------------------------------------------------------
417
418	IF(apdiss) THEN
419
420
421	c calcul de l'energie cinetique avant dissipation
422	call covcont(llm,ucov,vcov,ucont,vcont)
423	call enercin(vcov,ucov,vcont,ucont,ecin0)
424
425	c dissipation
426	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
427	ucov=ucov+dudis
428	vcov=vcov+dvdis
429	c teta=teta+dtetadis
430
431
432	c------------------------------------------------------------------------
433	if (dissip_conservative) then
434	C On rajoute la tendance due a la transform. Ec -> E therm. cree
435	C lors de la dissipation
436	call covcont(llm,ucov,vcov,ucont,vcont)
437	call enercin(vcov,ucov,vcont,ucont,ecin)
438	dtetaecdt= (ecin0-ecin)/ pk
439	c teta=teta+dtetaecdt
440	dtetadis=dtetadis+dtetaecdt
441	endif
442	teta=teta+dtetadis
443	c------------------------------------------------------------------------
444
445
446	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
447	c ... Calcul de la valeur moyenne, unique de h aux poles .....
448	c
449
450	DO l = 1, llm
451	DO ij = 1,iim
452	tppn(ij) = aire( ij ) * teta( ij ,l)
453	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
454	ENDDO
455	tpn = SSUM(iim,tppn,1)/apoln
456	tps = SSUM(iim,tpps,1)/apols
457
458	DO ij = 1, iip1
459	teta( ij ,l) = tpn
460	teta(ij+ip1jm,l) = tps
461	ENDDO
462	ENDDO
463
464	DO ij = 1,iim
465	tppn(ij) = aire( ij ) * ps ( ij )
466	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
467	ENDDO
468	tpn = SSUM(iim,tppn,1)/apoln
469	tps = SSUM(iim,tpps,1)/apols
470
471	DO ij = 1, iip1
472	ps( ij ) = tpn
473	ps(ij+ip1jm) = tps
474	ENDDO
475
476
477	END IF
478
479	c ajout debug
480	c IF( lafin ) then
481	c abort_message = 'Simulation finished'
482	c call abort_gcm(modname,abort_message,0)
483	c ENDIF
484
485	c ********************************************************************
486	c ********************************************************************
487	c .... fin de l'integration dynamique et physique pour le pas itau ..
488	c ********************************************************************
489	c ********************************************************************
490
491	c preparation du pas d'integration suivant ......
492
493	IF ( .NOT.purmats ) THEN
494	c ........................................................
495	c .............. schema matsuno + leapfrog ..............
496	c ........................................................
497
498	IF(forward. OR. leapf) THEN
499	itau= itau + 1
500	iday= day_ini+itau/day_step
501	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
502	IF(time.GT.1.) THEN
503	time = time-1.
504	iday = iday+1
505	ENDIF
506	ENDIF
507
508
509	IF( itau. EQ. itaufinp1 ) then
510	if (flag_verif) then
511	write(79,*) 'ucov',ucov
512	write(80,*) 'vcov',vcov
513	write(81,*) 'teta',teta
514	write(82,*) 'ps',ps
515	write(83,*) 'q',q
516	WRITE(85,*) 'q1 = ',q(:,:,1)
517	WRITE(86,*) 'q3 = ',q(:,:,3)
518	endif
519
520	abort_message = 'Simulation finished'
521
522	call abort_gcm(modname,abort_message,0)
523	ENDIF
524	c-----------------------------------------------------------------------
525	c ecriture du fichier histoire moyenne:
526	c -------------------------------------
527
528	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
529	IF(itau.EQ.itaufin) THEN
530	iav=1
531	ELSE
532	iav=0
533	ENDIF
534	#ifdef CPP_IOIPSL
535	CALL writedynav(histaveid, nqmx, itau,vcov ,
536	, ucov,teta,pk,phi,q,masse,ps,phis)
537	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
538	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
539	#endif
540
541	ENDIF
542
543	c-----------------------------------------------------------------------
544	c ecriture de la bande histoire:
545	c ------------------------------
546
547	IF( MOD(itau,iecri ).EQ.0) THEN
548	c IF( MOD(itau,iecri*day_step).EQ.0) THEN
549
550	nbetat = nbetatdem
551	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
552	unat=0.
553	do l=1,llm
554	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
555	vnat(:,l)=vcov(:,l)/cv(:)
556	enddo
557	#ifdef CPP_IOIPSL
558	c CALL writehist(histid,histvid, nqmx,itau,vcov,
559	c s ucov,teta,phi,q,masse,ps,phis)
560	#else
561	#include "write_grads_dyn.h"
562	#endif
563
564
565	ENDIF
566
567	IF(itau.EQ.itaufin) THEN
568
569
570	#ifdef CPP_IOIPSL
571	CALL dynredem1("restart.nc",0.0,
572	, vcov,ucov,teta,q,nqmx,masse,ps)
573	#endif
574
575	CLOSE(99)
576	ENDIF
577
578	c-----------------------------------------------------------------------
579	c gestion de l'integration temporelle:
580	c ------------------------------------
581
582	IF( MOD(itau,iperiod).EQ.0 ) THEN
583	GO TO 1
584	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
585
586	IF( forward ) THEN
587	c fin du pas forward et debut du pas backward
588
589	forward = .FALSE.
590	leapf = .FALSE.
591	GO TO 2
592
593	ELSE
594	c fin du pas backward et debut du premier pas leapfrog
595
596	leapf = .TRUE.
597	dt = 2.*dtvr
598	GO TO 2
599	END IF
600	ELSE
601
602	c ...... pas leapfrog .....
603
604	leapf = .TRUE.
605	dt = 2.*dtvr
606	GO TO 2
607	END IF
608
609	ELSE
610
611	c ........................................................
612	c .............. schema matsuno ...............
613	c ........................................................
614	IF( forward ) THEN
615
616	itau = itau + 1
617	iday = day_ini+itau/day_step
618	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
619
620	IF(time.GT.1.) THEN
621	time = time-1.
622	iday = iday+1
623	ENDIF
624
625	forward = .FALSE.
626	IF( itau. EQ. itaufinp1 ) then
627	abort_message = 'Simulation finished'
628	call abort_gcm(modname,abort_message,0)
629	ENDIF
630	GO TO 2
631
632	ELSE
633
634	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
635	IF(itau.EQ.itaufin) THEN
636	iav=1
637	ELSE
638	iav=0
639	ENDIF
640	#ifdef CPP_IOIPSL
641	CALL writedynav(histaveid, nqmx, itau,vcov ,
642	, ucov,teta,pk,phi,q,masse,ps,phis)
643	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
644	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
645	#endif
646
647	ENDIF
648
649	IF(MOD(itau,iecri ).EQ.0) THEN
650	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
651	nbetat = nbetatdem
652	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
653	unat=0.
654	do l=1,llm
655	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
656	vnat(:,l)=vcov(:,l)/cv(:)
657	enddo
658	#ifdef CPP_IOIPSL
659	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
660	c , ucov,teta,phi,q,masse,ps,phis)
661	#else
662	#include "write_grads_dyn.h"
663	#endif
664
665
666	ENDIF
667
668	#ifdef CPP_IOIPSL
669	IF(itau.EQ.itaufin)
670	. CALL dynredem1("restart.nc",0.0,
671	. vcov,ucov,teta,q,nqmx,masse,ps)
672	#endif
673
674	forward = .TRUE.
675	GO TO 1
676
677	ENDIF
678
679	END IF
680
681	STOP
682	END

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

Download in other formats:

Original Format