
module bgrid_horiz_mod

!-----------------------------------------------------------------------
!
!        allocates memory and initializes grid constants
!              for the FMS B-grid dynamical core
!
!-----------------------------------------------------------------------

use mpp_domains_mod, only: mpp_domains_init,           &
                           mpp_define_domains,         &
                           domain1D, domain2D,         &
                           mpp_get_global_domain,      &
                           mpp_get_data_domain,        &
                           mpp_get_compute_domain,     &
                           mpp_get_compute_domains,    &
                           mpp_get_domain_components,  &
                           mpp_get_pelist,             &
                           mpp_define_layout,          &
                           mpp_get_layout,             &
                           CYCLIC_GLOBAL_DOMAIN,       &
                           mpp_global_sum,             &
                           BITWISE_EXACT_SUM

use  constants_mod, only: RADIUS
use        fms_mod, only: error_mesg, FATAL, NOTE,  &
                          mpp_pe, mpp_root_pe, mpp_npes, &
                          write_version_number, stdlog

implicit none
private

!-----------------------------------------------------------------------
!------- public interfaces -------

public  horiz_grid_init,       &
        get_horiz_grid_bound,  &
        get_horiz_grid_size,   &
        update_np, update_sp,  &
        TGRID, VGRID

!-----------------------------------------------------------------------
!             public derived data types
!-----------------------------------------------------------------------

public bgrid_type
public bgrid_interp_type
public horiz_grid_type

!-----------------------------------------------------------------------
!
!    NOTE: all horizontal indexing references global indices
!                    ( 1:nlon, 1:nlat )

type bgrid_type
   integer :: is,  ie,  js,  je         ! compute domain indices
   integer :: isd, ied, jsd, jed        ! data    domain indices
   integer :: isg, ieg, jsg, jeg        ! global  domain indices
   real, pointer, dimension(:)   :: blong=>NULL(), blatg=>NULL()  ! global grid edges
   real, pointer, dimension(:)   ::    dx=>NULL(),   rdx=>NULL(), &
                                     area=>NULL(), rarea=>NULL(), &
                                      tph=>NULL(),   tlm=>NULL()
   real, pointer, dimension(:,:) ::   aph=>NULL(),   alm=>NULL()
   real                          :: dy, rdy
   real                          :: areasum
   type(domain2D) :: Domain, Domain_nohalo
end type

!    is   = starting x-axis index for the compute domain
!    ie   = ending   x-axis index for the compute domain
!    js   = starting y-axis index for the compute domain
!    je   = ending   y-axis index for the compute domain
!    isd  = starting x-axis index for the data    domain
!    ied  = ending   x-axis index for the data    domain
!    jsd  = starting y-axis index for the data    domain
!    jed  = ending   y-axis index for the data    domain
!    isg  = starting x-axis index for the global  domain
!    ieg  = ending   x-axis index for the global  domain
!    jsg  = starting y-axis index for the global  domain
!    jeg  = ending   y-axis index for the global  domain
!
!    dx    = grid spacing for x-axis (in meters)
!    rdx   = reciprocal of dx (1/m)
!    dy    = grid spacing for y-axis (in meters)
!    rdy   = reciprocal of dy (1/m)
!    area  = area of a grid box (in m2)
!    rarea = reciprocal of area (1/m2)
!    areasum = bit-reproducible global sum of area
!
!    tph  = latitude at the center of grid box (in radians)
!    tlm  = longitude at the center of grid box (in radians)
!    aph  = actual latitude at the center of grid box (in radians)
!    alm  = actual longitude at the center of grid box (in radians)
!
!    blong = longitude grid box boundaries along the global x/longitude axis (in radians)
!    blatg = latitude  grid box boundaries along the global y/latitude  axis (in radians)
!
!    Domain        = domain2D variable with halo size = 1
!    Domain_nohalo = domain2D variable with halo size = 0
!                    used for outputing diagnostic fields
!
!-----------------------------------------------------------------------
! interpolation weights: used to move data between model grids

type bgrid_interp_type
   real, pointer, dimension(:,:) :: tmpwts=>NULL(), velwts=>NULL(), nowts=>NULL()
end type bgrid_interp_type

! tmpwts, velwts = interpolation weights for area-weighted 4-point averages
! nowts          = interpolation weights for simple 4-point averages

!-----------------------------------------------------------------------

type horiz_grid_type
   type(bgrid_type)        :: Tmp, Vel
   type(bgrid_interp_type) :: Interp
   integer :: nlon, nlat, isize, jsize
   integer :: ilb, iub, jlb, jub
   logical :: channel, double_periodic, decompx, decompy
   real    :: dlmd, dphd, dlm, dph
   real    :: sb, nb, wb, eb
   real, pointer, dimension(:,:) :: sinphv=>NULL(), tanphv=>NULL()
end type horiz_grid_type

