      subroutine setophys
c
c  Sets up ocean circulation model.  This is a combination of old
c  setogrd.F and initoce.F, without unused options and without
c  reading the restart files, and finally without mapping, which is
c  done elsewhere.
c      
c#include "gloparam.F.h"
#include "oparams.F.h"
#include "ogrid.F.h"
#include "ostate.F.h"
#include "mapping.h"
#include "radpth.h"
      real Wk1(len)
c
c==>  Grid data file info  
      real xse( IM+1,JM+1)     ! Long of se corner of grid box (degrees)
      real yse( IM+1,JM+1)     ! Lat of se corner of grid box (degrees)
c      real chi( IM,JM )       ! angle between coord and zonal dir (Radians)
      real Um ( IM,JM )        ! U metric term
      real Vm ( IM,JM )        ! V metric term
      
      real    zx0    ! Location of center of restart grid box (1,1)
      real    zy0    ! Location of center of restart grid box (1,1)
      real    zdx    ! Separation of restart grid boxes
      real    zdy    ! Separation of restart grid boxes
      
      real  xc(IM,JM)
      real  yc(IM,JM)
      real  dxg( IM, JM+1)
      real  dyg( IM+1,JM )
      real  dxgu( IM, JM)
      real  dygu( IM,JM )
      real  cos_chig( IM, JM)
      real  sin_chig( IM,JM )
      real  ag( IM, JM)
      real  agu( IM, JM)
      real  dx1,dx2,dx3
      real  tag

c==> Functions
      real     SSUM
      real     G_CIRCLE_DIST
      integer  ISMIN, ISMAX
c
      
c---------------------------------------------------------------------
      print *,'#################################################'
      print *,'# Start SETOPHYS: '

c
         zx0  = along0 + 0.5*dlong
         zy0  = alat0 + 0.5*dlat
         zdx  = dlong
         zdy  = dlat
         write(6,*)'zx0,zy0,zdx,zdy = ',zx0,zy0,zdx,zdy
         do j=1,JM+1
         do i=1,IM+1
          xse(i,j) = zx0 + ( i - 1.5 ) * zdx
          yse(i,j) = zy0 + ( j - 1.5 ) * zdy
         enddo
         enddo
         write(6,*)'xse1,xse2 = ',xse(1,1),xse(2,1)
         write(6,*)'yse1,yse2 = ',yse(1,1),yse(1,2)
         
         do j=1,JM
         do i=1,IM
c           chi(i,j) = 0.
           Um(i,j)  = TAN( yse(i+1,j+1) * DTOR ) / RADIUS
           Vm(i,j)  = 0.
         enddo
         enddo
        
      print *,'# setting up ocean grid'

! The input grid has lat and long of corners defined 
!  This is IM+1 by JM+1
!  Get the lengths of each side
!  Get the center long and lats
!  Get the areas of each box
!  Get the distances between centers


!  Find grid box centers from corners
!  Watch out for 360 degree longitude shifts

      do j=1,JM
      do i=1,IM
        dx1 = xse(i+1,j  ) - xse(i,j)
        dx3 = xse(i  ,j+1) - xse(i,j)
        dx2 = xse(i+1,j+1) - xse(i,j)
        dx1 = AMOD( dx1 + 540. , 360. ) - 180.
        dx2 = AMOD( dx2 + 540. , 360. ) - 180.
        dx2 = AMOD( dx3 + 540. , 360. ) - 180.
        xc(i,j) = xse(i,j) + (  dx1 + dx2 + dx3 )*.25
        yc(i,j) = (yse(i,j) + yse(i+1,j) + yse(i,j+1) + yse(i+1,j+1))*.25
      enddo
      enddo
      
      do j=1,JM
      do i=1,IM
c        cos_chig(i,j) = COS( chi(i,j) )
c        sin_chig(i,j) = SIN( chi(i,j) )
        cos_chig(i,j) = 1.0
        sin_chig(i,j) = 0.0
      enddo
      enddo
      
! Length of South box edge
      do j=1,JM+1
      do i=1,IM
       dxg(i,j) = G_CIRCLE_DIST( xse(i  ,j)*DTOR, yse(i  ,j)*DTOR,
     &                           xse(i+1,j)*DTOR, yse(i+1,j)*DTOR,
     &                           RADIUS )
      enddo
      enddo
      
      i = ISMIN( IM*(JM+1), dxg, 1 )
      j = ISMAX(IM*(JM+1), dxg, 1 )
      print *,'#   Min, max of dxg = ',dxg(i,1),dxg(j,1)
      
