!### 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"
!
!---------------------------------------------------------------------------
!                    TM5-MP                                                   !
!---------------------------------------------------------------------------
!BOP
!
! !MODULE:      GEOMETRY
!
! !DESCRIPTION: geometry related routines.
!\\
!\\
! !INTERFACE: 
!
MODULE GEOMETRY
  !
  ! !USES:
  !
  use GO  , only : gol, goErr, goPr
  use dims, only : lm

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: geomtryv
  public :: geomtryh
  public :: calc_dxy
  !
  ! !REMARKS:
  !
  !EOP
  !-----------------------------------------------------------------------

CONTAINS

  !-----------------------------------------------------------------------
  !                    TM5-MP                                               !
  !-----------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    GEOMTRYH
  !
  ! !DESCRIPTION: fill area variables: region_dat(region)%dxyp and areag
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE GEOMTRYH( region )
    !
    ! !USES:
    !
    use binas,         only: ae, pi
    use dims,          only: dx, gtor, xref, dy, yref, im, jm, ybeg, areag
    use global_data,   only: region_dat
    use tm5_distgrid,  only: dgrid, Get_DistGrid
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)   :: region
    !
    ! !REVISION HISTORY: 
    !    mh, 27-jun-1989                        
    !    mh, 26-sep-1992                       
    !    aj, 23-may-1995             
    !    mk, 5-nov-1999 zoom version 
    !    9 Nov 2011 - P. Le Sager - adapted for TM5-MP
    !
    ! !REMARKS:
    ! (1) areag is never used, but saved in several output files
    ! (2) FIXME ZOOM : NOT TESTED FOR REGION > 1 
    !
    !EOP
    !---------------------------------------------------------------------
    !BOC

    real, pointer   ::  dxyp(:)
    integer         ::  j, i0, i1, j0, j1
    real            ::  dxx,dyy,lat,area, deltaX, yb

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

    dxyp => region_dat(region)%dxyp

    call Get_DistGrid( dgrid(region), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 )
    
    !-------------------------------------------
    ! Standard horizontal geometry parent region (always global)
    !-------------------------------------------

    dxx = dx*gtor/xref(region)
    dyy = dy*gtor/yref(region)

    deltaX = (i1-i0+1)

    ! Cannot do that for bitwise reproductibility :
    !    yb = ybeg(region) + ( j0 - 1 ) * dy
    !    lat = yb*gtor
    ! Need to loop globally, and start at j=1:

    lat = ybeg(region)*gtor
    area =0.0

    do j=1,jm(region)
       if (j>=j0.and.j<=j1) then
          dxyp(j) = dxx * (sin(lat+dyy)-sin(lat))*ae**2
          area = area + dxyp(j)*deltaX
       end if
       lat = lat+dyy
    end do

    areag(region) = area

    nullify(dxyp)

  END SUBROUTINE GEOMTRYH
  !EOC

  !-------------------------------------------------------------------------
  !                    TM5                                                 !
  !-------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CALC_DXY
  !
  ! !DESCRIPTION: for a 1x1 grid, covering the [-90., -90+nlat] latitude
  !               range, compute the area of each grid cells.
  ! 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CALC_DXY( dxy, nlat )
    !
    ! !USES:
    !
    use binas,  only : ae, pi
    use dims,   only : gtor, nlon360
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)   ::  nlat      ! number of 1 degree zonal bands
    !
    ! !OUTPUT PARAMETERS:
    !
    real, intent(out)     ::  dxy(nlat) ! area for each zonal band
    !
    ! !REVISION HISTORY: 
    !    9 Nov 2011 - P. Le Sager -
    !
    ! !REMARKS:
    ! (1) this is called only once in the initexit/start routine, and with
    !       nlat=180, to fill the dims:dxy11 variable
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real      ::  dxx, dyy, lat
    integer   ::  j

    dxx = 1.0*gtor
    dyy = 1.0*gtor
    lat = -90.0*gtor

    do j=1,nlat
       dxy(j) = dxx * (sin(lat+dyy)-sin(lat))*ae**2
       lat = lat+dyy
    end do

  END SUBROUTINE CALC_DXY
  !EOC

  !-------------------------------------------------------
  !                    TM5                               !
  !-------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    GEOMTRYV
  !
  ! !DESCRIPTION: define the vertical geometry of the tm 
  !               model grid v9.1knmi (aj, 30-8-1995) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE GEOMTRYV()
    !
    ! !USES:
    !
    use binas      , only : grav
    use const_ec_v , only : a_ec, b_ec
    use dims       , only : echlevs, at, bt, lm
    !
    ! !REVISION HISTORY: 
    !
    ! !REMARKS:
    ! (1) used only by meteo_init_grid
    !
    !EOP
    !-----------------------------------------------------
    !BOC
    integer ::  l

    ! 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
  !EOC

END MODULE GEOMETRY