!    Tmp = constants for the temperature/tracer/mass grid
!    Vel = constants for the u/v wind component grid
!
!    sinphv = sine of Vel%aph
!    tanphv = tangent of Vel%aph
!
!    nlon = number of grid points along the global x-axis (no halo points)
!    nlat = number of grid points along the global y-axis (no halo points)
!
!    isize = number of grid points along the x-axis for the current processor
!             (includes halo points)
!    jsize = number of grid points along the y-axis for the current processor
!             (includes halo points)
!
!    ilb  = lower bound x-axis
!    iub  = upper bound x-axis
!    jlb  = lower bound y-axis
!    jub  = upper bound y-axis
!
!    dlm  = grid spacing for x-axis (in radians)
!    dph  = grid spacing for y-axis (in radians)
!    dlmd = grid spacing for x-axis (in degrees of longitude)
!    dphd = grid spacing for y-axis (in degrees of latitude)
!
!    sb = southern boundary of temperature grid (in radians)
!    nb = northern boundary of temperature grid (in radians)
!    wb =  western boundary of temperature grid (in radians)
!    eb =  eastern boundary of temperature grid (in radians)
!
!    channel = channel model, uses grid boundaries: ed,wb,sb,nb.
!              also boundary condition at N/S walls is modified
!              may select f-plane option or double periodic options
!    double_periodic = f-plane channel with periodic boundaries in x and y.
!
!    decompx = x-axis is decomposed across more than one processor
!    decompy = y-axis is decomposed across more than one processor
!-----------------------------------------------------------------------
!-------- public parameters -------------

   integer, parameter :: TGRID = 51, VGRID = 52

!-----------------------------------------------------------------------
!-------- private data ------------

   real,    parameter :: eps=0.0001

!-------- internal parameters --------

!----- untested/unsupported options -----
!         use at you own risk 

   real :: wbd = 0.0  ! western edge in degrees of the first
                      ! longitude (i=1) of temperature grid boxes

   real :: ebd = 360.0 ! eastern edge in degrees of the last
                       ! longitude (i=nlon) of temperature grid boxes

   real :: sbd = -90.0  ! southern edge in degrees of the first
                        ! latitude row (j=1) of temperature grid boxes

   real :: nbd =  90.0  ! northern edge in degrees of the last latitude
                        ! row (j=nlat) of temperature grid boxes

   logical :: do_channel = .false. ! if TRUE, then use sbd,nbd as southern and northern 
                                   ! boundaries of a channel.

   logical :: do_fplane_approx = .false. ! f-plane channel without spherical geometry
                                         ! used in conjunction with do_channel option
   real    :: fphd = 45.                 ! coriolis latitude used for f-plane channel

   logical :: do_double_periodic = .false.  ! domain is cyclic in X and Y
                                            ! this option must be used with do_fplane_approx = true
                                            ! Tmp and Vel grids will have the same size global compute domain
!-----------------------------------------------------------------------

   real :: tph0d = 0.  ! grid/globe transformation (unsupported)
   real :: tlm0d = 0.  ! (set tph0d=tlm0d=0 for no transformation)

!-----------------------------------------------------------------------
 character(len=128) :: version = '$Id: bgrid_horiz.F90,v 19.0 2012/01/06 19:55:11 fms Exp $'
 character(len=128) :: tagname = '$Name: tikal $'
 logical :: do_vers = .true.
!-----------------------------------------------------------------------

contains

!#######################################################################

subroutine horiz_grid_init ( Hgrid, nlon, nlat, layout )

!-----------------------------------------------------------------------
!   Hgrid     = horizontal grid constants
!   nlon,nlat = the global horizontal grid resolution
!               number of longitude and latitude grid boxes, respectively
!   layout    = domain decomposition (num X pes by num Y pes),
!                default decomposition along y-axis then x-axis
!-----------------------------------------------------------------------

     type(horiz_grid_type), intent(inout) :: Hgrid
     integer,               intent(in)    :: nlon, nlat
     integer, optional,     intent(in)    :: layout(2)

!-----------------------------------------------------------------------
!------------------- local/private declarations ------------------------

real,    allocatable :: tlmi(:), tphj(:), slat(:), dxj(:), wt(:,:)
integer, allocatable :: xrows(:), yrows(:)

real    :: hpi, dtr
integer :: i, j, ilb, iub, jlb, jub, npes, pe, yflags, nlatv
integer :: is, ie, hs, he, vs, ve
integer :: isd, ied, hsd, hed, vsd, ved
integer :: isg, ieg, hsg, heg, vsg, veg
integer :: domain_layout(2)
integer :: logunit
logical :: global_x, global_y
type(domain1D) :: Domx, Domy

!------------- parallel interface --------------------------------------
!        ---- domain decomposition -----

      call mpp_domains_init

      npes = mpp_npes()
      logunit = stdlog()

! write version info to logfile
     if (do_vers) then
       call write_version_number(version, tagname)
       do_vers = .false.
     endif

