SUBROUTINE ana_pair (ng, tile, model)
!
!! svn $Id: $
!!======================================================================
!! Copyright (c) 2002-2007 The ROMS/TOMS Group                         !
!!   Licensed under a MIT/X style license                              !
!!   See License_ROMS.txt                                              !
!!                                                                     !
!=======================================================================
!                                                                      !
!  This routine sets surface air pressure (mb) using an analytical     !
!  expression.                                                         !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_forces
      USE mod_ncparam
# ifdef RISSAGA
!  Need grid coordinates in meters and pi
      USE mod_grid
      USE mod_scalars
# endif
!
! Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model

#include "tile.h"
!
      CALL ana_pair_tile (ng, model, Istr, Iend, Jstr, Jend,            &
     &                    LBi, UBi, LBj, UBj,                           &
# ifdef RISSAGA
     &                    GRID(ng) % xr, GRID(ng) % yr,                 &
# endif
     &                    FORCES(ng) % Pair)
!
! Set analytical header file name used.
!
      IF (Lanafile) THEN
        ANANAME(17)='ROMS/Functionals/ana_pair.h'
      END IF

      RETURN
      END SUBROUTINE ana_pair
!
!***********************************************************************
      SUBROUTINE ana_pair_tile (ng, model, Istr, Iend, Jstr, Jend,      &
     &                          LBi, UBi, LBj, UBj,                     &
# ifdef RISSAGA
     &                          xr, yr,                                 &
# endif
     &                          Pair )
!***********************************************************************
!
      USE mod_param
# ifdef RISSAGA
      USE mod_grid
      USE mod_scalars
# endif
!
#if defined EW_PERIODIC || defined NS_PERIODIC
      USE exchange_2d_mod, ONLY : exchange_r2d_tile
#endif
#ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange2d
#endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model, Iend, Istr, Jend, Jstr
      integer, intent(in) :: LBi, UBi, LBj, UBj
!
#ifdef ASSUMED_SHAPE
      real(r8), intent(out) :: Pair(LBi:,LBj:)
# ifdef RISSAGA
      real(r8), intent(in) :: xr(LBi:,LBj:)
      real(r8), intent(in) :: yr(LBi:,LBj:)
# endif
#else
      real(r8), intent(out) :: Pair(LBi:UBi,LBj:UBj)
# ifdef RISSAGA
      real(r8), intent(in) :: xr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: yr(LBi:UBi,LBj:UBj)
# endif
#endif
!
!  Local variable declarations.
!
#ifdef DISTRIBUTE
# ifdef EW_PERIODIC
      logical :: EWperiodic=.TRUE.
# else
      logical :: EWperiodic=.FALSE.
# endif
# ifdef NS_PERIODIC
      logical :: NSperiodic=.TRUE.
# else
      logical :: NSperiodic=.FALSE.
# endif
#endif
      integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
      integer :: i, j

#ifdef RISSAGA
	real(r8) :: dP, dir, spd, Rfront, pha
#endif

#include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Set analytical surface air pressure (mb).
!  (1 mb = 100 Pa = 1 hPa,  1 bar = 1.0e+5 N/m2 = 1.0e+5 dynes/cm2).
!-----------------------------------------------------------------------
!

#if defined RISSAGA
!
# ifdef ATM_PRESS
!
!  Propagating air pressure disturbance
!  speed of pressure front is           spd (m/s)
!  air pressure anomaly across front is dP (millibars)
!  length scale of frontal width is     Rfront (meters) 
!
!  The following are the first set of parameters tested
!     dP  = 7.0_r8 
!     dir = 20.0_r8*pi/180.0_r8
!     spd = 30.0_r8
!     Rfront = 18.0E+03_r8

!  User defined values set in ocean_rissaga.in
!  This allows us to run multiple parameter choices without having
!  to change this file and recompile roms
!
      dP     = user(1)
      dir    = user(2)*pi/180.0_r8  ! note dir was input in degrees
      spd    = user(3)
      Rfront = user(4)

      DO j=JstrR,JendR
        DO i=IstrR,IendR

!  Propagating phase argument is pha (direction dir, speed spd)
!  xr,yr are coordinates in meters from Mallorca grid file
!  Specify a sensible origin for xr,yr so that the front enters
!  the model domain soon after initialization.
!  Here we use x0 = 176.2e3 and y0 = 4632e3 which is 20 km southwest
!  of the southwest corner of the domain.
!
!  With this phase argument, we can use any Function(pha)
!  that describes the form of the air pressure perturbation as 
!  a function of distance (in meters) at a fixed time. This shape 
!  perturbation will then translate through the domain over time
!
!  Here the shape function is a tanh with length scale Lfront
!
           pha = (xr(i,j)-176.2E+03_r8)*COS(dir) +                     &
     &          (yr(i,j)-4632.0E+03_r8)*SIN(dir) -                     &
     &          spd * time(ng)
           pair(i,j) = 1020.0_r8 - dP * 0.5*(1.0_r8+TANH(pha/Rfront))
        END DO
      END DO
# else
      ana_pair.h: should not be here if not using atm_press option
# endif
#else
      ana_pair.h: no values provided for Pair.
#endif

#if defined EW_PERIODIC || defined NS_PERIODIC
      CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend,               &
     &                        LBi, UBi, LBj, UBj,                       &
     &                        Pair)
#endif
#ifdef DISTRIBUTE
      CALL mp_exchange2d (ng, model, 1, Istr, Iend, Jstr, Jend,         &
     &                    LBi, UBi, LBj, UBj,                           &
     &                    NghostPoints, EWperiodic, NSperiodic,         &
     &                    Pair)
#endif
      RETURN
      END SUBROUTINE ana_pair_tile