source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn3d_common/prather.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 17.5 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt)
5      IMPLICIT NONE
6
7c=======================================================================
8c   Adaptation LMDZ:  A.Armengaud (LGGE)
9c   ----------------
10c
11c   ************************************************
12c   Transport des traceurs par la methode de prather
13c   Ref : 
14c
15c   ************************************************
16c   q,w,pext,pbaru et pbarv : arguments d'entree  pour le s-pg
17c
18c=======================================================================
19
20
21!-----------------------------------------------------------------------
22!   INCLUDE 'dimensions.h'
23!
24!   dimensions.h contient les dimensions du modele
25!   ndm est tel que iim=2**ndm
26!-----------------------------------------------------------------------
27
28      INTEGER iim,jjm,llm,ndm
29
30      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
31
32!-----------------------------------------------------------------------
33!
34! $Header$
35!
36!
37!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
38!                 veillez  n'utiliser que des ! pour les commentaires
39!                 et  bien positionner les & des lignes de continuation
40!                 (les placer en colonne 6 et en colonne 73)
41!
42!
43!-----------------------------------------------------------------------
44!   INCLUDE 'paramet.h'
45
46      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
47      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
48      INTEGER  ijmllm,mvar
49      INTEGER jcfil,jcfllm
50
51      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
52     &    ,jjp1=jjm+1-1/jjm)
53      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
54      PARAMETER( kftd  = iim/2 -ndm )
55      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
56      PARAMETER( ip1jmi1= ip1jm - iip1 )
57      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
58      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
59      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
60
61!-----------------------------------------------------------------------
62!
63! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
64!
65!-----------------------------------------------------------------------
66! INCLUDE comconst.h
67
68      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
69     &                 iflag_top_bound,mode_top_bound
70      COMMON/comconstr/dtvr,daysec,                                     &
71     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
72     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
73     & ,dissip_pupstart  ,tau_top_bound,                                &
74     & daylen,molmass, ihf
75      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
76
77      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
78      REAL dtvr ! dynamical time step (in s)
79      REAL daysec !length (in s) of a standard day
80      REAL pi    ! something like 3.14159....
81      REAL dtphys ! (s) time step for the physics
82      REAL dtdiss ! (s) time step for the dissipation
83      REAL rad ! (m) radius of the planet
84      REAL r ! Reduced Gas constant r=R/mu
85             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
86      REAL cpp   ! Cp
87      REAL kappa ! kappa=R/Cp
88      REAL cotot
89      REAL unsim ! = 1./iim
90      REAL g ! (m/s2) gravity
91      REAL omeg ! (rad/s) rotation rate of the planet
92! Dissipation factors, for Earth model:
93      REAL dissip_factz,dissip_zref !dissip_deltaz
94! Dissipation factors, for other planets:
95      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
96      REAL dissip_pupstart
97      INTEGER iflag_top_bound,mode_top_bound
98      REAL tau_top_bound
99      REAL daylen ! length of solar day, in 'standard' day length
100      REAL molmass ! (g/mol) molar mass of the atmosphere
101
102      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
103      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
104
105
106!-----------------------------------------------------------------------
107!
108! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
109!
110!-----------------------------------------------------------------------
111!   INCLUDE 'comvert.h'
112
113      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
114     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
115     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
116
117      common/comverti/disvert_type, pressure_exner
118
119      real ap     ! hybrid pressure contribution at interlayers
120      real bp     ! hybrid sigma contribution at interlayer
121      real presnivs ! (reference) pressure at mid-layers
122      real dpres
123      real pa     ! reference pressure (Pa) at which hybrid coordinates
124                  ! become purely pressure
125      real preff  ! reference surface pressure (Pa)
126      real nivsigs
127      real nivsig
128      real aps    ! hybrid pressure contribution at mid-layers
129      real bps    ! hybrid sigma contribution at mid-layers
130      real scaleheight ! atmospheric (reference) scale height (km)
131      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
132                     ! preff and scaleheight
133
134      integer disvert_type ! type of vertical discretization:
135                           ! 1: Earth (default for planet_type==earth),
136                           !     automatic generation
137                           ! 2: Planets (default for planet_type!=earth),
138                           !     using 'z2sig.def' (or 'esasig.def) file
139
140      logical pressure_exner
141!     compute pressure inside layers using Exner function, else use mean
142!     of pressure values at interfaces
143
144 !-----------------------------------------------------------------------
145!
146! $Header$
147!
148!CDK comgeom2
149      COMMON/comgeom/                                                   &
150     & cu(iip1,jjp1),cv(iip1,jjm),unscu2(iip1,jjp1),unscv2(iip1,jjm)  , &
151     & aire(iip1,jjp1),airesurg(iip1,jjp1),aireu(iip1,jjp1)           , &
152     & airev(iip1,jjm),unsaire(iip1,jjp1),apoln,apols                 , &
153     & unsairez(iip1,jjm),airuscv2(iip1,jjm),airvscu2(iip1,jjm)       , &
154     & aireij1(iip1,jjp1),aireij2(iip1,jjp1),aireij3(iip1,jjp1)       , &
155     & aireij4(iip1,jjp1),alpha1(iip1,jjp1),alpha2(iip1,jjp1)         , &
156     & alpha3(iip1,jjp1),alpha4(iip1,jjp1),alpha1p2(iip1,jjp1)        , &
157     & alpha1p4(iip1,jjp1),alpha2p3(iip1,jjp1),alpha3p4(iip1,jjp1)    , &
158     & fext(iip1,jjm),constang(iip1,jjp1), rlatu(jjp1),rlatv(jjm),      &
159     & rlonu(iip1),rlonv(iip1),cuvsurcv(iip1,jjm),cvsurcuv(iip1,jjm)  , &
160     & cvusurcu(iip1,jjp1),cusurcvu(iip1,jjp1)                        , &
161     & cuvscvgam1(iip1,jjm),cuvscvgam2(iip1,jjm),cvuscugam1(iip1,jjp1), &
162     & cvuscugam2(iip1,jjp1),cvscuvgam(iip1,jjm),cuscvugam(iip1,jjp1) , &
163     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2                , &
164     & unsair_gam1(iip1,jjp1),unsair_gam2(iip1,jjp1)                  , &
165     & unsairz_gam(iip1,jjm),aivscu2gam(iip1,jjm),aiuscv2gam(iip1,jjm)  &
166     & , xprimu(iip1),xprimv(iip1)
167
168
169      REAL                                                               &
170     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,apoln,apols,unsaire &
171     & ,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4     , &
172     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 , &
173     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     , &
174     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1           , &
175     & unsapolnga2,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2     , &
176     & unsairz_gam,aivscu2gam,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu    , &
177     & cusurcvu,xprimu,xprimv
178
179c   Arguments:
180c   ----------
181      INTEGER iq,nt
182      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
183      REAL masse(iip1,jjp1,llm)
184      REAL q( iip1,jjp1,llm,0:9)
185      REAL w( ip1jmp1,llm )
186      integer ordre,ilim
187
188c   Local:
189c   ------
190      LOGICAL limit
191      real zq(iip1,jjp1,llm)
192      REAL sm ( iip1,jjp1, llm )
193      REAL s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
194      REAL sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )
195      REAL sxx( iip1,jjp1,llm)
196      REAL sxy( iip1,jjp1,llm)
197      REAL sxz( iip1,jjp1,llm)
198      REAL syy( iip1,jjp1,llm )
199      REAL syz( iip1,jjp1,llm )
200      REAL szz( iip1,jjp1,llm ),zz
201      INTEGER i,j,l,indice
202      real sxn(iip1),sxs(iip1)
203
204      real sinlon(iip1),sinlondlon(iip1)
205      real coslon(iip1),coslondlon(iip1)
206      real qmin,qmax
207      save qmin,qmax
208      save sinlon,coslon,sinlondlon,coslondlon
209      real dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
210      real masn,mass
211c
212      REAL      SSUM
213      integer ismax,ismin
214      EXTERNAL  SSUM, convflu,ismin,ismax
215      logical first
216      save first
217      EXTERNAL advxp,advyp,advzp 
218
219
220      data first/.true./
221      data qmin,qmax/-1.e33,1.e33/
222
223
224c==========================================================================
225c==========================================================================
226c     MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
227c==========================================================================
228c==========================================================================
229      REAL dt
230c==========================================================================
231      limit = .TRUE.
232 
233      if(first) then
234         print*,'SCHEMA PRATHER'
235         first=.false.
236         do i=2,iip1
237            coslon(i)=cos(rlonv(i))
238            sinlon(i)=sin(rlonv(i))
239            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
240            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
241         enddo
242         coslon(1)=coslon(iip1)
243         coslondlon(1)=coslondlon(iip1)
244         sinlon(1)=sinlon(iip1)
245         sinlondlon(1)=sinlondlon(iip1)
246
247        DO l = 1,llm
248        DO j = 1,jjp1
249        DO i = 1,iip1
250        q( i,j,l,1 )=0.
251        q( i,j,l,2)=0.
252        q( i,j,l,3)=0.
253        q( i,j,l,4)=0.
254        q( i,j,l,5)=0.
255        q( i,j,l,6)=0.
256        q( i,j,l,7)=0.
257        q( i,j,l,8)=0.
258        q( i,j,l,9)=0.
259        ENDDO
260        ENDDO
261        ENDDO
262      endif
263c   Fin modif Fred
264
265c *** On calcule la masse d'air en kg
266
267       DO l = 1,llm
268        DO j = 1,jjp1
269         DO i = 1,iip1
270         sm( i,j,llm+1-l ) =masse(i,j,l)
271         ENDDO
272        ENDDO
273       ENDDO
274
275c *** q contient les qqtes de traceur avant l'advection 
276
277c *** Affectation des tableaux S a partir de Q
278 
279       DO l = 1,llm
280        DO j = 1,jjp1
281         DO i = 1,iip1
282       s0( i,j,l) = q ( i,j,llm+1-l,0 )*sm(i,j,l)
283       sx( i,j,l) = q( i,j,llm+1-l,1 )*sm(i,j,l)
284       sy( i,j,l) = q( i,j,llm+1-l,2)*sm(i,j,l)
285       sz( i,j,l) = q( i,j,llm+1-l,3)*sm(i,j,l)
286       sxx( i,j,l) = q( i,j,llm+1-l,4)*sm(i,j,l)
287       sxy( i,j,l) = q( i,j,llm+1-l,5)*sm(i,j,l)
288       sxz( i,j,l) = q( i,j,llm+1-l,6)*sm(i,j,l)
289       syy( i,j,l) = q( i,j,llm+1-l,7)*sm(i,j,l)
290       syz( i,j,l) = q( i,j,llm+1-l,8)*sm(i,j,l)
291       szz( i,j,l) = q( i,j,llm+1-l,9)*sm(i,j,l)
292         ENDDO
293        ENDDO
294       ENDDO
295c *** Appel des subroutines d'advection en X, en Y et en Z
296c *** Advection avec "time-splitting"
297     
298c-----------------------------------------------------------
299       do indice =1,nt
300       call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
301     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
302        end do
303        do l=1,llm
304        do i=1,iip1
305        sy(i,1,l)=0.
306        sy(i,jjp1,l)=0.
307        enddo
308        enddo
309c---------------------------------------------------------
310       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
311     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
312c---------------------------------------------------------
313
314c---------------------------------------------------------
315       do j=1,jjp1
316          do i=1,iip1
317             sz(i,j,1)=0.
318             sz(i,j,llm)=0.
319             sxz(i,j,1)=0.
320             sxz(i,j,llm)=0.
321             syz(i,j,1)=0.
322             syz(i,j,llm)=0.
323             szz(i,j,1)=0.
324             szz(i,j,llm)=0.
325          enddo
326       enddo
327       call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz
328     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
329        do l=1,llm
330        do i=1,iip1
331        sy(i,1,l)=0.
332        sy(i,jjp1,l)=0.
333        enddo
334        enddo
335
336c---------------------------------------------------------
337
338c---------------------------------------------------------
339       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
340     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
341c---------------------------------------------------------
342       DO l = 1,llm
343        DO j = 1,jjp1
344             s0( iip1,j,l)=s0( 1,j,l )
345             sx( iip1,j,l)=sx( 1,j,l )
346             sy( iip1,j,l)=sy( 1,j,l )
347             sz( iip1,j,l)=sz( 1,j,l )
348             sxx( iip1,j,l)=sxx( 1,j,l )
349             sxy( iip1,j,l)=sxy( 1,j,l)
350             sxz( iip1,j,l)=sxz( 1,j,l )
351             syy( iip1,j,l)=syy( 1,j,l )
352             syz( iip1,j,l)=syz( 1,j,l)
353             szz( iip1,j,l)=szz( 1,j,l )
354        ENDDO
355       ENDDO
356       do indice=1,nt
357       call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
358     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
359        end do
360c---------------------------------------------------------
361c---------------------------------------------------------
362c ***   On repasse les S dans la variable qpr
363c ***   On repasse les S dans la variable q directement 14/10/94
364
365       DO  l = 1,llm
366        DO  j = 1,jjp1
367         DO  i = 1,iip1
368      q( i,j,llm+1-l,0 )=s0( i,j,l )/sm(i,j,l)
369      q( i,j,llm+1-l,1 ) = sx( i,j,l )/sm(i,j,l)
370      q( i,j,llm+1-l,2 ) = sy( i,j,l )/sm(i,j,l)
371      q( i,j,llm+1-l,3 ) = sz( i,j,l )/sm(i,j,l)
372      q( i,j,llm+1-l,4 ) = sxx( i,j,l )/sm(i,j,l)
373      q( i,j,llm+1-l,5 ) = sxy( i,j,l )/sm(i,j,l)
374      q( i,j,llm+1-l,6 ) = sxz( i,j,l )/sm(i,j,l)
375      q( i,j,llm+1-l,7 ) = syy( i,j,l )/sm(i,j,l)
376      q( i,j,llm+1-l,8 ) = syz( i,j,l )/sm(i,j,l)
377      q( i,j,llm+1-l,9 ) = szz( i,j,l )/sm(i,j,l)
378      ENDDO
379      ENDDO
380      ENDDO
381
382c---------------------------------------------------------
383c      go to  777
384c   filtrages aux poles
385
386c Traitements specifiques au pole
387
388c   filtrages aux poles
389         DO l=1,llm
390c   filtrages aux poles
391         masn=ssum(iim,sm(1,1,l),1)
392         mass=ssum(iim,sm(1,jjp1,l),1)
393         qpn=ssum(iim,s0(1,1,l),1)/masn
394         qps=ssum(iim,s0(1,jjp1,l),1)/mass
395         dqzpn=ssum(iim,sz(1,1,l),1)/masn
396         dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
397         do i=1,iip1
398          q( i,1,llm+1-l,3)=dqzpn
399          q( i,jjp1,llm+1-l,3)=dqzps
400          q( i,1,llm+1-l,0)=qpn
401          q( i,jjp1,llm+1-l,0)=qps
402         enddo
403c       enddo
404c         print*,'qpn',qpn,'qps',qps
405c          print*,'dqzpn',dqzpn,'dqzps',dqzps
406c       enddo
407           dyn1=0.
408           dys1=0.
409           dyn2=0.
410           dys2=0.
411        do i=1,iim
412        zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
413        dyn1=dyn1+sinlondlon(i)*zz
414        dyn2=dyn2+coslondlon(i)*zz
415        zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
416        dys1=dys1+sinlondlon(i)*zz
417        dys2=dys2+coslondlon(i)*zz
418        enddo
419         do i=1,iim
420         q(i,1,llm+1-l,2)=
421     $   (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
422         q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)
423     $          +q(i,1,llm+1-l,2)
424         q(i,jjp1,llm+1-l,2)=
425     $   (sinlon(i)*dys1+coslon(i)*dys2)/2.
426         q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
427     $      -q(i,jjp1,llm+1-l,2)
428         enddo
429      q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
430      q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
431      do i=1,iim
432      sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
433      sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
434      enddo
435      sxn(iip1)=sxn(1)
436      sxs(iip1)=sxs(1)
437      do i=1,iim
438      q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
439      q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
440      END DO
441      q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
442      q(1,jjp1,llm+1-l,1)=
443     $   q(iip1,jjp1,llm+1-l,1)
444        enddo
445         do l=1,llm
446           do i=1,iim
447            q( i,1,llm+1-l,4)=0.
448            q( i,jjp1,llm+1-l,4)=0.
449            q( i,1,llm+1-l,5)=0.
450            q( i,jjp1,llm+1-l,5)=0.
451            q( i,1,llm+1-l,6)=0.
452            q( i,jjp1,llm+1-l,6)=0.
453            q( i,1,llm+1-l,7)=0.
454            q( i,jjp1,llm+1-l,7)=0.
455            q( i,1,llm+1-l,8)=0.
456            q( i,jjp1,llm+1-l,8)=0.
457            q( i,1,llm+1-l,9)=0.
458            q( i,jjp1,llm+1-l,9)=0.
459          enddo
460         ENDDO
461
462777      continue
463c
464c   bouclage en longitude
465      do l=1,llm
466      do j=1,jjp1
467      q(iip1,j,l,0)=q(1,j,l,0)
468      q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)
469      q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)
470      q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)
471      q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)
472      q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)
473      q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)
474      q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)
475      q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)
476      q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)
477      q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)
478      enddo
479      enddo
480        DO l = 1,llm
481         DO j = 2,jjm
482           DO i = 1,iip1
483         IF (q(i,j,l,0).lt.0.)  THEN
484         PRINT*,'------------ BIP-----------'
485         PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),
486     $          q(i,j-1,l,0)
487         PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
488         PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),
489     $   q(i,j-1,l,2)   
490         PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
491c                    PRINT*,' PBL EN SORTIE D'' ADVZP'
492                     q(i,j,l,0)=0.
493c                  STOP
494               ENDIF
495           ENDDO
496         ENDDO
497         do j=1,jjp1,jjm
498         do i=1,iip1
499               IF (q(i,j,l,0).lt.0.)  THEN
500               PRINT*,'------------ BIP 2-----------'
501         PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)
502         PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
503         PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)
504         PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
505
506                     q(i,j,l,0)=0.
507c                  STOP
508               ENDIF
509         enddo
510         enddo
511        ENDDO
512      RETURN
513      END
Note: See TracBrowser for help on using the repository browser.