! error checks and messages

   Hgrid % channel = do_channel
   Hgrid % double_periodic = do_double_periodic

   if (Hgrid%channel) then
       write (logunit,'(a)') 'Channel model option has been selected.'
   endif

   if (do_fplane_approx) then
       if (.not.Hgrid%channel) call error_mesg ('horiz_grid_init', &
             'f-plane approximation cannot be used without do_channel=TRUE', FATAL)
       write (logunit,'(a)') '... the f-plane approximation will be used'
   endif

   if (do_double_periodic) then
       ! error check
       if (.not.do_fplane_approx) call error_mesg ('horiz_grid_init', &
              'double periodic option cannot be used without do_fplane_approx=TRUE', FATAL)
       yflags = CYCLIC_GLOBAL_DOMAIN
       nlatv  = nlat
       write (logunit,'(a)') '... and double periodic boundary conditions will be used'
   else
       yflags = 0
       nlatv  = nlat-1
   endif

! print channel options used
   if (Hgrid%channel) then
       write (logunit,'(4x,a4,f10.3,4x,a4,f10.3)') 'WBD=',wbd,'EBD=',ebd
       write (logunit,'(4x,a4,f10.3,4x,a4,f10.3)') 'SBD=',sbd,'NBD=',nbd
       if (do_fplane_approx) then
           write (logunit,'(4x,a4,f10.3)') 'FPHD=',fphd
       endif
   endif
 
