!#################################################################
!
! * geometry related routines *
!
!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################

module geometry

  use GO  , only : gol, goErr, goPr
  use dims, only : lm
  
  implicit none
  
  
  ! --- in/out ------------------------------------

  private

  public :: geomtryh, geomtryv, calc_dxy


contains
  
     
  !--------------------------------------------------------------------
  ! construct arrays of horizontal model geometry
  !
  !                             mh, 27-jun-1989
  !                             mh, 26-sep-1992
  !                             v 9.1
  ! Changed vertical definition to hybrid levels according to v 9.1 in new
  ! subroutine geomtryv
  !
  !                                       aj, 23-may-1995
  !                                       mk, 5-nov-1999 zoom version
  !--------------------------------------------------------------------

  subroutine geomtryh( region )

    use binas,       only: ae, pi
    use dims,        only: dx, gtor, xref, dy, yref, im, jm, ybeg, areag
    use global_data, only: region_dat

    ! --- in/out ----------------------------------

    integer, intent(in)   :: region
    
    ! --- local -----------------------------------

    real, pointer   ::  dxyp(:)
    integer         ::  j
    real            ::  dxx,dyy,lat,area
    
    ! --- begin ------------------------------------

    dxyp => region_dat(region)%dxyp

    !-------------------------------------------
    ! Standard horizontal geometry parent region (always global)
    !-------------------------------------------

    dxx = dx*gtor/xref(region)
    dyy = dy*gtor/yref(region)
    lat = ybeg(region)*gtor
    area =0.0
    do j=1,jm(region)
       dxyp(j) = dxx * (sin(lat+dyy)-sin(lat))*ae**2
       lat = lat+dyy
       area = area + dxyp(j)*im(region)
    end do

    !print *, 'geomtryh: global area ',4*pi*ae**2
    !print *, 'geomtryh: area region ',region,area

    areag(region) = area

    nullify(dxyp)

  end subroutine geomtryh


  !
  ! compute area of the model grid cells
  !

  subroutine calc_dxy( dxy, nlat )

    use binas,  only : ae, pi
    use dims,   only : gtor, nlon360
    
    ! --- in/out -----------------------------

    integer, intent(in)   ::  nlat
    real, intent(out)     ::  dxy(nlat)
    
    ! --- local -------------------------------

    real                  ::  dxx,dyy,lat,areag2
    integer               ::  j
    
    ! --- begin --------------------------------

    dxx = 1.0*gtor
    dyy = 1.0*gtor
    lat = -90.0*gtor
    areag2 =0.0
    do j=1,nlat
       dxy(j) = dxx * (sin(lat+dyy)-sin(lat))*ae**2
       lat = lat+dyy
       areag2 = areag2 + dxy(j)*nlon360
    end do
    
    !print *, 'calc_dxy: global area                   ',4*pi*ae**2
    !print *, 'calc_dxy: area calculated for 1x1 grid..',areag2

  end subroutine calc_dxy


  !--------------------------------------------------------------------
  !  define the vertical geometry of the tm model grid v9.1knmi
  !                                   aj, 30-8-1995
  !--------------------------------------------------------------------

  subroutine geomtryv(region)

    use binas      , only : grav
    use const_ec_v , only : a_ec, b_ec
    use dims       , only : echlevs, at, bt, lm

    ! --- in/out -----------------------

    integer,intent(in) :: region

    ! --- local ------------------------
    
    integer            ::  l

    ! --- begin --------------------------

    ! hybride coeff at TM levels:
    do l = 1, lm(1)+1
      at(l) = a_ec(echlevs(l-1))
      bt(l) = b_ec(echlevs(l-1))
    end do

  end subroutine geomtryv


end module geometry