SUBROUTINE pltall(ndim),16
!*********************************************************************
! plot ndim-variables in x-y planes for all grids
!*********************************************************************
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'grddsc.inc'
!*********************************************************************
! plot x-y planes of ndim-D variables on grid mptr
!*********************************************************************
CALL xcltyp(0)
CALL xcfull
CALL xhilit(0)
mptr=1
!*********************************************************************
! goto 44
! plot vector component ivar=2 u
!*********************************************************************
ivar = 2
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltvxy
(ndim,2,5,11, 1, 2)
CALL xframe
CALL pltvxy(ndim,2,5,11, 1, nzz/2)
CALL xframe
!*********************************************************************
! plot vector component ivar=5 v
!*********************************************************************
ivar = 5
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltvxy(ndim,2,5,11, 2, 2)
CALL xframe
CALL pltvxy(ndim,2,5,11, 2, nzz/2)
CALL xframe
!*********************************************************************
! plot scalar w
!*********************************************************************
ivar=8
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltsxy(ndim,ivar,2)
CALL xframe
CALL pltsxy(ndim,ivar,nzz/2)
CALL xframe
!*********************************************************************
! plot scalar ptprt
!*********************************************************************
ivar=11
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltsxy(ndim,ivar,2)
CALL xframe
CALL pltsxy(ndim,ivar,nzz/2)
CALL xframe
!*********************************************************************
! plot scalar pprt
!*********************************************************************
ivar=14
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltsxy(ndim,ivar,2)
CALL xframe
CALL pltsxy(ndim,ivar,nzz/2)
CALL xframe
!*********************************************************************
! plot scalar qvprt
!*********************************************************************
ivar=17
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltsxy(ndim,ivar,2)
CALL xframe
CALL pltsxy(ndim,ivar,nzz/2)
CALL xframe
!*********************************************************************
! plot scalar qcprt
!*********************************************************************
ivar=20
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltsxy(ndim,ivar,2)
CALL xframe
CALL pltsxy(ndim,ivar,nzz/2)
CALL xframe
!*********************************************************************
! plot scalar qrprt
!*********************************************************************
ivar=23
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL pltsxy(ndim,ivar,2)
CALL xframe
CALL pltsxy(ndim,ivar,nzz/2)
CALL xframe
!*********************************************************************
RETURN
END SUBROUTINE pltall
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE gtrcrd(mptr,ndim,ivar,x,y,nx1,ny1) 1,1
!*********************************************************************
! To return the coordinate arrays x and y of variable ivar on grid mptr
! relative to its own grid origin
! nx1, ny1 give the dimension of the portion of array actually used.
!*********************************************************************
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'grddsc.inc'
REAL :: x(1), y(1)
CALL chkdim(ndim,'GTRCRD')
nx=node(5,mptr)
ny=node(6,mptr)
cosr=rnode(21,mptr)
sinr=rnode(22,mptr)
dx =rnode(9 ,mptr)
dy =rnode(10,mptr)
xor=rnode(1,mptr)
yor=rnode(2,mptr)
IF(ndim == 2) THEN
IF(ivar > nxy2d)THEN
WRITE(6,5) ivar,nxy2d
WRITE(6,'(1x,a)') 'Job stopped in GTRCRD!'
STOP
END IF
xshift=stgxy(1,ivar)
yshift=stgxy(2,ivar)
nx1=nx-idmxy(1,ivar)
ny1=ny-idmxy(2,ivar)
ELSE
xshift=stgxyz(1,ivar)
yshift=stgxyz(2,ivar)
nx1=nx-idmxyz(1,ivar)
ny1=ny-idmxyz(2,ivar)
IF(ivar > nxyz3d)THEN
WRITE(6,10) ivar,nxyz3d
WRITE(6,'(1x,a)') 'Job stopped in GTRCRD!'
STOP
END IF
END IF
DO j=1,ny
DO i=1,nx
ij=i+(j-1)*nx
x(ij)= (i-1+xshift)*dx
y(ij)= (j-1+yshift)*dy
END DO
END DO
5 FORMAT &
('Error in GTACRD: variable number for 2-d x-y array ivar=' &
,i3,' > nxy2d=',i3)
10 FORMAT('Error in GTACRD: variable number for 3-d array ivar=',i3, &
' > nxyz3d=',i3)
RETURN
END SUBROUTINE gtrcrd
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE gtacrd(mptr,ndim,ivar,x,y,nx1,ny1) 2,1
!*********************************************************************
! To return the absolute coordinate arrays x and y of variable ivar
! on grid mptr relative to the absolute origin.
! nx1, ny1 give the dimension of the portion of array actually used.
!*********************************************************************
INCLUDE 'nodal.inc'
INCLUDE 'grddsc.inc'
REAL :: x(1),y(1)
CALL chkdim(ndim,'GTRCRD')
nx=node(5,mptr)
ny=node(6,mptr)
cosr=rnode(21,mptr)
sinr=rnode(22,mptr)
dx =rnode(9 ,mptr)
dy =rnode(10,mptr)
xor=rnode(1,mptr)
yor=rnode(2,mptr)
IF(ndim == 2) THEN
IF(ivar > nxy2d)THEN
WRITE(6,5) ivar,nxy2d
WRITE(6,'(1x,a)') 'Job stopped in GTACRD!'
STOP
END IF
xshift=stgxy(1,ivar)
yshift=stgxy(2,ivar)
nx1=nx-idmxy(1,ivar)
ny1=ny-idmxy(2,ivar)
ELSE
IF(ivar > nxyz3d)THEN
WRITE(6,10) ivar,nxyz3d
WRITE(6,'(1x,a)') 'Job stopped in GTACRD!'
STOP
END IF
xshift=stgxyz(1,ivar)
yshift=stgxyz(2,ivar)
nx1=nx-idmxyz(1,ivar)
ny1=ny-idmxyz(2,ivar)
END IF
DO j=1,ny
DO i=1,nx
xx = (i-1+xshift)*dx
yy = (j-1+yshift)*dy
ij=i+(j-1)*nx
x(ij) = xor+xx*cosr-yy*sinr
y(ij) = yor+xx*sinr+yy*cosr
END DO
END DO
5 FORMAT &
('Error in GTACRD: variable number for 2-d x-y array ivar=', &
i3,' > nxy2d=',i3)
10 FORMAT('Error in GTACRD: variable number for 3-d array ivar=', &
i3,' > nxyz3d=',i3)
RETURN
END SUBROUTINE gtacrd
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE gtxypl( a,tmp,m,n,l,k ) 2
!*********************************************************************
! copy a x-y plane from a 3-d field
!*********************************************************************
DIMENSION a(m,n,l),tmp(m,n)
DO j=1,n
DO i=1,m
tmp(i,j) = a(i,j,k)
END DO
END DO
RETURN
END SUBROUTINE gtxypl
!*********************************************************************
!*********************************************************************
!*********************************************************************
INTEGER FUNCTION nxtgrd(mptr0)
!*********************************************************************
! return a grid pointer that is next to mptr0 on the same level or the
! first one on the level that is one-level finer than the level of mptr0
!*********************************************************************
INCLUDE 'nodal.inc'
mptr=mptr0
IF(mptr == 0) RETURN
level=node(4,mptr)
5 mptr=node(10,mptr)
IF(mptr /= 0) GO TO 20
level = level+1
IF(level > lfine)THEN
nxtgrd=0
GO TO 999
END IF
nxtgrd=lstart(level)
GO TO 999
20 nxtgrd=mptr
999 RETURN
END FUNCTION nxtgrd
!*********************************************************************
!*********************************************************************
!*********************************************************************
INTEGER FUNCTION lstgrd(mptr0)
!*********************************************************************
! return a grid pointer that is next to mptr0 on the same level or the
! first one on the level that is one-level coarser than the level of mptr0
!*********************************************************************
INCLUDE 'nodal.inc'
mptr=mptr0
IF(mptr == 0) RETURN
level=node(4,mptr)
5 mptr=node(10,mptr)
IF(mptr /= 0) GO TO 20
level = level-1
IF(level <= 0)THEN
lstgrd=0
GO TO 999
END IF
lstgrd=lstart(level)
GO TO 999
20 lstgrd=mptr
999 RETURN
END FUNCTION lstgrd
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE pltsxy(ndim,ivar, kk) 12,5
!*********************************************************************
! plot x-y slices of ndim-D scalar ivar at k=kk.
!*********************************************************************
INCLUDE 'nodal.inc'
INCLUDE 'agrialloc.inc'
CHARACTER (LEN=132) :: title
pi=4.0*ATAN(1.0)
CALL xqmap (x1,x2,y1,y2)
CALL xqlmsk(nmask)
!*********************************************************************
! plot a over all frame and axes.
!*********************************************************************
mptr=mstart
nx=node(5,mptr)
ny=node(6,mptr)
nz=node(14,mptr)
pi=4.0*ATAN(1.0)
xor=rnode(1, mptr)
yor=rnode(2, mptr)
dx=rnode(9 , mptr)
dy=rnode(10, mptr)
xl=(nx-1)*dx
yl=(ny-1)*dy
IF(yl > xl) THEN
xm = (1.0-0.8*xl/yl)/2
CALL xpspac( xm, 1.0-xm, 0.1, 0.9)
ELSE
ym = (1.0-0.8*yl/xl)/2
CALL xpspac(0.1,0.9, ym, 1.0-ym)
END IF
CALL xmap (0.0, xl, 0.0, yl)
CALL xqpspc(p1,p2,p3,p4)
CALL xqmap(xa,xb,ya,yb)
CALL xaxsca(0.0, xl, dx, 0.0, yl, dy)
WRITE(title,'(i1,''-D VAR NO.'',I2,'' K='',I3)')ndim,ivar,kk
CALL xchsiz(0.025*yl)
CALL xcharc(0.5*xl, 1.05*yl, title(1:21))
PRINT*
PRINT*,'Plotting ',title(1:21), ' for grid ', mptr
!*********************************************************************
! loop over all grids
!*********************************************************************
mptr = lstart(lfine)
100 CONTINUE
nx=node(5,mptr)
ny=node(6,mptr)
nz=node(14,mptr)
CALL xqlmsk(level)
CALL xunmsk(1)
IF(mptr /= mstart) THEN
CALL xpenup(rnode(7,mptr),rnode(8,mptr))
DO i=1,7,2
CALL xpendn(rnode(i,mptr),rnode(i+1,mptr))
END DO
END IF
CALL xrsmsk(level)
IF(ndim == 3) in1 = igtxyz(mptr,ivar,1 )
IF(ndim == 2) in1 = igtnxy(mptr,ivar,1 )
n3rd=nz
IF(ndim == 2) n3rd=1
!
needsp = nx*ny
ntem1=igetsp(needsp)
ntem2=igetsp(needsp)
ntem3=igetsp(needsp)
CALL gtxypl( CALL gtacrd(mptr,ndim,ivar,a(ntem2),a(ntem3),nxa,nya)
IF(mptr == lstart(lfine) ) THEN
zinc=1.0
mode=1
ELSE
mode=2 ! now use the zinc of first grid lstart(lfine).
END IF
IF(mptr == mstart)THEN
CALL xconts2( zinc,mode)
ELSE
CALL xconts2( zinc,mode)
xor=rnode(1, mptr)
yor=rnode(2, mptr)
dx =rnode(9 , mptr)
dy =rnode(10, mptr)
xl=(nx-1)*dx
yl=(ny-1)*dy
cosr=rnode(21,mptr)
sinr=rnode(22,mptr)
ang = ASIN (sinr )*180/pi
xo = xor+0.5*(dx*cosr-dy*sinr)
yo = yor+0.5*(dx*sinr+dy*cosr)
CALL xrmask(xo,yo,xl-dx,yl-dx,ang)
END IF
CALL resett
mptr= lstgrd(mptr)
IF(mptr /= 0)GO TO 100
CALL xmap(x1,x2,y1,y2)
CALL xrsmsk(nmask)
RETURN
END SUBROUTINE pltsxy
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE pltvxy(ndim,ivar1,ivar2,iscalr, kk, iplot) 4,5
!*********************************************************************
! plot contour map of one of two components (ivar1, ivar2) of
! a ndim-D vector projected to absolute coordinate system, i.e.
! that of the base grid. The coordinates for scalar variable iscalr
! are used.
!
! iplot =1, ivar1 isplotted, else ivar2 is plotted.
!*********************************************************************
INCLUDE 'nodal.inc'
INCLUDE 'agrialloc.inc'
CHARACTER (LEN=132) :: title
pi=4.0*ATAN(1.0)
CALL xqmap (x1,x2,y1,y2)
CALL xqlmsk(nmask)
IF(iplot == 1) THEN
ivar=ivar1
ELSE
ivar=ivar2
END IF
mptr=mstart
nx=node(5,mptr)
ny=node(6,mptr)
nz=node(14,mptr)
pi=4.0*ATAN(1.0)
xor=rnode(1, mptr)
yor=rnode(2, mptr)
dx=rnode(9 , mptr)
dy=rnode(10, mptr)
xl=(nx-1)*dx
yl=(ny-1)*dy
IF(yl > xl) THEN
xm = (1.0-0.8*xl/yl)/2
CALL xpspac( xm, 1.0-xm, 0.1, 0.9)
ELSE
ym = (1.0-0.8*yl/xl)/2
CALL xpspac(0.1,0.9, ym, 1.0-ym)
END IF
CALL xmap (0.0, xl, 0.0, yl)
CALL xqpspc(p1,p2,p3,p4)
CALL xqmap(xa,xb,ya,yb)
CALL xaxsca(0.0, xl, dx, 0.0, yl, dy)
WRITE(title,'(i1,''-D VAR NO.'',I2,'' K='',I3)')ndim,ivar,kk
CALL xchsiz(0.025*yl)
CALL xcharc(0.5*xl, 1.05*yl, title(1:21))
PRINT*
PRINT*,'Plotting ',title(1:21), ' for grid ', mptr
mptr = lstart(lfine)
100 CONTINUE
nx=node(5,mptr)
ny=node(6,mptr)
nz=node(14,mptr)
CALL xqlmsk(level)
CALL xunmsk(1)
IF(mptr /= mstart) THEN
CALL xpenup(rnode(7,mptr),rnode(8,mptr))
DO i=1,7,2
CALL xpendn(rnode(i,mptr),rnode(i+1,mptr))
END DO
END IF
CALL xrsmsk(level)
IF(ndim == 3) in1 = igtxyz(mptr,ivar1,1 )
IF(ndim == 2) in1 = igtnxy(mptr,ivar1,1 )
IF(ndim == 3) in2 = igtxyz(mptr,ivar2,1 )
IF(ndim == 2) in2 = igtnxy(mptr,ivar2,1 )
n3rd=nz
IF(ndim == 2) n3rd=1
!
needsp = nx*ny
ix=igetsp(needsp)
iy=igetsp(needsp)
iu=igetsp(needsp)
iv=igetsp(needsp)
CALL tranuv(mptr,a(in1),a(in2),a(iu),a(iv),nx,ny,n3rd,MAX(1,kk))
CALL gtacrd(mptr,ndim,iscalr,a(ix),a(iy),nxa,nya)
IF(mptr == lstart(lfine) ) THEN
zinc=1.0
mode=1
ELSE
mode=2 ! now use the zinc of first grid lstart(lfine).
END IF
IF(iplot == 1) THEN
iptr=iu
ELSE
iptr=iv
END IF
IF(mptr == mstart)THEN
CALL xconts2( zinc,mode)
ELSE
CALL xconts2( zinc,mode)
xor=rnode(1, mptr)
yor=rnode(2, mptr)
dx =rnode(9 , mptr)
dy =rnode(10, mptr)
xl=(nx-1)*dx
yl=(ny-1)*dy
cosr=rnode(21,mptr)
sinr=rnode(22,mptr)
ang = ASIN (sinr )*180/pi
xo = xor+0.5*(dx*cosr-dy*sinr)
yo = yor+0.5*(dx*sinr+dy*cosr)
CALL xrmask(xo,yo,xl-dx,yl-dx,ang)
END IF
CALL resett
mptr= lstgrd(mptr)
IF(mptr /= 0)GO TO 100
CALL xmap(x1,x2,y1,y2)
CALL xrsmsk(nmask)
RETURN
END SUBROUTINE pltvxy
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE pltgrd(mptr, ndim),14
!*********************************************************************
! plot contour maps of ndim-D fields on grid mptr. Only selectedx-y slices
! are plotted.
!*********************************************************************
INCLUDE 'nodal.inc'
INCLUDE 'agrigrid.inc'
INCLUDE 'grddsc.inc'
!*********************************************************************
! plot x-y planes of ndim-D variables on grid mptr
!
! goto 111
!*********************************************************************
ivar=1
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL plotxy(mptr,ndim,ivar,1)
CALL xframe
CALL plotxy(mptr,ndim,ivar,nzz)
CALL xframe
!
ivar=3
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL plotxy(mptr,ndim,ivar,1)
CALL xframe
CALL plotxy(mptr,ndim,ivar,nzz)
CALL xframe
ivar=5
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL plotxy(mptr,ndim,ivar,1)
CALL xframe
CALL plotxy(mptr,ndim,ivar,nzz)
CALL xframe
111 CONTINUE
ivar=7
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL plotxy(mptr,ndim,ivar,1)
CALL xframe
CALL plotxy(mptr,ndim,ivar,nzz)
CALL xframe
RETURN
ivar=9
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL plotxy(mptr,ndim,ivar,1)
CALL xframe
CALL plotxy(mptr,ndim,ivar,nzz)
CALL xframe
RETURN
ivar=11
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL plotxy(mptr,ndim,ivar,1)
CALL xframe
CALL plotxy(mptr,ndim,ivar,nzz)
CALL xframe
ivar=12
nzz = node(14,mptr)-idmxyz(3,ivar)
CALL plotxy(mptr,ndim,ivar,1)
CALL xframe
CALL plotxy(mptr,ndim,ivar,nzz)
CALL xframe
RETURN
END SUBROUTINE pltgrd
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE plotxy(mptr,ndim,ivar, kk) 14,4
!*********************************************************************
! plot x-y slices ofr ndim-d field ivar on grid mptr.
! no coordinate transformation is done.
!*********************************************************************
INCLUDE 'nodal.inc'
INCLUDE 'agrialloc.inc'
CHARACTER (LEN=132) :: title
!
nx=node(5,mptr)
ny=node(6,mptr)
nz=node(14,mptr)
pi=4.0*ATAN(1.0)
xor=rnode(1, mptr)
yor=rnode(2, mptr)
dx=rnode(9 , mptr)
dy=rnode(10, mptr)
xl=(nx-1)*dx
yl=(ny-1)*dy
CALL xqmap (x1,x2,y1,y2)
IF(yl > xl) THEN
CALL xpspac(0.1+(yl-xl)/(yl*2),0.9-(yl-xl)/(yl*2),0.1,0.8)
ELSE
CALL xpspac(0.1,0.9,0.1+(xl-yl)/(xl*2),0.9-(xl-yl)/(xl*2))
END IF
CALL xmap (0.0, xl, 0.0, yl)
CALL xaxsca(0.0, xl, dx, 0.0, yl, dy)
IF(ndim == 3) in1 = igtxyz(mptr,ivar,1 )
IF(ndim == 2) in1 = igtnxy(mptr,ivar,1 )
n3rd=nz
IF(ndim == 2) n3rd=1
!
needsp = nx*ny
ntem1=igetsp(needsp)
ntem2=igetsp(needsp)
ntem3=igetsp(needsp)
CALL gtxypl(! call gtacrd(mptr,ndim,ivar,a(ntem2),a(ntem3),nxa,nya)
CALL gtrcrd(mptr,ndim,ivar,a(ntem2),a(ntem3),nxa,nya)
CALL xconts1( WRITE(title,'(''VAR NO.'',I2,'' K='',I3,'' FOR GRID '',I2)') &
ivar, kk, mptr
CALL xchsiz(0.025*yl)
CALL xcharc(0.5*xl, 1.05*yl, title(1:27))
CALL xmap (x1,x2,y1,y2)
WRITE(title,'(''X-Y PLANE FOR K='',I2,'' OF '', &
i1,''-d var. no.'',i3,'' of grid'',i2)') kk,ndim,ivar,mptr
PRINT*,title(1:60)
! call WRIGAR(A(ntem1),1,nx,1,ny,1,1,1,nx,1,ny,1,1,
! : title(1:60),0.0 ,1)
CALL resett
RETURN
END SUBROUTINE plotxy
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE xconts1(z,x,y,md,m,n, limit) 1,1
!*********************************************************************
! plot a contour map for field Z.
! if limit=.false., max, min, inc are not plotted.
!*********************************************************************
REAL :: z(md,n),x(md,n),y(md,n),cl(150)
INTEGER :: iwrk(10000)
LOGICAL :: limit
IF(m*n > 10000)THEN
PRINT*,'Working array IW in plotting routine xconts1', &
' defined too small'
WRITE(6,'(a,i6,a,i6)')'Space needed:',m*n,', allocated:',10000
PRINT*,'No further plotting performed in this routine.'
RETURN
END IF
zmin=z(1,1)
zmax=zmin
DO j=1,n
DO i=1,m
zmax=MAX(zmax,z(i,j))
zmin=MIN(zmin,z(i,j))
END DO
END DO
PRINT*,'zmax, zmin ', zmax, zmin
cl(1)=0.0
IF(zmax-zmin > 1.0E-8)THEN
cl(2)=cl(1)+ xfinc(zmax-zmin)
IF(cl(2)-cl(1) == 0.0) cl(2)=cl(1)+1.0
mode=1
CALL xconta(z,x,y,iwrk,md,m,n,cl,ncl,mode)
ELSE
cl(2)=1.0
ncl=2
END IF
IF(.NOT.limit) RETURN
! zmax =cl(ncl)
! zmin =cl( 1)
zinc=cl(MIN(2,ncl))-cl(1)
CALL xclimt1(zmax,zmin,zinc,0.0)
RETURN
END SUBROUTINE xconts1
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE xconts2(z,x,y,md,m,n, limit, zinc0, model) 4,1
!*********************************************************************
! produce contour plot for field Z with coordinate X(i,j),Y(i,j).
! limit = .t. then zmin,zmax & zinc are plotted above the top
! model = 1, zinc is determined inside and on output zinc0 = zinc.
! = 2, zinc = zinc0.
!*********************************************************************
REAL :: z(md,n),x(md,n),y(md,n),cl(150)
INTEGER :: iwrk(10000)
LOGICAL :: limit
IF(m*n > 10000)THEN
PRINT*,'Working array IW in plotting routine xconts1', &
' defined too small'
WRITE(6,'(a,i6,a,i6)')'Space needed:',m*n,', allocated:',10000
PRINT*,'No further plotting performed in this routine.'
RETURN
END IF
zmin=z(1,1)
zmax=zmin
DO j=1,n
DO i=1,m
zmax=MAX(zmax,z(i,j))
zmin=MIN(zmin,z(i,j))
END DO
END DO
PRINT*,'zmax, zmin ', zmax, zmin
cl(1)=0.0
IF(zmax-zmin > 1.0E-8)THEN
IF(model /= 2)THEN
cl(2)=cl(1)+ xfinc(zmax-zmin)
IF(cl(2)-cl(1) == 0.0) cl(2)=cl(1)+1.0
mode=1
ELSE
cl(2)=zinc0
mode=2
END IF
CALL xconta(z,x,y,iwrk,md,m,n,cl,ncl,mode)
ELSE
cl(2)=1.0
ncl=2
END IF
! zmax =cl(ncl)
! zmin =cl( 1)
zinc=cl(MIN(2,ncl))-cl(1)
zinc0 = zinc
IF(.NOT.limit) RETURN
CALL xclimt1(zmax,zmin,zinc,0.0)
RETURN
END SUBROUTINE xconts2
!*********************************************************************
!*********************************************************************
!*********************************************************************
REAL FUNCTION xfinc(x)
!*********************************************************************
! to automatically divide domain (0,x) to a number of subdomain
! with intervale xfinc which is >=4 and =<16 for fold=1.0)
!*********************************************************************
ipower= INT( LOG10(x) )
d= INT(x/(10.0**ipower))
fold=1.0
xfinc=0.1*x
IF(d >= 0.0.AND.d < 3.0) THEN
xfinc=2.0*10.0**(ipower-1)
ELSE IF(d >= 3.0.AND.d < 7.0) THEN
xfinc=5.0*10.0**(ipower-1)*fold
ELSE IF(d >= 7.0.AND.d < 10.) THEN
xfinc=1.0*10.0** ipower*fold
END IF
IF(xfinc == 0.0) xfinc=x*0.1
RETURN
END FUNCTION xfinc
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE xclimt1(fmax,fmin,finc ) 2
!
CHARACTER (LEN=150) :: ch
CALL xqmap(xl,xr,yb,yt)
CALL xqrang( xrg, yrg )
CALL xqchmg( siz0 )
siz=0.02*SQRT( xrg*yrg)
CALL xchmag(siz)
WRITE(ch,'(''MIN='',G9.3E2,'' MAX='',G9.3E2, &
'' inc='',g9.3E2)')fmin,fmax,finc
CALL xcharc( 0.5*(xl+xr), yt+0.02*(yt-yb),ch(1:41) )
CALL xchmag( siz0)
RETURN
END SUBROUTINE xclimt1
!*********************************************************************
!*********************************************************************
!*********************************************************************
SUBROUTINE tranuv(mptr,u,v,uc,vc,nx,ny,nz,kk) 1
!*********************************************************************
! transform vector field (u,v) on grid mptr to absolute coordinate
! and copy x-y slice at k=kk to 2-D field (uc,vc).
! these vectors are located in the center of staggered grid cells.
!*********************************************************************
REAL :: u(nx,ny,nz),v(nx,ny,nz), uc(nx,ny),vc(nx,ny)
INCLUDE 'nodal.inc'
cosg=rnode(21,mptr)
sing=rnode(22,mptr)
k=MAX(kk,1)
DO j=1,ny-1
DO i=1,nx-1
uc(i,j)=(cosg*(u(i,j,k)+u(i+1,j,k))- &
sing*(v(i,j,k)+v(i,j+1,k)))*0.5
vc(i,j)=(sing*(u(i,j,k)+u(i+1,j,k))+ &
cosg*(v(i,j,k)+v(i,j+1,k)))*0.5
END DO
END DO
RETURN
END SUBROUTINE tranuv
FUNCTION abssuv(array,i0,i1,j0,j1,k0,k1)
REAL :: array(*)
lmn=(i1-i0+1)*(j1-j0+1)*(k1-k0+1)
abssuv=0.0
DO ijk=1,lmn
abssuv=abssuv+ABS(array(ijk))
END DO
RETURN
END FUNCTION abssuv