!---- set-up x- & y-axis decomposition -----

 domain_layout = (/ 0, 0 /)
 if (present(layout)) domain_layout = layout
 if (domain_layout(1)+domain_layout(2) == 0) then
     call mpp_define_layout ( (/1,nlon,1,nlat/), npes, domain_layout )
 else
     if (domain_layout(1) == 0) domain_layout(1) = npes/domain_layout(2)
     if (domain_layout(2) == 0) domain_layout(2) = npes/domain_layout(1)
 endif

  if ( domain_layout(1)*domain_layout(2) /= npes ) call error_mesg &
              ('horiz_grid_init', 'number of processors requested not &
                                  &compatible with grid', FATAL )

! flag to indicate axis decomposition
  Hgrid % decompx = domain_layout(1) .gt. 1
  Hgrid % decompy = domain_layout(2) .gt. 1


!    ---- mass/temperature grid domain with halo = 0,1 ----

     call mpp_define_domains ( (/1,nlon,1,nlat/), domain_layout,   &
                               Hgrid % Tmp % Domain,               &
                               xflags = CYCLIC_GLOBAL_DOMAIN,      &
                               yflags = yflags,                    &
                               xhalo = 1, yhalo = 1,               &
          name = 'ATMOSPHERIC (B-GRID) MODEL, temperature grid, halo=1,')
     call mpp_define_domains ( (/1,nlon,1,nlat/), domain_layout,   &
                               Hgrid % Tmp % Domain_nohalo,        &
                               xflags = CYCLIC_GLOBAL_DOMAIN,      &
                               yflags = yflags,                    &
                               xhalo = 0, yhalo = 0                )

!   ---- compute exact decomposition ----
!   ---- compute 2d layout of PEs ----

     allocate ( xrows(domain_layout(1)), yrows(domain_layout(2)) )
     call mpp_get_domain_components ( Hgrid%Tmp%Domain, Domx, Domy )
     call mpp_get_compute_domains   ( Domx, size=xrows )
     call mpp_get_compute_domains   ( Domy, size=yrows )

!    ---- velocity grid may have one less latitude row ----

     if (nlatv == nlat-1) then
         yrows(domain_layout(2)) = yrows(domain_layout(2)) - 1
     endif

!    ---- velocity grid domain with halo = 0,1 ----

     call mpp_define_domains ( (/1,nlon,1,nlatv/), domain_layout,  &
                               Hgrid % Vel % Domain,               &
                               xflags = CYCLIC_GLOBAL_DOMAIN,      &
                               yflags = yflags,                    &
                               xhalo = 1, yhalo = 1,               &
                               xextent = xrows, yextent = yrows    )
     call mpp_define_domains ( (/1,nlon,1,nlatv/), domain_layout,  &
                               Hgrid % Vel % Domain_nohalo,        &
                               xflags = CYCLIC_GLOBAL_DOMAIN,      &
                               yflags = yflags,                    &
                               xhalo = 0, yhalo = 0,               &
                               xextent = xrows, yextent = yrows    )

     deallocate ( xrows, yrows )


!------------- indices for global compute domain -----------------------

     call mpp_get_global_domain ( Hgrid%Tmp%Domain, isg, ieg, hsg, heg )
     call mpp_get_global_domain ( Hgrid%Vel%Domain, isg, ieg, vsg, veg )

     Hgrid % Tmp % isg = isg;   Hgrid % Tmp % ieg = ieg
     Hgrid % Tmp % jsg = hsg;   Hgrid % Tmp % jeg = heg
     Hgrid % Vel % isg = isg;   Hgrid % Vel % ieg = ieg
     Hgrid % Vel % jsg = vsg;   Hgrid % Vel % jeg = veg

!------------- indices for data domain -----------------------

     call mpp_get_data_domain ( Hgrid%Tmp%Domain, isd, ied, hsd, hed )
     call mpp_get_data_domain ( Hgrid%Vel%Domain, isd, ied, vsd, ved )

     Hgrid % Tmp % isd = isd;   Hgrid % Tmp % ied = ied
     Hgrid % Tmp % jsd = hsd;   Hgrid % Tmp % jed = hed
     Hgrid % Vel % isd = isd;   Hgrid % Vel % ied = ied
     Hgrid % Vel % jsd = vsd;   Hgrid % Vel % jed = ved

!------------- indices for computational domain ------------------------

     call mpp_get_compute_domain ( Hgrid%Tmp%Domain, is, ie, hs, he )
     call mpp_get_compute_domain ( Hgrid%Vel%Domain, is, ie, vs, ve )

     Hgrid % Tmp % is = is;   Hgrid % Tmp % ie = ie
     Hgrid % Tmp % js = hs;   Hgrid % Tmp % je = he
     Hgrid % Vel % is = is;   Hgrid % Vel % ie = ie
     Hgrid % Vel % js = vs;   Hgrid % Vel % je = ve

!------------- indices including halo regions --------------------------

     call mpp_get_data_domain ( Hgrid%Tmp%Domain, ilb, iub, jlb, jub )

     Hgrid % ilb = ilb;  Hgrid % iub = iub
     Hgrid % jlb = jlb;  Hgrid % jub = jub

!------------ other values ---------------------------------------------

     Hgrid % nlon = nlon  ! global resolution
     Hgrid % nlat = nlat

     Hgrid % isize = Hgrid % iub - Hgrid % ilb + 1  ! array size for current PE
     Hgrid % jsize = Hgrid % jub - Hgrid % jlb + 1

!-----------------------------------------------------------------------
!    -------- allocate space for arrays ---------

    call alloc_array_space ( ilb, iub, jlb, jub, Hgrid%Tmp )
    call alloc_array_space ( ilb, iub, jlb, jub, Hgrid%Vel )
                 
    allocate ( Hgrid % sinphv(ilb:iub,jlb:jub), &
               Hgrid % tanphv(ilb:iub,jlb:jub)  )

    allocate ( slat (hs-1:he+2), wt (ilb:iub,jlb:jub) ) 

!-----------------------------------------------------------------------
!--------------derived geometrical constants----------------------------

      hpi = acos(0.0)
      dtr = hpi/90.

      if (Hgrid%channel) then
          Hgrid % dphd = (nbd-sbd)/real(nlat)
          Hgrid % dlmd = (ebd-wbd)/real(nlon)
          Hgrid % sb = sbd*dtr
          Hgrid % nb = nbd*dtr
          Hgrid % wb = wbd*dtr
          Hgrid % eb = ebd*dtr
          global_x = abs(ebd-wbd-360.).lt.eps
          global_y = abs(nbd-sbd-180.).lt.eps
      else
          ! global code used for bit-reprocibility
          Hgrid % dphd = 180./real(nlat)
          Hgrid % dlmd = 360./real(nlon)
          Hgrid % sb =  -hpi
          Hgrid % nb =  +hpi
          Hgrid % wb =   0.0
          Hgrid % eb = 360.0*dtr
          global_x = .true.
          global_y = .true.
      endif

      Hgrid % dlm = Hgrid % dlmd*dtr
      Hgrid % dph = Hgrid % dphd*dtr

!     --- dy is the same on both grids ---

      Hgrid % Tmp %  dy = RADIUS * Hgrid % dph
      Hgrid % Vel %  dy = RADIUS * Hgrid % dph
      Hgrid % Tmp % rdy = 1.0 / Hgrid % Tmp % dy
      Hgrid % Vel % rdy = 1.0 / Hgrid % Vel % dy

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!---- get grid box boundaries for computational global domain ----

     Hgrid % Tmp % blong(isg-1:ieg+2) = get_grid (nlon+3, Hgrid%wb-Hgrid%dlm, &
                                                                   Hgrid%dlm  )
   ! global cyclic continuity (may not be necessary)
     if (global_x) then
        Hgrid % Tmp % blong(isg-1) = Hgrid % Tmp % blong(ieg)-4.*hpi
        Hgrid % Tmp % blong(ieg+1) = Hgrid % Tmp % blong(isg)+4.*hpi
        Hgrid % Tmp % blong(ieg+2) = Hgrid % Tmp % blong(isg+1)+4.*hpi
     endif

     Hgrid % Tmp % blatg(hsg+1:heg) = get_grid (nlat-1, Hgrid%sb+Hgrid%dph, &
                                                                 Hgrid%dph  )
   ! polar halos
     Hgrid % Tmp % blatg(heg+2) = Hgrid % Tmp % blatg(heg)
     Hgrid % Tmp % blatg(heg+1) = Hgrid % nb
     Hgrid % Tmp % blatg(hsg)   = Hgrid % sb
     Hgrid % Tmp % blatg(hsg-1) = Hgrid % Tmp % blatg(hsg+1)

     do i = is-1, ie+1
        Hgrid % Tmp % tlm(i) = 0.5*(Hgrid%Tmp%blong(i)+Hgrid%Tmp%blong(i+1))
     enddo

     do j = hs-1, he+1
        Hgrid % Tmp % tph(j) = 0.5*(Hgrid%Tmp%blatg(j)+Hgrid%Tmp%blatg(j+1))
     enddo

     if (do_fplane_approx) then
!    --- no spherical geometry ---
         Hgrid % Tmp % dx(hs-1:he+1) = RADIUS * Hgrid % dlm
     else
!    --- sphere ---
         slat(hs-1:he+2) = sin(Hgrid%Tmp%blatg(hs-1:he+2))
         Hgrid % Tmp % dx (hs-1:he+1) = RADIUS * Hgrid % dlm / Hgrid % dph *  &
                                        abs(slat(hs:he+2)-slat(hs-1:he+1))
     endif

!-------------initialize lat/lon at velocity points---------------------

     Hgrid % Vel % blong(isg-1:ieg+2) = get_grid (nlon+3, Hgrid%wb-0.5*Hgrid%dlm, &
                                                                       Hgrid%dlm  )
   ! global cyclic continuity
     if (global_x) then
        Hgrid % Vel % blong(isg-1) = Hgrid % Vel % blong(ieg)-4.*hpi
        Hgrid % Vel % blong(ieg+1) = Hgrid % Vel % blong(isg)+4.*hpi
        Hgrid % Vel % blong(ieg+2) = Hgrid % Vel % blong(isg+1)+4.*hpi
     endif

     Hgrid % Vel % blatg(vsg:veg+1) = get_grid (nlat, Hgrid%sb+0.5*Hgrid%dph, &
                                                                   Hgrid%dph  )
   ! polar halos
     Hgrid % Vel % blatg(vsg-1) = Hgrid % sb
     Hgrid % Vel % blatg(veg+2) = Hgrid % nb

     do i = is-1, ie+1
        Hgrid % Vel % tlm(i) = 0.5*(Hgrid%Vel%blong(i)+Hgrid%Vel%blong(i+1))
     enddo
     do j = vs-1, ve+1
        Hgrid % Vel % tph(j) = 0.5*(Hgrid%Vel%blatg(j)+Hgrid%Vel%blatg(j+1))
     enddo

     if (do_fplane_approx) then
!    --- no spherical geometry ---
         Hgrid % Vel % dx = RADIUS * Hgrid % dlm
     else
!    --- sphere ---
         slat(vs-1:ve+2) = sin(Hgrid%Vel%blatg(vs-1:ve+2))
         Hgrid % Vel % dx (vs-1:ve+1) = RADIUS * Hgrid % dlm / Hgrid % dph *  &
                                        (slat(vs:ve+2)-slat(vs-1:ve+1))
     endif
     Hgrid % Vel % rdx(vs-1:ve+1) = 1.0 / Hgrid % Vel % dx(vs-1:ve+1)

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!------------- grid box areas ------------------------------------------

      Hgrid % Tmp % area = Hgrid % Tmp % dx * Hgrid % Tmp % dy
      Hgrid % Vel % area = Hgrid % Vel % dx * Hgrid % Vel % dy

! Note: area of the pole on velocity grid is not important because:
!         1) pole is not in compute domain
!         2) current polar boundary condition is u=v=0
!         3) fluxes between pole and sub-pole row are set to zero