! Length of west box edge
      do j=1,JM
      do i=1,IM+1
       dyg(i,j) = G_CIRCLE_DIST( xse(i,j  )*DTOR, yse(i,j  )*DTOR,
     &                           xse(i,j+1)*DTOR, yse(i,j+1)*DTOR,
     &                           RADIUS )
      enddo
      enddo
      i = ISMIN( (IM+1)*JM, dyg, 1 )
      j = ISMAX( (IM+1)*JM, dyg, 1 )
      print *,'#   Min, max of dyg = ',dyg(i,1),dyg(j,1)

! EW dist between box centers 
      do j=1,JM
      do i=1,IM-1
       dxgu(i,j) = G_CIRCLE_DIST( xc(i  ,j)*DTOR, yc(i  ,j)*DTOR,
     &                            xc(i+1,j)*DTOR, yc(i+1,j)*DTOR,
     &                           RADIUS )
      enddo
      enddo
      do j=1,JM
       dxgu(IM,j) = G_CIRCLE_DIST( xc(IM ,j)*DTOR, yc(IM ,j)*DTOR,
     &                             xc(1  ,j)*DTOR, yc(1  ,j)*DTOR,
     &                           RADIUS )
      enddo
      
! NS dist between box centers 
      do j=1,JM-1
      do i=1,IM
       dygu(i,j) = G_CIRCLE_DIST( xc(i,j  )*DTOR, yc(i,j  )*DTOR,
     &                            xc(i,j+1)*DTOR, yc(i,j+1)*DTOR,
     &                           RADIUS )
      enddo
      enddo
      
! Area of grid boxes
      do j=1,JM
      do i=1,IM
        ag(i,j) = 0.25*(dxg(i,j) + dxg(i,j+1)) * (dyg(i,j)+dyg(i+1,j))      
      enddo
      enddo
      
      i = ISMIN( IM*JM, ag, 1 )
      j = ISMAX( IM*JM, ag, 1 )
      print *,'#   Min, max of ag  = ',ag(i,1),ag(j,1)

! Area of u grid boxes

      do j=1,JM-1
      do i=1,IM
        ip = MOD( i, IM ) + 1
        agu(i,j) = 0.25*(ag(i,j) + ag(i,j+1) + ag(ip,j) + ag(ip,j+1))
      enddo
      enddo
      i = ISMIN( IM*(JM-1), agu, 1 )
      j = ISMAX( IM*(JM-1), agu, 1 )
      print *,'#   Min, max of agu  = ',agu(i,1),agu(j,1)

      l = nwater+nland
      do i=1,l
       if (ieast(i).lt.1 .or. ieast(i).gt.l) then
            print *,'# Error in output from MAPPING, ieast failure:'
            print *,'#     ',i,ieast(i),' E'
            stop
       endif
       if (iwest(i).lt.1 .or. iwest(i).gt.l) then
            print *,'# Error in output from MAPPING, iwest failure:'
             print *,'#     ',i,iwest(i),' W'
            stop
       endif
       if (inorth(i).lt.1 .or. inorth(i).gt.l) then
            print *,'# Error in output from MAPPING, inorth failure:'
            print *,'#     ',i,inorth(i),' N'
            stop
       endif
       if (isouth(i).lt.1 .or. isouth(i).gt.l) then
            print *,'# Error in output from MAPPING, isouth failure:'
            print *,'#     ',i,isouth(i),' S'
            stop
       endif
       
      enddo
      
      do i=1,LEN
       Bh(i) = 0
       Bu(i) = 0
       Area_h(i) = 0.
       Area_h_r(i) = 0.0
       Area_u(i) = 0.
       Area_u_r(i) = 0.0
       oc_lat(i) = -999.
       oc_lon(i) = -999.
       oc_lat_ne(i) = -999.
       oc_lon_ne(i) = -999.
       Dx_h(i) = 0.
       Dy_h(i) = 0.
       U_metric(i)  = 0.
       V_metric(i)  = 0.
       Coriolis(i)  = 0.
       Zbottom(i) = 0.0
       cos_chi(i)   = 0.
       sin_chi(i)   = 0.
      enddo
      
            
      do i=1,nwater
       Bh(i) = 1
      enddo


      do i=1,nwater
       Bu(i) = Bh(i)*Bh(ieast(i))*Bh(inorth(i))*Bh(ieast(inorth(i)))
       oc_lon(i)    = xc ( gridi(i)+1, gridj(i)+1 )
       oc_lat(i)    = yc ( gridi(i)+1, gridj(i)+1 )
       oc_lon_ne(i) = xse( gridi(i)+2, gridj(i)+2 )
       oc_lat_ne(i) = yse( gridi(i)+2, gridj(i)+2 )
       Area_h(i)    = ag ( gridi(i)+1, gridj(i)+1 ) * Bh(i)
       Dx_h(i)      = dxg( gridi(i)+1, gridj(i)+2 )
       Dy_h(i)      = dyg( gridi(i)+2, gridj(i)+1 )
       Area_u(i)    = agu( gridi(i)+1, gridj(i)+1 ) * Bu(i)
       U_metric(i)  = um ( gridi(i)+1, gridj(i)+1 ) * Bu(i)
       V_metric(i)  = vm ( gridi(i)+1, gridj(i)+1 ) * Bu(i)
       Coriolis(i)  = 2.*OMEGA*SIN( oc_lat_ne(i)*DTOR ) * Bu(i)
       Zbottom(i)   = Depth( gridi(i)+1, gridj(i)+1 ) * Bh(i)
       cos_chi(i)   = cos_chig( gridi(i)+1, gridj(i)+1 ) * Bu(i)
       sin_chi(i)   = sin_chig( gridi(i)+1, gridj(i)+1 ) * Bu(i)
       Area_h_r(i) = Bh(i) / ( Area_h(i) + 1.E-36)
       Area_u_r(i) = Bu(i) / ( Area_u(i) + 1.E-36)
      enddo
      totalah = SSUM(nwater,Area_h,1)
      tag     = SSUM(IM*JM , Ag,   1)