!--- reciprocal of area ----

      do j = jlb, jub
         if (Hgrid % Tmp % area(j) > 0.0) Hgrid % Tmp % rarea(j) = 1.0 / Hgrid % Tmp % area(j)
         if (Hgrid % Vel % area(j) > 0.0) Hgrid % Vel % rarea(j) = 1.0 / Hgrid % Vel % area(j)
      enddo

!--- compute bit-reproducible global sum of area ---

      do j = hsd, hed
         wt(:,j) = Hgrid%Tmp%area(j)
      enddo
      Hgrid % Tmp % areasum = mpp_global_sum (Hgrid%Tmp%Domain, wt, flags=BITWISE_EXACT_SUM)

      do j = vsd, ved  ! may not be full array extent
         wt(:,j) = Hgrid%Vel%area(j)
      enddo
      Hgrid % Vel % areasum = mpp_global_sum (Hgrid%Vel%Domain, wt(:,vsd:ved), flags=BITWISE_EXACT_SUM)

      deallocate ( slat, wt )


!--- initialization of weights for grid interpolation ---

     call bgrid_interp_init ( Hgrid )

!-----------------------------------------------------------------------
! unsupported option for transforming the position of poles
! tph0d,tlm0d are the lat,lon position in the transformed grid
! of the point [lat=0,lon=0] ..... I think ????

      if (tph0d > eps .or. tlm0d > eps) then
         ! compute "actual" lat/lon ( aph, aph )
         ! actual lat/lon would be used for the coriolis, radiation, and ...
         ! at temperature points
           call trans_latlon (tph0d,  tlm0d,                        &
                              Hgrid % Tmp % tlm, Hgrid % Tmp % tph, &
                              Hgrid % Tmp % alm, Hgrid % Tmp % aph)
         ! at velocity points
           call trans_latlon (tph0d,  tlm0d,                        &
                              Hgrid % Vel % tlm, Hgrid % Vel % tph, &
                              Hgrid % Vel % alm, Hgrid % Vel % aph)
      else
         ! no transformation
           do j = jlb, jub
           do i = ilb, iub
              Hgrid % Tmp % aph(i,j) = Hgrid % Tmp % tph(j)
              Hgrid % Vel % aph(i,j) = Hgrid % Vel % tph(j)
              Hgrid % Tmp % alm(i,j) = Hgrid % Tmp % tlm(i)
              Hgrid % Vel % alm(i,j) = Hgrid % Vel % tlm(i)
           enddo
           enddo
      endif

!------------- trigometric constants -----------------------------------

      if (do_fplane_approx) then
!     --- channel model ---
          Hgrid % sinphv = sin(fphd*dtr)
          Hgrid % tanphv = tan(fphd*dtr)
      else
!     --- sphere ---
          Hgrid % sinphv = sin(Hgrid % Vel % aph)
          Hgrid % tanphv = tan(Hgrid % Vel % aph)
      endif

!-----------------------------------------------------------------------

end subroutine horiz_grid_init

!##############################################################################
! initializes weights for interpolation between model grids

subroutine bgrid_interp_init ( Hgrid )
type(horiz_grid_type), intent(inout) :: Hgrid

real :: hpi
real :: sph  (1:2*Hgrid%nlat+1), area(-1:2*Hgrid%nlat+2)
real :: areah(Hgrid%Tmp%jsd:Hgrid%Tmp%jed), areav(Hgrid%Vel%jsd:Hgrid%Vel%jed)
integer :: i, j, k
! indexing for interpolation weights array
  integer, parameter :: SOUTH = 1, NORTH = 2

! allocate memory
 allocate ( Hgrid%Interp%tmpwts(Hgrid%jlb:Hgrid%jub,2), &
            Hgrid%Interp%velwts(Hgrid%jlb:Hgrid%jub,2), &
            Hgrid%Interp%nowts (Hgrid%jlb:Hgrid%jub,2)  )

! weights for simple 4-pt averages with equal weighting
  Hgrid%Interp%nowts = 0.25

! equal weighting for non-spherical grid
  if (do_fplane_approx) then
      Hgrid%Interp%tmpwts = 0.25
      Hgrid%Interp%velwts = 0.25
      return
  endif

! create weights that will conserve the quantity being interpolated

! create areas for grid boxes and interpolation
! need sin(lat) every 1/2 delta-lat
  hpi = acos(0.0)
  if (.not.Hgrid%channel) then
      ! original global code (used for bit-reproducibility)
      sph(1) = -1.; sph(2*Hgrid%nlat+1) = 1.
      do j = 2, 2*Hgrid%nlat
        sph(j) = sin(-hpi+real(j-1)*Hgrid%dph*0.5)
      enddo
  else
      sph(1) = sin(Hgrid%sb); sph(2*Hgrid%nlat+1) = sin(Hgrid%nb)
      do j = 2, 2*Hgrid%nlat
        sph(j) = sin(Hgrid%sb+real(j-1)*Hgrid%dph*0.5)
      enddo
  endif
! 1/4 area/radius^2 of grid boxes
  area(1:2*Hgrid%nlat) = 0.5*Hgrid%dlm*(sph(2:2*Hgrid%nlat+1)-sph(1:2*Hgrid%nlat))
! assume halo size of one (fill two 1/2 size grid boxes)
  area(-1:0) = area(2:1:-1)
  area(2*Hgrid%nlat+1:2*Hgrid%nlat+2) = area(2*Hgrid%nlat:2*Hgrid%nlat-1:-1)

! area of temperature grid boxes
  do j = Hgrid%Tmp%jsd, Hgrid%Tmp%jed
    Hgrid%Interp%tmpwts(j,SOUTH) = area(2*j-1)
    Hgrid%Interp%tmpwts(j,NORTH) = area(2*j)
    areah(j) = (Hgrid%Interp%tmpwts(j,SOUTH)+Hgrid%Interp%tmpwts(j,NORTH))*2.0
    ! compute final weights for 4-pt averages
    Hgrid%Interp%tmpwts(j,SOUTH) = Hgrid%Interp%tmpwts(j,SOUTH) / areah(j)
    Hgrid%Interp%tmpwts(j,NORTH) = Hgrid%Interp%tmpwts(j,NORTH) / areah(j)
  enddo