! Printout grid statistics

      print *,'#  POSEIDON GRID HAS BEEN SET UP FROM RESTART FILE'
      print *,'#'
      print 1, tag
      print 2, totalah
      print 3, 100.*totalah/tag
      print 4, 100.*totalah/( 4 * PI * RADIUS * RADIUS )
      print *,'#'
      
 1    format(' #   The total area on the rectangular grid is ',1PE20.12,
     &         ' sq. meters')
 2    format(' #   The total area in water is                ',1PE20.12)
 3    format(' #   Or ',F8.3,' percent of the grid area')
 4    format(' #   Or ',F8.3,' percent of Earth surface')
 
      
      print *,'#                   done with ocean grid setup'
      print *,'#################################################'

c  Now adjust D_min and D_max to accomodate the bottom topography
 
 
      if ( EXTMODE ) then
      write (6,*) '#  USING EXTERNAL MODE VERSION'
       do i=1,nwater
        D_max(i,km) = -Zbottom(i) ! Max depth (positive down)
        D_min(i,km) = -Zbottom(i) !  Z is positive UPwards
       enddo
       do k=km-1 ,1, -1
        do i=1,nwater
         D_max(i,k) = AMIN1(D_max1d(k), D_max(i,k+1) - hmin_b(k))
         D_min(i,k) = AMIN1(D_min1d(k), D_min(i,k+1) - hmin_b(k))
        enddo
       enddo
       do k=1,km
        do i=nwater+1,LEN
         D_max(i,k) = 0.
         D_min(i,k) = 0.
        enddo
       enddo
      else ! EXTMODE
       do k=1,km
        do i=1,nwater
         D_min(i,k) = D_min1d(k)
         D_max(i,k) = D_max1d(k)
        enddo
        do i=nwater+1,len
         D_min(i,k) = 0.
         D_max(i,k) = 0.
        enddo
       enddo ! k
 
      endif ! EXTMODE
c
c      do k=1,km
c       i=ISMIN( nwater, D_min(1,k), 1)
c       j=ISMAX( nwater, D_min(1,k), 1)
c       write (6,*) '# Min, Max of D_min(',k,') =',D_min(i,k),D_min(j,k)
c       i=ISMIN( nwater, D_max(1,k), 1)
c       j=ISMAX( nwater, D_max(1,k), 1)
c       write (6,*) '# Min, Max of D_max(',k,') =',D_max(i,k),D_max(j,k)
c       do i=1,nwater
c         Wk1(i) = D_max(i,k) + Zbottom(i)
c       enddo
c       i=ISMIN( nwater, Wk1, 1)
c       j=ISMAX( nwater, Wk1, 1)
c       write (6,*) '#   Min, Max below D_max(',k,') =',Wk1(i),Wk1(j)
c      enddo
c
c  Set up initial buoyancy       
        do k=1,KM
        do i=1,nwater
         B(i,k,1) = BYNCY( T(i,k,1), S(i,k,1), 0.)
         B(i,k,2) = BYNCY( T(i,k,2), S(i,k,2), 0.)
        enddo
        do i=nwater+1,nwater+nland
         B(i,k,1) = 0.
         B(i,k,2) = 0.
        enddo
        enddo 

      call PRNTAVGS( LEN, KM, nwater, extmode,
     &               H, U, V, T, S, B, Area_h, Zbottom )
c
c  Initialize entrainment velocity
      do k = 1,km
       do i = 1,len
        W_ent(i,k) = 0.0
       enddo
      enddo
c
      print *,'                   done with ocean initialization'
      return
      end