! area of velocity grid boxes
  do j = Hgrid%Vel%jsd, Hgrid%Vel%jed
    Hgrid%Interp%velwts(j,SOUTH) = area(2*j)
    Hgrid%Interp%velwts(j,NORTH) = area(2*j+1)
    areav(j) = (Hgrid%Interp%velwts(j,SOUTH)+Hgrid%Interp%velwts(j,NORTH))*2.0
    ! compute final weights for 4-pt averages
    Hgrid%Interp%velwts(j,SOUTH) = Hgrid%Interp%velwts(j,SOUTH) / areav(j)
    Hgrid%Interp%velwts(j,NORTH) = Hgrid%Interp%velwts(j,NORTH) / areav(j)
  enddo

end subroutine bgrid_interp_init

!#######################################################################

function get_grid (npts, start, space) result (grid)

integer, intent(in) :: npts
real,    intent(in) :: start, space
real                :: grid(npts)
integer :: j

!---- compute equally spaced grid ----

      do j = 1, npts
         grid(j) = start + real(j-1)*space
      enddo

end function get_grid

!#######################################################################

   subroutine trans_latlon (tph0d, tlm0d, tlm, tph, alm, aph)

   real, intent(in)  :: tph0d, tlm0d, tlm(:), tph(:)
   real, intent(out) :: alm(:,:), aph(:,:)

   real :: dtr, pi, tph0, tlm0, stph0, ctph0, ttph0, cc, ee
   real, dimension(size(tph)) :: stph, ctph
   integer :: i, j

!-----------------------------------------------------------------------

    ! scalars
      pi  = 2.*acos(0.0)
      dtr  = pi/180.
      tph0 = tph0d*dtr;  tlm0 = tlm0d*dtr
      stph0 = sin(tph0); ctph0 = cos(tph0)
      ttph0 = stph0/ctph0

!     ------- compute actual lat and lon -------
      stph = sin(tph)
      ctph = cos(tph)
      do j = 1, size(aph,2)
      do i = 1, size(aph,1)
         cc = ctph(j)*cos(tlm(i))
         aph(i,j) = asin(ctph0*stph(j)+stph0*cc)
         ee = cc/(ctph0*cos(aph(i,j))) - tan(aph(i,j))*ttph0
         ee = min(ee,1.)
         if (tlm(i) > pi) then
             alm(i,j) = tlm0-acos(ee)
         else
             alm(i,j) = tlm0+acos(ee)
         endif
      enddo
      enddo

!-----------------------------------------------------------------------

end subroutine trans_latlon

!#######################################################################

subroutine get_horiz_grid_bound ( Hgrid, grid, blon, blat, global )

!       returns the grid box boundaries for either
!           the compute or global grid
!
!   Hgrid = horizontal grid constants
!   grid  = grid identifier, possible values: TGRID, VGRID
!   blon  = longitude edges in radians
!   blat  = latitude edges in radians
!   global = values for compute(F) or global(T) grid?

   type (horiz_grid_type), intent(in)  :: Hgrid
   integer,                intent(in)  ::  grid
   real,                   intent(out) :: blon(:), blat(:)
   logical, optional,      intent(in)  :: global

    select case (grid)
       case (TGRID)
          call horiz_grid_bound ( Hgrid%Tmp, blon, blat, global )
       case (VGRID)
          call horiz_grid_bound ( Hgrid%Vel, blon, blat, global )
       case default
          call error_mesg ('get_horiz_grid_bound', 'invalid grid', FATAL)
    end select

!-----------------------------------------------------------------------

end subroutine get_horiz_grid_bound

!#######################################################################

subroutine get_horiz_grid_size ( Hgrid, grid, nlon, nlat, global )

!     returns the number of longitude and latitude grid boxes
!            for either the compute or global grid
!
!   Hgrid = horizontal grid constants
!   grid  = grid identifier, possible values: TGRID, VGRID
!   nlon  = number longitude grid boxes
!   nlat  = number latitude grid boxes
!   global = values for compute(F) or global(T) grid?

   type (horiz_grid_type), intent(in)  :: Hgrid
   integer,                intent(in)  ::  grid
   integer,                intent(out) :: nlon, nlat
   logical, optional,      intent(in)  :: global

    select case (grid)
       case (TGRID)
          call horiz_grid_size ( Hgrid%Tmp, nlon, nlat, global )
       case (VGRID)
          call horiz_grid_size ( Hgrid%Vel, nlon, nlat, global )
       case default
          call error_mesg ('get_horiz_grid_size', 'invalid grid', FATAL)
    end select

end subroutine get_horiz_grid_size

!#######################################################################

subroutine horiz_grid_bound ( Grid, blon, blat, global )

   type (bgrid_type), intent(in)  :: Grid
   real,              intent(out) :: blon(:), blat(:) 
   logical, optional, intent(in)  :: global

!      private routine that returns the grid box boundaries
!          for either the compute or global grid

    integer :: is, ie, js, je
    logical :: lglobal 

    lglobal = .false.;  if (present(global)) lglobal = global


    if (lglobal) then
      is = Grid % isg; ie = Grid % ieg  ! global grid
      js = Grid % jsg; je = Grid % jeg
    else    
      is = Grid % is ; ie = Grid % ie   ! compute grid
      js = Grid % js ; je = Grid % je
    endif   

     !----- define longitudinal grid box edges -----

      if (size(blon) /= ie-is+2) call error_mesg  &
                         ('get_horiz_grid_bound', &
                          'invalid argument dimension for blon', FATAL)

      blon = Grid % blong (is:ie+1)

     !----- define latitudinal grid box edges -----

      if (size(blat) /= je-js+2) call error_mesg  &
                         ('get_horiz_grid_bound', &
                          'invalid argument dimension for blat', FATAL)

      blat = Grid %  blatg (js:je+1)

!-----------------------------------------------------------------------

end subroutine horiz_grid_bound

!#######################################################################

subroutine horiz_grid_size ( Grid, nlon, nlat, global )

   type (bgrid_type), intent(in)  :: Grid
   integer,                intent(out) :: nlon, nlat
   logical, optional,      intent(in)  :: global

   logical :: lglobal

   lglobal = .false.;  if (present(global)) lglobal = global

  !---- return the size of the requested grid ----

   if (lglobal) then
      nlon = Grid % ieg - Grid % isg + 1   ! global grid
      nlat = Grid % jeg - Grid % jsg + 1
   else
      nlon = Grid % ie - Grid % is + 1     ! compute grid
      nlat = Grid % je - Grid % js + 1
   endif

end subroutine horiz_grid_size

!#######################################################################

function update_np (Hgrid, grid) result (answer)

!   Hgrid = horizontal grid constants
!   grid  = grid identifier, possible values: TGRID, VGRID
!
!   Returns TRUE if the northernmost latitude row
!   is adjacent to the north pole.

  type (horiz_grid_type), intent(in)  :: Hgrid
  integer,                intent(in)  ::  grid
  logical :: answer
  integer :: je, jeg

     answer = .false.
    !if (Hgrid%Tmp%jeg == Hgrid%Vel%jeg) return ! periodic in y
     if (Hgrid%double_periodic) return ! periodic in y

     select case (grid)
       case (TGRID)
          if ( Hgrid%Tmp%je == Hgrid%Tmp%jeg ) answer = .true.
       case(VGRID)
          if ( Hgrid%Vel%je == Hgrid%Vel%jeg ) answer = .true.
       case default
          call error_mesg ('update_np', 'invalid grid', FATAL)
     end select

end function update_np

!#######################################################################

function update_sp (Hgrid, grid) result (answer)

!   Hgrid = horizontal grid constants
!   grid  = grid identifier, possible values: TGRID, VGRID
!
!   Returns TRUE if the southernmost latitude row
!   is adjacent to the south pole.

  type (horiz_grid_type), intent(in)  :: Hgrid
  integer,                intent(in)  ::  grid
  logical :: answer
  integer :: js, jsg

     answer = .false.
    !if (Hgrid%Tmp%jeg == Hgrid%Vel%jeg) return ! periodic in y
     if (Hgrid%double_periodic) return ! periodic in y

     select case (grid)
       case (TGRID)
          if ( Hgrid%Tmp%js == Hgrid%Tmp%jsg ) answer = .true.
       case(VGRID)
          if ( Hgrid%Vel%js == Hgrid%Vel%jsg ) answer = .true.
       case default
          call error_mesg ('update_sp', 'invalid grid', FATAL)
     end select

end function update_sp

!#######################################################################

 subroutine alloc_array_space ( ilb, iub, jlb, jub, Grid )
 integer, intent(in) :: ilb, iub, jlb, jub
 type(bgrid_type), intent(inout) :: Grid
         
      allocate ( Grid % dx   (jlb:jub), &
                 Grid % rdx  (jlb:jub), &
                 Grid % tph  (jlb:jub), &
                 Grid % tlm  (ilb:iub), &
                 Grid % area (jlb:jub), &
                 Grid % rarea(jlb:jub), &
                 Grid % aph  (ilb:iub,jlb:jub), &
                 Grid % alm  (ilb:iub,jlb:jub)  )

      allocate ( Grid % blong (Grid%isg-1:Grid%ieg+2), &
                 Grid % blatg (Grid%jsg-1:Grid%jeg+2)  )

      Grid % dx   = 0.; Grid % rdx   = 0.
      Grid % tph  = 0.; Grid % tlm   = 0.
      Grid % area = 0.; Grid % rarea = 0.

 end subroutine alloc_array_space

!#######################################################################

end module bgrid_horiz_mod

