text:
MODULE dcmip_initial_conditions_test_1_2_3

  !=======================================================================
  !
  !  Functions for setting up initial conditions for the dynamical core tests:
  !
  !  11 - Deformational Advection Test
  !  12 - Hadley Cell Advection Test
  !  13 - Orography Advection Test
  !  20 - Impact of orography on a steady-state at rest
  !  21 and 22 - Non-Hydrostatic Mountain Waves Over A Schaer-Type Mountain without and with vertical wind shear
  !  31 - Non-Hydrostatic Gravity Waves
  !
  !  Given a point specified by:
  !   lon longitude (radians)
  !   lat latitude (radians)
  !   p/z pressure/height
  !  the functions will return:
  ! u zonal wind (m s^-1)
  ! v meridional wind (m s^-1)
  ! w vertical velocity (m s^-1)
  ! t temperature (K)
  ! phis  surface geopotential (m^2 s^-2)
  ! ps  surface pressure (Pa)
  ! rho density (kj m^-3)
  ! q specific humidity (kg/kg)
  ! qi  tracers (kg/kg)
  !     p       pressure if height based (Pa)
  !
  !
  !  Authors: James Kent, Paul Ullrich, Christiane Jablonowski
  !   (University of Michigan, dcmip@ucar.edu)
  !          version 4
  !          July/8/2012
  !
  ! Change log: (v3, June/8/2012, v4 July/8/2012, v5 July/20/2012)
  !
  ! v2: bug fixes in the tracer initialization for height-based models
  ! v3: test 3-1: the density is now initialized with the unperturbed background temperature (not the perturbed temperature)
  ! v3: constants converted to double precision
  ! v4: modified tracers in test 1-1, now with cutoff altitudes. Outside of the vertical domain all tracers are set to 0
  ! v4: modified cos-term in vertical velocity (now cos(2 pi t/tau)) in test 1-1, now completing one full up and down cycle
  ! v4: added subroutine test1_advection_orography for test 1-3
  ! v4: added subroutine test2_steady_state_mountain for test 2-0
  ! v4: modified parameter list for tests 2-1 and 2-2 (routine test2_schaer_mountain): addition of parameters hybrid_eta, hyam, hybm
  !     if the logical flag hybrid_eta is true then the pressure in pressure-based model with hybrid sigma-p (eta) coordinates is
  !     computed internally. In that case the hybrid coefficients hyam and hybm need to be supplied via the parameter list,
  !     otherwise they are not used.
  ! v5: Change in test 11 - change to u and w, cutoff altitudes (introduced in v4) are removed again
  ! v5: Change in test 12 - velocities multiplies by rho0/rho, different w0 and vertical location of the initial tracer
  !
  !
  !=======================================================================
  use mod_adm, only: &
     ADM_LOG_FID

  IMPLICIT NONE

  public :: test1_advection_deformation
  public :: test1_advection_hadley
  public :: test1_advection_orography
  public :: test2_steady_state_mountain
  public :: test2_schaer_mountain
  public :: test3_gravity_wave

!-----------------------------------------------------------------------
!     Physical Parameters
!-----------------------------------------------------------------------

  real(8), parameter :: a = 6371220.0d0,  & ! Earth's Radius (m)
        Rd  = 287.0d0,  & ! Ideal gas const dry air (J kg^-1 K^1)
        g = 9.80616d0,  & ! Gravity (m s^2)
        cp  = 1004.5d0, & ! Specific heat capacity (J kg^-1 K^1)
        pi  = 3.141592653589793238D0       ! pi

!-----------------------------------------------------------------------
!     Additional constants
!-----------------------------------------------------------------------

  real(8), parameter :: p0  = 100000.d0   ! reference pressure (Pa)

CONTAINS

!==========================================================================================
! TEST CASE 11 - PURE ADVECTION - 3D DEFORMATIONAL FLOW
!==========================================================================================

! The 3D deformational flow test is based on the deformational flow test of Nair and Lauritzen (JCP 2010),
! with a prescribed vertical wind velocity which makes the test truly 3D. An unscaled planet (with scale parameter
! X = 1) is selected.

SUBROUTINE test1_advection_deformation (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q,q1,q2,q3,q4)

IMPLICIT NONE
!-----------------------------------------------------------------------
!     input/output params parameters at given location
!-----------------------------------------------------------------------

  real(8), intent(in)  :: lon, &    ! Longitude (radians)
        lat, &    ! Latitude (radians)
        z   ! Height (m)

  real(8), intent(inout) :: p   ! Pressure  (Pa)

  integer,  intent(in) :: zcoords   ! 0 or 1 see below

  real(8), intent(out) :: u, &    ! Zonal wind (m s^-1)
        v, &    ! Meridional wind (m s^-1)
        w, &    ! Vertical Velocity (m s^-1)
        t, &    ! Temperature (K)
        phis, &   ! Surface Geopotential (m^2 s^-2)
        ps, &     ! Surface Pressure (Pa)
        rho, &    ! density (kg m^-3)
        q, &    ! Specific Humidity (kg/kg)
        q1, &     ! Tracer q1 (kg/kg)
        q2, &     ! Tracer q2 (kg/kg)
        q3, &     ! Tracer q3 (kg/kg)
        q4    ! Tracer q4 (kg/kg)

  ! if zcoords = 1, then we use z and output p
  ! if zcoords = 0, then we use p

!-----------------------------------------------------------------------
!     test case parameters
!-----------------------------------------------------------------------
  real(8), parameter ::   tau     = 12.d0 * 86400.d0, & ! period of motion 12 days
            u0      = (2.d0*pi*a)/tau,  & ! 2 pi a / 12 days
            k0  = (10.d0*a)/tau,  & ! Velocity Magnitude
            omega0  = (23000.d0*pi)/tau,  & ! Velocity Magnitude
                              T0      = 300.d0,   & ! temperature
                              H       = Rd * T0 / g,    & ! scale height
                              RR      = 1.d0/2.d0,    & ! horizontal half width divided by 'a'
                              ZZ      = 1000.d0,    & ! vertical half width
                              z0      = 5000.d0,    & ! center point in z
                              lambda0 = 5.d0*pi/6.d0,   & ! center point in longitudes
                              lambda1 = 7.d0*pi/6.d0,   & ! center point in longitudes
                              phi0    = 0.d0,     & ! center point in latitudes
                              phi1    = 0.d0

      real(8) :: height             ! The height of the model levels
      real(8) :: ptop             ! Model top in p
      real(8) :: sin_tmp, cos_tmp, sin_tmp2, cos_tmp2     ! Calculate great circle distances
      real(8) :: d1, d2, r, r2, d3, d4          ! For tracer calculations
      real(8) :: s, bs              ! Shape function, and parameter
      real(8) :: lonp             ! Translational longitude, depends on time
      real(8) :: ud             ! Divergent part of u
      real(8) :: time             ! Initially set to zero seconds, needs
                  ! to be modified when used in dycore

!-----------------------------------------------------------------------
!    HEIGHT AND PRESSURE
!-----------------------------------------------------------------------

  ! Height and pressure are aligned (p = p0 exp(-z/H))

  if (zcoords .eq. 1) then

    height = z
    p = p0 * exp(-z/H)

  else

    height = H * log(p0/p)

  endif

  ! Model top in p

  ptop    = p0*exp(-12000.d0/H)

!-----------------------------------------------------------------------
!    THE VELOCITIES ARE TIME DEPENDENT AND THEREFORE MUST BE UPDATED
!    IN THE DYNAMICAL CORE
!-----------------------------------------------------------------------

  ! These are initial conditions hence time = 0

  time = 0.d0

  ! Translational longitude = longitude when time = 0

  lonp = lon - 2.d0*pi*time/tau

  ! Shape function
!********
! change in version 5: shape function
!********
  bs = 0.2
  s = 1.0 + exp( (ptop-p0)/(bs*ptop) ) - exp( (p-p0)/(bs*ptop)) - exp( (ptop-p)/(bs*ptop))

  ! Zonal Velocity
!********
! change in version 5: ud
!********

  ud = (omega0*a)/(bs*ptop) * cos(lonp) * (cos(lat)**2.0) * cos(2.0*pi*time/tau) * &
    ( - exp( (p-p0)/(bs*ptop)) + exp( (ptop-p)/(bs*ptop))  )

  u = k0*sin(lonp)*sin(lonp)*sin(2.d0*lat)*cos(pi*time/tau) + u0*cos(lat) + ud

  ! Meridional Velocity

  v = k0*sin(2.d0*lonp)*cos(lat)*cos(pi*time/tau)

  ! Vertical Velocity - can be changed to vertical pressure velocity by
  ! omega = -(g*p)/(Rd*T0)*w
        !
!********
! change in version 4: cos(2.0*pi*time/tau) is now used instead of cos(pi*time/tau)
!********
!********
! change in version 5: shape function in w
!********
  w = -((Rd*T0)/(g*p))*omega0*sin(lonp)*cos(lat)*cos(2.0*pi*time/tau)*s

!-----------------------------------------------------------------------
!    TEMPERATURE IS CONSTANT 300 K
!-----------------------------------------------------------------------

  t = T0

!-----------------------------------------------------------------------
!    PHIS (surface geopotential)
!-----------------------------------------------------------------------

  phis = 0.d0

!-----------------------------------------------------------------------
!    PS (surface pressure)
!-----------------------------------------------------------------------

  ps = p0

!-----------------------------------------------------------------------
!    RHO (density)
!-----------------------------------------------------------------------

  rho = p/(Rd*t)

!-----------------------------------------------------------------------
!     initialize Q, set to zero
!-----------------------------------------------------------------------

  q = 0.d0

!-----------------------------------------------------------------------
!     initialize tracers
!-----------------------------------------------------------------------

  ! Tracer 1 - Cosine Bells

    ! To calculate great circle distance

  sin_tmp = sin(lat) * sin(phi0)
  cos_tmp = cos(lat) * cos(phi0)
  sin_tmp2 = sin(lat) * sin(phi1)
  cos_tmp2 = cos(lat) * cos(phi1)

    ! great circle distance without 'a'

  r  = ACOS (sin_tmp + cos_tmp*cos(lon-lambda0))
  r2  = ACOS (sin_tmp2 + cos_tmp2*cos(lon-lambda1))
  d1 = min( 1.d0, (r/RR)**2 + ((height-z0)/ZZ)**2 )
  d2 = min( 1.d0, (r2/RR)**2 + ((height-z0)/ZZ)**2 )

  q1 = 0.5d0 * (1.d0 + cos(pi*d1)) + 0.5d0 * (1.d0 + cos(pi*d2))

  ! Tracer 2 - Correlated Cosine Bells

  q2 = 0.9d0 - 0.8d0*q1**2

  ! Tracer 3 - Slotted Ellipse

    ! Make the ellipse

  if (d1 .le. RR) then
    q3 = 1.d0
  elseif (d2 .le. RR) then
    q3 = 1.d0
  else
    q3 = 0.1d0
  endif

    ! Put in the slot

  if (height .gt. z0 .and. abs(lat) .lt. 0.125d0) then

    q3 = 0.1d0

  endif

  ! Tracer 4: q4 is chosen so that, in combination with the other three tracer
        !           fields with weight (3/10), the sum is equal to one

  q4 = 1.d0 - 0.3d0*(q1+q2+q3)

!************
! change in version 4: added cutoff altitudes, tracers are set to zero outside this region
! prevents tracers from being trapped near the bottom and top of the domain
!************
        ! use a cutoff altitude for
  ! tracer 2 3 and 4
  ! Set them to zero outside `buffer zone'
!************
! change in version 5: change from v4 reversed, no cutoff altitudes due to continuity equation being satisfied
!       commented out below
!************

  !if (height .gt. (z0+1.25*ZZ) .or. height .lt. (z0-1.25*ZZ)) then

  ! q2 = 0.0
  ! q3 = 0.0
  ! q4 = 0.0

  !endif

END SUBROUTINE test1_advection_deformation

!==========================================================================================
! TEST CASE 12 - PURE ADVECTION - 3D HADLEY-LIKE FLOW
!==========================================================================================

SUBROUTINE test1_advection_hadley (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q,q1)

IMPLICIT NONE
!-----------------------------------------------------------------------
!     input/output params parameters at given location
!-----------------------------------------------------------------------

  real(8), intent(in)  :: lon, &    ! Longitude (radians)
        lat, &    ! Latitude (radians)
        z   ! Height (m)

  real(8), intent(inout) :: p   ! Pressure  (Pa)

  integer,  intent(in) :: zcoords   ! 0 or 1 see below

  real(8), intent(out) :: u, &    ! Zonal wind (m s^-1)
        v, &    ! Meridional wind (m s^-1)
        w, &    ! Vertical Velocity (m s^-1)
        t, &    ! Temperature (K)
        phis, &   ! Surface Geopotential (m^2 s^-2)
        ps, &     ! Surface Pressure (Pa)
        rho, &    ! density (kg m^-3)
        q, &    ! Specific Humidity (kg/kg)
        q1    ! Tracer q1 (kg/kg)

  ! if zcoords = 1, then we use z and output p
  ! if zcoords = 0, then we use p

!-----------------------------------------------------------------------
!     test case parameters
!-----------------------------------------------------------------------
  real(8), parameter ::   tau     = 1.d0 * 86400.d0,  & ! period of motion 1 day (in s)
            u0      = 40.d0,    & ! Zonal velocity magnitude (m/s)
            w0  = 0.15d0,   & ! Vertical velocity magnitude (m/s), changed in v5
                              T0      = 300.d0,   & ! temperature (K)
                              H       = Rd * T0 / g,    & ! scale height
                              K       = 5.d0,     & ! number of Hadley-like cells
                              z1      = 2000.d0,    & ! position of lower tracer bound (m), changed in v5
                              z2      = 5000.d0,    & ! position of upper tracer bound (m), changed in v5
                              z0      = 0.5d0*(z1+z2),  & ! midpoint (m)
                              ztop    = 12000.d0      ! model top (m)

      real(8) :: rho0             ! reference density at z=0 m
      real(8) :: height             ! Model level heights
      real(8) :: time             ! Initially set to zero seconds, needs
                  ! to be modified when used in dycore

!-----------------------------------------------------------------------
!    HEIGHT AND PRESSURE
!-----------------------------------------------------------------------

  ! Height and pressure are aligned (p = p0 exp(-z/H))

  if (zcoords .eq. 1) then

    height = z
    p = p0 * exp(-z/H)

  else

    height = H * log(p0/p)

  endif

!-----------------------------------------------------------------------
!    TEMPERATURE IS CONSTANT 300 K
!-----------------------------------------------------------------------

  t = T0

!-----------------------------------------------------------------------
!    PHIS (surface geopotential)
!-----------------------------------------------------------------------

  phis = 0.d0

!-----------------------------------------------------------------------
!    PS (surface pressure)
!-----------------------------------------------------------------------

  ps = p0

!-----------------------------------------------------------------------
!    RHO (density)
!-----------------------------------------------------------------------

  rho = p/(Rd*t)
  rho0 = p0/(Rd*t)

!-----------------------------------------------------------------------
!    THE VELOCITIES ARE TIME DEPENDENT AND THEREFORE MUST BE UPDATED
!    IN THE DYNAMICAL CORE
!-----------------------------------------------------------------------

  time = 0.d0

  ! Zonal Velocity

  u = u0*cos(lat)

  ! Meridional Velocity

!************
! change in version 5: multiply v and w by rho0/rho
!************

  v = -(rho0/rho) * (a*w0*pi)/(K*ztop) *cos(lat)*sin(K*lat)*cos(pi*height/ztop)*cos(pi*time/tau)

  ! Vertical Velocity - can be changed to vertical pressure velocity by
  ! omega = -g*rho*w

  w = (rho0/rho) *(w0/K)*(-2.d0*sin(K*lat)*sin(lat) + K*cos(lat)*cos(K*lat)) &
    *sin(pi*height/ztop)*cos(pi*time/tau)

!-----------------------------------------------------------------------
!     initialize Q, set to zero
!-----------------------------------------------------------------------

  q = 0.d0

!-----------------------------------------------------------------------
!     initialize tracers
!-----------------------------------------------------------------------

  ! Tracer 1 - Layer

  if (height .lt. z2 .and. height .gt. z1) then

    q1 = 0.5d0 * (1.d0 + cos( 2.d0*pi*(height-z0)/(z2-z1) ) )

  else

    q1 = 0.d0

  endif

END SUBROUTINE test1_advection_hadley

!==========================================================================================
! TEST CASE 13 - HORIZONTAL ADVECTION OF THIN CLOUD-LIKE TRACERS IN THE PRESENCE OF OROGRAPHY
!==========================================================================================

SUBROUTINE test1_advection_orography (lon,lat,p,z,zcoords,cfv,hybrid_eta,hyam,hybm,gc,u,v,w,t,phis,ps,rho,q,q1,q2,q3,q4)

IMPLICIT NONE
!-----------------------------------------------------------------------
!     input/output params parameters at given location
!-----------------------------------------------------------------------

  real(8), intent(in)  :: lon, &    ! Longitude (radians)
        lat, &    ! Latitude (radians)
        z, &    ! Height (m)
        hyam, &   ! A coefficient for hybrid-eta coordinate, at model level midpoint
        hybm, &   ! B coefficient for hybrid-eta coordinate, at model level midpoint
        gc    ! bar{z} for Gal-Chen coordinate

  logical, intent(in)  :: hybrid_eta      ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used
                                                ! if set to .true., then the pressure will be computed via the
                                                !    hybrid coefficients hyam and hybm, they need to be initialized
                                                ! if set to .false. (for pressure-based models): the pressure is already pre-computed
                                                !    and is an input value for this routine
                                                ! for height-based models: pressure will always be computed based on the height and
                                                !    hybrid_eta is not used

  ! Note that we only use hyam and hybm for the hybrid-eta coordiantes, and we only use
  ! gc for the Gal-Chen coordinates. If not required then they become dummy variables

  real(8), intent(inout) :: p   ! Pressure  (Pa)

  integer,  intent(in) :: zcoords   ! 0 or 1 see below
  integer,  intent(in) :: cfv     ! 0, 1 or 2 see below

  real(8), intent(out) :: u, &    ! Zonal wind (m s^-1)
        v, &    ! Meridional wind (m s^-1)
        w, &    ! Vertical Velocity (m s^-1)
        t, &    ! Temperature (K)
        phis, &   ! Surface Geopotential (m^2 s^-2)
        ps, &     ! Surface Pressure (Pa)
        rho, &    ! density (kg m^-3)
        q, &    ! Specific Humidity (kg/kg)
        q1, &     ! Tracer q1 (kg/kg)
        q2, &     ! Tracer q2 (kg/kg)
        q3, &     ! Tracer q3 (kg/kg)
        q4    ! Tracer q4 (kg/kg)

  ! if zcoords = 1, then we use z and output p
  ! if zcoords = 0, then we use p

  ! if cfv = 0 we assume that our horizontal velocities are not coordinate following
  ! if cfv = 1 then our velocities follow hybrid eta coordinates and we need to specify w
  ! if cfv = 2 then our velocities follow Gal-Chen coordinates and we need to specify w

  ! In hybrid-eta coords: p = hyam p0 + hybm ps
  ! In Gal-Chen coords: z = zs + (gc/ztop)*(ztop - zs)

        ! if other orography-following coordinates are used, the w wind needs to be newly derived for them

!-----------------------------------------------------------------------
!     test case parameters
!-----------------------------------------------------------------------
  real(8), parameter ::   tau     = 12.d0 * 86400.d0, & ! period of motion 12 days (s)
            u0      = 2.d0*pi*a/tau,  & ! Velocity Magnitude (m/s)
                              T0      = 300.d0,   & ! temperature (K)
                              H       = Rd * T0 / g,    & ! scale height (m)
                              alpha   = pi/6.d0,    & ! rotation angle (radians), 30 degrees
                              lambdam = 3.d0*pi/2.d0,   & ! mountain longitude center point (radians)
                              phim    = 0.d0,     & ! mountain latitude center point (radians)
                              h0      = 2000.d0,    & ! peak height of the mountain range (m)
                              Rm      = 3.d0*pi/4.d0,   & ! mountain radius (radians)
                              zetam   = pi/16.d0,   & ! mountain oscillation half-width (radians)
                              lambdap = pi/2.d0,    & ! cloud-like tracer longitude center point (radians)
                              phip    = 0.d0,       & ! cloud-like tracer latitude center point (radians)
                              Rp      = pi/4.d0,    & ! cloud-like tracer radius (radians)
                              zp1     = 3050.d0,    & ! midpoint of first (lowermost) tracer (m)
                              zp2     = 5050.d0,    & ! midpoint of second tracer (m)
                              zp3     = 8200.d0,    & ! midpoint of third (topmost) tracer (m)
                              dzp1    = 1000.d0,    & ! thickness of first (lowermost) tracer (m)
                              dzp2    = 1000.d0,    & ! thickness of second tracer (m)
                              dzp3    = 400.d0,   & ! thickness of third (topmost) tracer (m)
                              ztop    = 12000.d0      ! model top (m)

      real(8) :: height             ! Model level heights (m)
      real(8) :: r              ! Great circle distance (radians)
      real(8) :: rz             ! height differences
      real(8) :: zs             ! Surface elevation (m)

!-----------------------------------------------------------------------
!    PHIS (surface geopotential)
!-----------------------------------------------------------------------

  r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) )

  if (r .lt. Rm) then

    zs = (h0/2.d0)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)**2.d0

  else

    zs = 0.d0

  endif

  phis = g*zs

!-----------------------------------------------------------------------
!    PS (surface pressure)
!-----------------------------------------------------------------------

  ps = p0 * exp(-zs/H)

!-----------------------------------------------------------------------
!    HEIGHT AND PRESSURE
!-----------------------------------------------------------------------

  ! Height and pressure are aligned (p = p0 exp(-z/H))

  if (zcoords .eq. 1) then

    height = z
    p = p0 * exp(-z/H)

  else

                if (hybrid_eta) p = hyam*p0 + hybm*ps   ! compute the pressure based on the surface pressure and hybrid coefficients
    height = H * log(p0/p)

  endif

!-----------------------------------------------------------------------
!    THE VELOCITIES ARE TIME INDEPENDENT
!-----------------------------------------------------------------------

  ! Zonal Velocity

  u = u0*(cos(lat)*cos(alpha)+sin(lat)*cos(lon)*sin(alpha))

  ! Meridional Velocity

  v = -u0*(sin(lon)*sin(alpha))

!-----------------------------------------------------------------------
!    TEMPERATURE IS CONSTANT 300 K
!-----------------------------------------------------------------------

  t = T0

!-----------------------------------------------------------------------
!    RHO (density)
!-----------------------------------------------------------------------

  rho = p/(Rd*t)

!-----------------------------------------------------------------------
!     initialize Q, set to zero
!-----------------------------------------------------------------------

  q = 0.d0

!-----------------------------------------------------------------------
!    VERTICAL VELOCITY IS TIME INDEPENDENT
!-----------------------------------------------------------------------

  ! Vertical Velocity - can be changed to vertical pressure velocity by
  ! omega = -(g*p)/(Rd*T0)*w

  ! NOTE that if orography-following coordinates are used then the vertical
  ! velocity needs to be translated into the new coordinate system due to
  ! the variation of the height along coordinate surfaces
  ! See section 1.3 and the appendix of the test case document

  if (cfv .eq. 0) then

    ! if the horizontal velocities do not follow the vertical coordinate

    w = 0.d0

  elseif (cfv .eq. 1) then

    ! if the horizontal velocities follow hybrid eta coordinates then
    ! the perceived vertical velocity is

    call test1_advection_orograph_hybrid_eta_velocity(w)

  elseif (cfv .eq. 2) then

    ! if the horizontal velocities follow Gal Chen coordinates then
    ! the perceived vertical velocity is

    call test1_advection_orograph_Gal_Chen_velocity(w)

!        else
!               compute your own vertical velocity if other orography-following
!               vertical coordinate is used
!
  endif

!-----------------------------------------------------------------------
!     initialize tracers
!-----------------------------------------------------------------------

  ! Tracer 1 - Cloud Layer

  r = acos( sin(phip)*sin(lat) + cos(phip)*cos(lat)*cos(lon - lambdap) )

  rz = abs(height - zp1)

  if (rz .lt. 0.5d0*dzp1 .and. r .lt. Rp) then

    q1 = 0.25d0*(1.d0+cos(2.d0*pi*rz/dzp1))*(1.d0+cos(pi*r/Rp))

  else

    q1 = 0.d0

  endif

  rz = abs(height - zp2)

  if (rz .lt. 0.5d0*dzp2 .and. r .lt. Rp) then

    q2 = 0.25d0*(1.d0+cos(2.d0*pi*rz/dzp2))*(1.d0+cos(pi*r/Rp))

  else

    q2 = 0.d0

  endif

  rz = abs(height - zp3)

  if (rz .lt. 0.5d0*dzp3 .and. r .lt. Rp) then

    q3 = 1.d0

  else

    q3 = 0.d0

  endif

  q4 = q1 + q2 + q3

  CONTAINS

  !-----------------------------------------------------------------------
  !    SUBROUTINE TO CALCULATE THE PERCEIVED VERTICAL VELOCITY
  !       UNDER HYBRID-ETA COORDINATES
  !-----------------------------------------------------------------------

  SUBROUTINE test1_advection_orograph_hybrid_eta_velocity(w)
  IMPLICIT NONE
  real(8), intent(out) :: w

  real(8) ::  press, &    ! hyam *p0 + hybm *ps
      r, &      ! Great Circle Distance
      dzsdx, &    ! Part of surface height derivative
      dzsdlambda, &     ! Derivative of zs w.r.t lambda
      dzsdphi, &    ! Derivative of zs w.r.t phi
      dzdlambda, &    ! Derivative of z w.r.t lambda
      dzdphi, &     ! Derivative of z w.r.t phi
      dpsdlambda, &     ! Derivative of ps w.r.t lambda
      dpsdphi     ! Derivative of ps w.r.t phi

    ! Calculate pressure and great circle distance to mountain center

    press = hyam*p0 + hybm*ps

    r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) )

    ! Derivatives of surface height

    if (r .lt. Rm) then
      dzsdx = -h0*pi/(2.d0*Rm)*sin(pi*r/Rm)*cos(pi*r/zetam)**2 - &
        (h0*pi/zetam)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)*sin(pi*r/zetam)
    else
      dzsdx = 0.d0
    endif

    ! Prevent division by zero

    if (1.d0-cos(r)**2 .gt. 0.d0) then
      dzsdlambda = dzsdx * (cos(phim)*cos(lat)*sin(lon-lambdam)) &
          /sqrt(1.d0-cos(r)**2)
      dzsdphi    = dzsdx * (-sin(phim)*cos(lat) + cos(phim)*sin(lat)*cos(lon-lambdam)) &
          /sqrt(1.d0-cos(r)**2)
    else
      dzsdlambda = 0.d0
      dzsdphi    = 0.d0
    endif

    ! Derivatives of surface pressure

    dpsdlambda = -(g*p0/(Rd*T0))*exp(-g*zs/(Rd*T0))*dzsdlambda
    dpsdphi    = -(g*p0/(Rd*T0))*exp(-g*zs/(Rd*T0))*dzsdphi

    ! Derivatives of coordinate height

    dzdlambda = -(Rd*T0/(g*press))*hybm*dpsdlambda
    dzdphi    = -(Rd*T0/(g*press))*hybm*dpsdphi

    ! Prevent division by zero

    if (abs(lat) .lt. pi/2.d0) then
      w = - (u/(a*cos(lat)))*dzdlambda - (v/a)*dzdphi
    else
      w = 0.d0
    endif

  END SUBROUTINE test1_advection_orograph_hybrid_eta_velocity

  !-----------------------------------------------------------------------
  !    SUBROUTINE TO CALCULATE THE PERCEIVED VERTICAL VELOCITY
  !       UNDER GAL-CHEN COORDINATES
  !-----------------------------------------------------------------------

  SUBROUTINE test1_advection_orograph_Gal_Chen_velocity(w)
  IMPLICIT NONE
  real(8), intent(out) :: w

  real(8) ::  r, &      ! Great Circle Distance
      dzsdx, &    ! Part of surface height derivative
      dzsdlambda, &     ! Derivative of zs w.r.t lambda
      dzsdphi, &    ! Derivative of zs w.r.t phi
      dzdlambda, &    ! Derivative of z w.r.t lambda
      dzdphi      ! Derivative of z w.r.t phi

    ! Calculate great circle distance to mountain center

    r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) )

    ! Derivatives of surface height

    if (r .lt. Rm) then
      dzsdx = -h0*pi/(2.d0*Rm)*sin(pi*r/Rm)*cos(pi*r/zetam)**2 - &
        (h0*pi/zetam)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)*sin(pi*r/zetam)
    else
      dzsdx = 0.d0
    endif

    ! Prevent division by zero

    if (1.d0-cos(r)**2 .gt. 0.d0) then
      dzsdlambda = dzsdx * (cos(phim)*cos(lat)*sin(lon-lambdam)) &
          /sqrt(1.d0-cos(r)**2)
      dzsdphi    = dzsdx * (-sin(phim)*cos(lat) + cos(phim)*sin(lat)*cos(lon-lambdam)) &
          /sqrt(1.d0-cos(r)**2)
    else
      dzsdlambda = 0.d0
      dzsdphi    = 0.d0
    endif

    ! Derivatives of coordinate height

    dzdlambda = (1.d0-gc/ztop)*dzsdlambda
    dzdphi    = (1.d0-gc/ztop)*dzsdphi

    ! Prevent division by zero

    if (abs(lat) .lt. pi/2.d0) then
      w = - (u/(a*cos(lat)))*dzdlambda - (v/a)*dzdphi
    else
      w = 0.d0
    endif

  END SUBROUTINE test1_advection_orograph_Gal_Chen_velocity

END SUBROUTINE test1_advection_orography

!==========================================================================================
! TEST CASE 2X - IMPACT OF OROGRAPHY ON A NON-ROTATING PLANET
!==========================================================================================
! The tests in section 2-x examine the impact of 3D Schaer-like circular mountain profiles on an
! atmosphere at rest (2-0), and on flow fields with wind shear (2-1) and without vertical wind shear (2-2).
! A non-rotating planet is used for all configurations. Test 2-0 is conducted on an unscaled regular-size
! planet and primarily examines the accuracy of the pressure gradient calculation in a steady-state
! hydrostatically-balanced atmosphere at rest. This test is especially appealing for models with
! orography-following vertical coordinates. It increases the complexity of test 1-3, that investigated
! the impact of the same Schaer-type orographic profile on the accuracy of purely-horizontal passive
! tracer advection.
!
! Tests 2-1 and 2-2 increase the complexity even further since non-zero flow fields are now prescribed
! with and without vertical wind shear. In order to trigger non-hydrostatic responses the two tests are
! conducted on a reduced-size planet with reduction factor $X=500$ which makes the horizontal and
! vertical grid spacing comparable. This test clearly discriminates between non-hydrostatic and hydrostatic
! models since the expected response is in the non-hydrostatic regime. Therefore, the flow response is
! captured differently by hydrostatic models.

!=========================================================================
! Test 2-0:  Steady-State Atmosphere at Rest in the Presence of Orography
!=========================================================================
SUBROUTINE test2_steady_state_mountain (lon,lat,p,z,zcoords,hybrid_eta,hyam,hybm,u,v,w,t,phis,ps,rho,q)

IMPLICIT NONE
!-----------------------------------------------------------------------
!     input/output params parameters at given location
!-----------------------------------------------------------------------

  real(8), intent(in)  :: lon, &    ! Longitude (radians)
        lat, &    ! Latitude (radians)
        z, &    ! Height (m)
        hyam, &   ! A coefficient for hybrid-eta coordinate, at model level midpoint
        hybm    ! B coefficient for hybrid-eta coordinate, at model level midpoint

  logical, intent(in)  :: hybrid_eta      ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used
                                                ! if set to .true., then the pressure will be computed via the
                                                !    hybrid coefficients hyam and hybm, they need to be initialized
                                                ! if set to .false. (for pressure-based models): the pressure is already pre-computed
                                                !    and is an input value for this routine
                                                ! for height-based models: pressure will always be computed based on the height and
                                                !    hybrid_eta is not used

  real(8), intent(inout) :: p   ! Pressure  (Pa)

  integer,  intent(in) :: zcoords   ! 0 or 1 see below

  real(8), intent(out) :: u, &    ! Zonal wind (m s^-1)
        v, &    ! Meridional wind (m s^-1)
        w, &    ! Vertical Velocity (m s^-1)
        t, &    ! Temperature (K)
        phis, &   ! Surface Geopotential (m^2 s^-2)
        ps, &     ! Surface Pressure (Pa)
        rho, &    ! density (kg m^-3)
        q     ! Specific Humidity (kg/kg)

  ! if zcoords = 1, then we use z and output p
  ! if zcoords = 0, then we compute or use p
        !
  ! In hybrid-eta coords: p = hyam p0 + hybm ps
        !
        ! The grid-point based initial data are computed in this routine.

!-----------------------------------------------------------------------
!     test case parameters
!-----------------------------------------------------------------------
  real(8), parameter ::   T0      = 300.d0,   & ! temperature (K)
                                gamma   = 0.0065d0,             &       ! temperature lapse rate (K/m)
                              lambdam = 3.d0*pi/2.d0,   & ! mountain longitude center point (radians)
                              phim    = 0.d0,     & ! mountain latitude center point (radians)
                              h0      = 2000.d0,    & ! peak height of the mountain range (m)
                              Rm      = 3.d0*pi/4.d0,   & ! mountain radius (radians)
                              zetam   = pi/16.d0,   & ! mountain oscillation half-width (radians)
                              ztop    = 12000.d0      ! model top (m)

      real(8) :: height             ! Model level heights (m)
      real(8) :: r              ! Great circle distance (radians)
      real(8) :: zs             ! Surface elevation (m)
      real(8) :: exponent                                               ! exponent: g/(Rd * gamma)
      real(8) :: exponent_rev                                           ! reversed exponent

!-----------------------------------------------------------------------
!    compute exponents
!-----------------------------------------------------------------------
     exponent     = g/(Rd*gamma)
     exponent_rev = 1.d0/exponent

!-----------------------------------------------------------------------
!    PHIS (surface geopotential)
!-----------------------------------------------------------------------

  r = acos( sin(phim)*sin(lat) + cos(phim)*cos(lat)*cos(lon - lambdam) )

  if (r .lt. Rm) then

    zs = (h0/2.d0)*(1.d0+cos(pi*r/Rm))*cos(pi*r/zetam)**2.d0   ! mountain height

  else

    zs = 0.d0

  endif

  phis = g*zs

!-----------------------------------------------------------------------
!    PS (surface pressure)
!-----------------------------------------------------------------------

  ps = p0 * (1.d0 - gamma/T0*zs)**exponent

!-----------------------------------------------------------------------
!    HEIGHT AND PRESSURE
!-----------------------------------------------------------------------

  ! Height and pressure are aligned (p = p0 * (1.d0 - gamma/T0*z)**exponent)

  if (zcoords .eq. 1) then

    height = z
    p = p0 * (1.d0 - gamma/T0*z)**exponent

  else

                if (hybrid_eta) p = hyam*p0 + hybm*ps              ! compute the pressure based on the surface pressure and hybrid coefficients
    height = T0/gamma * (1.d0 - (p/p0)**exponent_rev)  ! compute the height at this pressure

  endif

!-----------------------------------------------------------------------
!    THE VELOCITIES ARE ZERO (STATE AT REST)
!-----------------------------------------------------------------------

  ! Zonal Velocity

  u = 0.d0

  ! Meridional Velocity

  v = 0.d0

        ! Vertical Velocity

        w = 0.d0

!-----------------------------------------------------------------------
!    TEMPERATURE WITH CONSTANT LAPSE RATE
!-----------------------------------------------------------------------

  t = T0 - gamma*height

!-----------------------------------------------------------------------
!    RHO (density)
!-----------------------------------------------------------------------

  rho = p/(Rd*t)

!-----------------------------------------------------------------------
!     initialize Q, set to zero
!-----------------------------------------------------------------------

  q = 0.d0

END SUBROUTINE test2_steady_state_mountain

!=====================================================================================
! Tests 2-1 and 2-2:  Non-hydrostatic Mountain Waves over a Schaer-type Mountain
!=====================================================================================

SUBROUTINE test2_schaer_mountain (lon,lat,p,z,zcoords,hybrid_eta,hyam,hybm,shear,u,v,w,t,phis,ps,rho,q)

IMPLICIT NONE
!-----------------------------------------------------------------------
!     input/output params parameters at given location
!-----------------------------------------------------------------------

  real(8), intent(in)  :: lon, &    ! Longitude (radians)
        lat, &    ! Latitude (radians)
        z,   &    ! Height (m)
        hyam, &   ! A coefficient for hybrid-eta coordinate, at model level midpoint
        hybm    ! B coefficient for hybrid-eta coordinate, at model level midpoint

  logical, intent(in)  :: hybrid_eta      ! flag to indicate whether the hybrid sigma-p (eta) coordinate is used
                                                ! if set to .true., then the pressure will be computed via the
                                                !    hybrid coefficients hyam and hybm, they need to be initialized
                                                ! if set to .false. (for pressure-based models): the pressure is already pre-computed
                                                !    and is an input value for this routine
                                                ! for height-based models: pressure will always be computed based on the height and
                                                !    hybrid_eta is not used

  real(8), intent(inout) :: p   ! Pressure  (Pa)

  integer,  intent(in) :: zcoords, &  ! 0 or 1 see below
        shear   ! 0 or 1 see below

  real(8), intent(out) :: u, &    ! Zonal wind (m s^-1)
        v, &    ! Meridional wind (m s^-1)
        w, &    ! Vertical Velocity (m s^-1)
        t, &    ! Temperature (K)
        phis, &   ! Surface Geopotential (m^2 s^-2)
        ps, &     ! Surface Pressure (Pa)
        rho, &    ! density (kg m^-3)
        q     ! Specific Humidity (kg/kg)

  ! if zcoords = 1, then we use z and output p
  ! if zcoords = 0, then we either compute or use p

  ! if shear = 1, then we use shear flow
  ! if shear = 0, then we use constant u

!-----------------------------------------------------------------------
!     test case parameters
!-----------------------------------------------------------------------
  real(8), parameter ::   X       = 500.d0,   & ! Reduced Earth reduction factor
            Om      = 0.d0,     & ! Rotation Rate of Earth
                              as      = a/X,      & ! New Radius of small Earth
            ueq     = 20.d0,    & ! Reference Velocity
                              Teq     = 300.d0,   & ! Temperature at Equator
            Peq     = 100000.d0,    & ! Reference PS at Equator
                              ztop    = 30000.d0,   & ! Model Top
        lambdac = pi/4.d0,    &   ! Lon of Schar Mountain Center
        phic    = 0.d0,     &   ! Lat of Schar Mountain Center
        h0      = 250.d0,     &   ! Height of Mountain
        d       = 5000.d0,    &   ! Mountain Half-Width
        xi      = 4000.d0,    &   ! Mountain Wavelength
        cs      = 0.00025d0       ! Wind Shear (shear=1)

      real(8) :: height             ! Model level heights
      real(8) :: sin_tmp, cos_tmp         ! Calculation of great circle distance
      real(8) :: r              ! Great circle distance
      real(8) :: zs             ! Surface height
      real(8) :: c              ! Shear

!-----------------------------------------------------------------------
!    PHIS (surface geopotential)
!-----------------------------------------------------------------------

  sin_tmp = sin(lat) * sin(phic)
  cos_tmp = cos(lat) * cos(phic)

  ! great circle distance with 'a/X'

  r  = as * ACOS (sin_tmp + cos_tmp*cos(lon-lambdac))
  zs   = h0 * exp(-(r**2)/(d**2))*(cos(pi*r/xi)**2)
  phis = g*zs

!-----------------------------------------------------------------------
!    SHEAR FLOW OR CONSTANT FLOW
!-----------------------------------------------------------------------

  if (shear .eq. 1) then

    c = cs

  else

    c = 0.d0

  endif

!-----------------------------------------------------------------------
!    TEMPERATURE
!-----------------------------------------------------------------------

  t = Teq *(1.d0 - (c*ueq*ueq/(g))*(sin(lat)**2) )

!-----------------------------------------------------------------------
!    PS (surface pressure)
!-----------------------------------------------------------------------

  ps = peq*exp( -(ueq*ueq/(2.d0*Rd*Teq))*(sin(lat)**2) - phis/(Rd*t)    )

!-----------------------------------------------------------------------
!    HEIGHT AND PRESSURE
!-----------------------------------------------------------------------

  if (zcoords .eq. 1) then

    height = z
    p = peq*exp( -(ueq*ueq/(2.d0*Rd*Teq))*(sin(lat)**2) - g*height/(Rd*t)    )

  else
                if (hybrid_eta) p = hyam*p0 + hybm*ps ! compute the pressure based on the surface pressure and hybrid coefficients
    height = (Rd*t/(g))*log(peq/p) - (t*ueq*ueq/(2.d0*Teq*g))*(sin(lat)**2)

  endif

!-----------------------------------------------------------------------
!    THE VELOCITIES
!-----------------------------------------------------------------------

  ! Zonal Velocity

  u = ueq * cos(lat) * sqrt( (2.d0*Teq/(t))*c*height + t/(Teq) )

  ! Meridional Velocity

  v = 0.d0

  ! Vertical Velocity = Vertical Pressure Velocity = 0

  w = 0.d0

!-----------------------------------------------------------------------
!    RHO (density)
!-----------------------------------------------------------------------

  rho = p/(Rd*t)

!-----------------------------------------------------------------------
!     initialize Q, set to zero
!-----------------------------------------------------------------------

  q = 0.d0

END SUBROUTINE test2_schaer_mountain

!==========================================================================================
! TEST CASE 3 - GRAVITY WAVES
!==========================================================================================

! The non-hydrostatic gravity wave test examines the response of models to short time-scale wavemotion
! triggered by a localized perturbation. The formulation presented in this document is new,
! but is based on previous approaches by Skamarock et al. (JAS 1994), Tomita and Satoh (FDR 2004), and
! Jablonowski et al. (NCAR Tech Report 2008)

!==========
! Test 3-1
!==========
SUBROUTINE test3_gravity_wave (lon,lat,p,z,zcoords,u,v,w,t,phis,ps,rho,q)

IMPLICIT NONE
!-----------------------------------------------------------------------
!     input/output params parameters at given location
!-----------------------------------------------------------------------

  real(8), intent(in)  :: lon, &    ! Longitude (radians)
        lat, &    ! Latitude (radians)
        z   ! Height (m)

  real(8), intent(inout) :: p   ! Pressure  (Pa)

  integer,  intent(in) :: zcoords   ! 0 or 1 see below

  real(8), intent(out) :: u, &    ! Zonal wind (m s^-1)
        v, &    ! Meridional wind (m s^-1)
        w, &    ! Vertical Velocity (m s^-1)
        t, &    ! Temperature (K)
        phis, &   ! Surface Geopotential (m^2 s^-2)
        ps, &     ! Surface Pressure (Pa)
        rho, &    ! density (kg m^-3)
        q     ! Specific Humidity (kg/kg)

  ! if zcoords = 1, then we use z and output z
  ! if zcoords = 0, then we use p

!-----------------------------------------------------------------------
!     test case parameters
!-----------------------------------------------------------------------
  real(8), parameter ::   X       = 125.d0,   & ! Reduced Earth reduction factor
            Om      = 0.d0,     & ! Rotation Rate of Earth
                              as      = a/X,      & ! New Radius of small Earth
            u0      = 20.d0,    & ! Reference Velocity
                              Teq     = 300.d0,   & ! Temperature at Equator
            Peq     = 100000.d0,    & ! Reference PS at Equator
                              ztop    = 10000.d0,   & ! Model Top
        lambdac = 2.d0*pi/3.d0,   &   ! Lon of Pert Center
        d       = 5000.d0,    &   ! Width for Pert
        phic    = 0.d0,     &   ! Lat of Pert Center
        delta_theta = 1.d0,     &   ! Max Amplitude of Pert
        Lz      = 20000.d0,     &   ! Vertical Wavelength of Pert
        N       = 0.01d0,     &   ! Brunt-Vaisala frequency
        N2      = N*N,      & ! Brunt-Vaisala frequency Squared
        bigG    = (g*g)/(N2*cp)     ! Constant

      real(8) :: height             ! Model level height
      real(8) :: sin_tmp, cos_tmp         ! Calculation of great circle distance
      real(8) :: r, s             ! Shape of perturbation
      real(8) :: TS               ! Surface temperature
      real(8) :: t_mean, t_pert           ! Mean and pert parts of temperature
      real(8) :: theta_pert           ! Pot-temp perturbation

!-----------------------------------------------------------------------
!    THE VELOCITIES
!-----------------------------------------------------------------------

  ! Zonal Velocity

  u = u0 * cos(lat)

  ! Meridional Velocity

  v = 0.d0

  ! Vertical Velocity = Vertical Pressure Velocity = 0

  w = 0.d0

!-----------------------------------------------------------------------
!    PHIS (surface geopotential)
!-----------------------------------------------------------------------

  phis = 0.d0

!-----------------------------------------------------------------------
!    SURFACE TEMPERATURE
!-----------------------------------------------------------------------

  TS = bigG + (Teq-bigG)*exp( -(u0*N2/(4.d0*g*g))*(u0+2.d0*om*as)*(cos(2.d0*lat)-1.d0)    )

!-----------------------------------------------------------------------
!    PS (surface pressure)
!-----------------------------------------------------------------------

  ps = peq*exp( (u0/(4.0*bigG*Rd))*(u0+2.0*Om*as)*(cos(2.0*lat)-1.0)  ) &
    * (TS/Teq)**(cp/Rd)

!-----------------------------------------------------------------------
!    HEIGHT AND PRESSURE AND MEAN TEMPERATURE
!-----------------------------------------------------------------------

  if (zcoords .eq. 1) then

    height = z
    p = ps*( (bigG/TS)*exp(-N2*height/g)+1.d0 - (bigG/TS)  )**(cp/Rd)

  else

    height = (-g/N2)*log( (TS/bigG)*( (p/ps)**(Rd/cp) - 1.d0  ) + 1.d0 )

  endif

  t_mean = bigG*(1.d0 - exp(N2*height/g))+ TS*exp(N2*height/g)

!-----------------------------------------------------------------------
!    rho (density), unperturbed using the background temperature t_mean
!-----------------------------------------------------------------------

!***********
!       change in version 3: density is now initialized with unperturbed background temperature,
!                            temperature perturbation is added afterwards
!***********
  rho = p/(Rd*t_mean)

!-----------------------------------------------------------------------
!    POTENTIAL TEMPERATURE PERTURBATION,
!    here: converted to temperature and added to the temperature field
!    models with a prognostic potential temperature field can utilize
!    the potential temperature perturbation theta_pert directly and add it
!    to the background theta field (not included here)
!-----------------------------------------------------------------------

  sin_tmp = sin(lat) * sin(phic)
  cos_tmp = cos(lat) * cos(phic)

    ! great circle distance with 'a/X'

  r  = as * ACOS (sin_tmp + cos_tmp*cos(lon-lambdac))

  s = (d**2)/(d**2 + r**2)

  theta_pert = delta_theta*s*sin(2.d0*pi*height/Lz)

  t_pert = theta_pert*(p/p0)**(Rd/cp)

  t = t_mean + t_pert

!-----------------------------------------------------------------------
!     initialize Q, set to zero
!-----------------------------------------------------------------------

  q = 0.d0

END SUBROUTINE test3_gravity_wave

END MODULE dcmip_initial_conditions_test_1_2_3

max_stars_repo_name: hisashiyashiro/nicamdc

max_issues_repo_name: nicamdc-dev/nicamdc

max_forks_repo_name: nicamdc-dev/nicamdc


text:
      subroutine wrcrc(r)
C...wrcrc      01/06/75   01/01/80
      common/gfgel/gf(28)
      common/gfbftr/xtr,ytr
      gf(10)=100.
      i=abs(gf(9))
      gf(i)=2.
      gf(i+1)=xtr
      gf(i+2)=ytr
      gf(i+3)=r
      gf(i+6)=0.
      gf(i+7)=360.
      if(gf(9).gt.0.)call strmod(0)
      return
      end

ext: f

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
!
! TauCoeff_Compare
!
! Program to compare TauCoeff data read from either netCDF or
! CRTM Binary format files.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 02-Apr-2003
!                       paul.vandelst@ssec.wisc.edu
!

PROGRAM TauCoeff_Compare

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module usage
  USE File_Utility,       ONLY: File_Exists
  USE Message_Handler,    ONLY: SUCCESS, FAILURE, INFORMATION, &
                                Program_Message, Display_Message
  USE TauCoeff_Define,    ONLY: TauCoeff_type, &
                                Destroy_TauCoeff, Equal_TauCoeff
  USE TauCoeff_Binary_IO, ONLY: Read_TauCoeff_Binary
  USE TauCoeff_netCDF_IO, ONLY: Read_TauCoeff_netCDF
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'TauCoeff_Compare'
  CHARACTER( * ), PARAMETER :: PROGRAM_RCS_ID = &

  ! Allowable formats.
  ! Format codes MUST start at 1 and increment by 1.
  INTEGER, PARAMETER :: N_FORMATS = 2
  INTEGER, PARAMETER :: NETCDF_FORMAT = 1
  INTEGER, PARAMETER :: BINARY_FORMAT = 2
  INTEGER, PARAMETER, DIMENSION(N_FORMATS) :: &
    FORMAT_CODE = (/ NETCDF_FORMAT, &
                     BINARY_FORMAT /)
  CHARACTER(*), PARAMETER, DIMENSION(N_FORMATS) :: &
    FORMAT_NAME = (/ 'netCDF', &
                     'Binary' /)

  ! ---------
  ! Variables
  ! ---------
  INTEGER :: Error_Status
  CHARACTER(256) :: File1_Name
  CHARACTER(256) :: File2_Name
  INTEGER :: File1_Format
  INTEGER :: File2_Format
  TYPE(TauCoeff_type) :: TauCoeff1
  TYPE(TauCoeff_type) :: TauCoeff2

  ! Output prgram header
  CALL Program_Message( PROGRAM_NAME, &
                        'Program to compare TauCoeff data read from either '//&
                        'netCDF or CRTM Binary format files.', &
                        '$Revision$' )

  ! Get the filenames and formats
  CALL Get_File_Name_and_Format( 'FIRST',  File1_Name, File1_Format )
  CALL Get_File_Name_and_Format( 'SECOND', File2_Name, File2_Format )

  ! Read the files
  CALL Read_the_File( 'FIRST',  File1_Name, File1_Format, TauCoeff1 )
  CALL Read_the_File( 'SECOND', File2_Name, File2_Format, TauCoeff2 )

  ! Compare the TauCoeff structures
  WRITE( *, '( /5x, "Comparing the two TauCoeff structures ..." )' )
  Error_Status = Equal_TauCoeff( TauCoeff1, TauCoeff2 )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Differences found in TauCoeff structure comparison.', &
                          Error_Status )
  ELSE
    CALL Display_Message( PROGRAM_NAME, &
                          'TauCoeff structures are equal.', &
                          INFORMATION )
  END IF

  ! Destroy the structures
  Error_Status = Destroy_TauCoeff( TauCoeff1 )
  Error_Status = Destroy_TauCoeff( TauCoeff2 )

CONTAINS

  SUBROUTINE Get_File_Name_and_Format(Rank, Filename, FileFormat)
    CHARACTER(*), INTENT(IN)  :: Rank
    CHARACTER(*), INTENT(OUT) :: Filename
    INTEGER,      INTENT(OUT) :: FileFormat
    INTEGER :: n, IO_Status

    ! Get the filename
    WRITE(*, FMT    ='( /5x, "Enter the ",a," TauCoeff file: " )', &
             ADVANCE='NO' ) TRIM(Rank)
    READ(*, '(a)') Filename
    Filename = ADJUSTL(Filename)
    IF ( .NOT. File_Exists( TRIM(Filename) ) ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'File '//TRIM(Filename)//' not found.', &
                            FAILURE )
      STOP
    END IF

    ! Get the file format type
    WRITE(*, FMT='(/5x, "Select the format type for ", a)') TRIM(Filename)
    DO n = 1, N_FORMATS
      WRITE(*, FMT='(10x, i1, ") ", a)') FORMAT_CODE(n), &
                                         TRIM(FORMAT_NAME(n))
    END DO
    WRITE(*, FMT='(5x, "Enter choice: ")', &
             ADVANCE='NO')
    READ(*, FMT   ='(i10)', &
            IOSTAT=IO_Status) FileFormat

    ! Invalid input
    IF ( IO_Status /= 0 ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Input error for '//TRIM(Filename)//' format type.', &
                            FAILURE )
      STOP
    END IF

    ! Invalid value
    IF ( .NOT. ANY( FORMAT_CODE == FileFormat ) ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Invalid format type for '//TRIM(Filename)//'.', &
                            FAILURE )
      STOP
    END IF
  END SUBROUTINE Get_File_Name_and_Format

  SUBROUTINE Read_the_File(Rank, Filename, FileFormat, TauCoeff)
    CHARACTER(*),        INTENT(IN)     :: Rank
    CHARACTER(*),        INTENT(IN)     :: Filename
    INTEGER,             INTENT(IN)     :: FileFormat
    TYPE(TauCoeff_type), INTENT(IN OUT) :: TauCoeff
    INTEGER :: Error_Status

    WRITE(*, '(/5x, "Reading ",a," TauCoeff datafile ",a," ...")') TRIM(Rank), TRIM(Filename)

    SELECT CASE ( FileFormat )
      CASE ( NETCDF_FORMAT )
        Error_Status = Read_TauCoeff_netCDF( Filename, &
                                             TauCoeff )
      CASE ( BINARY_FORMAT )
        Error_Status = Read_TauCoeff_Binary( Filename, &
                                             TauCoeff )
    END SELECT

    IF ( Error_Status /= SUCCESS ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Error reading '//FORMAT_NAME(FileFormat)//&
                            ' TauCoeff file '//TRIM(Filename), &
                            Error_Status )
      STOP
    END IF

  END SUBROUTINE Read_the_File

END PROGRAM TauCoeff_Compare

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
program icdegrade
  implicit none

  integer::i1,i2,i3,i,j,k,narg,iargc,i_file
  integer(kind=4)::np1,np2,np3
  integer(kind=4)::np1o2,np2o2,np3o2
  real::dx,dx2,x1o,x2o,x3o,astart,omegam,omegav,h0
  real,dimension(:,:,:),allocatable::f,f2
  character*80::input,output
  character*80,dimension(2)::filename
  logical::ok

  narg = iargc()
  IF(narg .NE. 2)THEN
     write(*,*)'You should type: a.out input output'
     write(*,*)'where directory input should contain GRAFIC files'
     write(*,*)'and directory output should be empty'
     STOP
  END IF

  CALL getarg(1,input)
  CALL getarg(2,output)

  !  SAFETY CONDITION
  if (input == output) then
     write(*,*)'If input and output directories are the same'
     write(*,*)'input files will be erased by output ones'
     write(*,*)'so type DIFFERENT directories !!!!'
     stop
  endif

  ! READING INPUT FILES
  filename(1) =TRIM(input)//'/ic_refmap'
  filename(2) =TRIM(output)//'/ic_refmap'
  INQUIRE(file=filename(1),exist=ok)
  if(ok)then

     open(10,file=filename(1),form='unformatted')
     read (10)np1,np2,np3,dx,x1o,x2o,x3o,astart,omegam,omegav,h0
     close(10)
     write(*,*)'Input array size is:',np1,np2,np3
     allocate(f(np1,np2,1:2))
     allocate(f2(np1/2,np2/2,1:1))

     np1o2=np1/2
     np2o2=np2/2
     np3o2=np3/2
     dx2=2.*dx

     write(*,*)'Reading input file '//TRIM(filename(1))
     open(10,file=filename(1),form='unformatted')
     rewind(10)
     read (10)np1,np2,np3,dx,x1o,x2o,x3o,astart,omegam,omegav,h0

     write(*,*)'Writing ouput file '//TRIM(filename(2))
     open(11,file=filename(2),form='unformatted')
     rewind 11
     write (11)np1o2,np2o2,np3o2,dx2,x1o,x2o,x3o,astart,omegam,omegav,h0

     write(*,*)'Degrading initial conditions...'

     ! Loop over planes
     do i3=1,np3,2

        ! READING INPUT DATA
        read(10) ((f(i1,i2,1),i1=1,np1),i2=1,np2)
        read(10) ((f(i1,i2,2),i1=1,np1),i2=1,np2)

        !  DEGRADING INITIAL CONDITIONS
        do i1=1,np1o2
           do i2=1,np2o2
              i=i1-1
              j=i2-1
              f2(i1,i2,1)=max(f(2*i+1,2*j+1,1),f(2*i+2,2*j+1,1),f(2*i+1,2*j+2,1),f(2*i+1,2*j+1,2), &
                   &          f(2*i+2,2*j+2,1),f(2*i+1,2*j+2,2),f(2*i+2,2*j+1,2),f(2*i+2,2*j+2,2))
           end do
        end do

        !   WRITING OUTPUT FILES
        write(11) ((f2(i1,i2,1),i1=1,np1o2),i2=1,np2o2)

     enddo

     close(10)
     close(11)
     write(*,*)'done'

     deallocate(f)
     deallocate(f2)
  endif

  filename(1) =TRIM(input)//'/ic_pvar_00001'
  filename(2) =TRIM(output)//'/ic_pvar_00001'
  INQUIRE(file=filename(1),exist=ok)
  if(ok)then

     open(10,file=filename(1),form='unformatted')
     read (10)np1,np2,np3,dx,x1o,x2o,x3o,astart,omegam,omegav,h0
     close(10)
     write(*,*)'Input array size is:',np1,np2,np3
     allocate(f(np1,np2,1:2))
     allocate(f2(np1/2,np2/2,1:1))

     np1o2=np1/2
     np2o2=np2/2
     np3o2=np3/2
     dx2=2.*dx

     write(*,*)'Reading input file '//TRIM(filename(1))
     open(10,file=filename(1),form='unformatted')
     rewind(10)
     read (10)np1,np2,np3,dx,x1o,x2o,x3o,astart,omegam,omegav,h0

     write(*,*)'Writing ouput file '//TRIM(filename(2))
     open(11,file=filename(2),form='unformatted')
     rewind 11
     write (11)np1o2,np2o2,np3o2,dx2,x1o,x2o,x3o,astart,omegam,omegav,h0

     write(*,*)'Degrading initial conditions...'

     ! Loop over planes
     do i3=1,np3,2

        ! READING INPUT DATA
        read(10) ((f(i1,i2,1),i1=1,np1),i2=1,np2)
        read(10) ((f(i1,i2,2),i1=1,np1),i2=1,np2)

        !  DEGRADING INITIAL CONDITIONS
        do i1=1,np1o2
           do i2=1,np2o2
              i=i1-1
              j=i2-1
              f2(i1,i2,1)=max(f(2*i+1,2*j+1,1),f(2*i+2,2*j+1,1),f(2*i+1,2*j+2,1),f(2*i+1,2*j+1,2), &
                   &          f(2*i+2,2*j+2,1),f(2*i+1,2*j+2,2),f(2*i+2,2*j+1,2),f(2*i+2,2*j+2,2))
           end do
        end do

        !   WRITING OUTPUT FILES
        write(11) ((f2(i1,i2,1),i1=1,np1o2),i2=1,np2o2)

     enddo

     close(10)
     close(11)
     write(*,*)'done'

     deallocate(f)
     deallocate(f2)

  endif

end program icdegrade

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


program mandel
    implicit none
    !file
    real :: x,y,h,a,b,atemp
    integer :: i,j,k,n,m
    n = 500
    m = 100
    h = 4. / n
    do i=0,n
        do j=0,n
            x = j * h - 2.
            y = i * h - 2.

            a = 0.
            b = 0.

            do k = 1,m
                if (a * a + b * b > 4.) exit
                atemp = a * a - b * b + x
                b = 2 * a * b + y
                a = atemp
            end do
            write (*,'(I4)',advance='no') k
        end do
        write (*,*) ''
    end do
end program


! { dg-do compile }
! { dg-options "-fdec" }
block data
   structure /s1/
   end structure
end block data


text: Pay no attention to the man behind the curtain. He joined the wiki in September of 2005.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE RMESH (NPER, MXND, X, Y, NID, XN, YN, NUID, LXK, KXL,
     &   NXL, LXN, M1, M2, KKK, KKKOLD, NNN, NNNOLD, LLL, LLLOLD,
     &   IAVAIL, NAVAIL, ERR)
C***********************************************************************

C  SUBROUTINE RMESH = GENERATES AN INITIAL LOGICAL RECTANGULAR MESH
C                     WHOSE PERIMETER IS  (X (I), Y (I), I=1, N).

C***********************************************************************

C  VARIABLES USED:
C     X     = X VALUES AROUND THE PERIMETER
C     Y     = Y VALUES AROUND THE PERIMETER
C     NID   = PERIMETER NODE UNIQUE ID'S
C     N     = NUMBER OF PERIMETER NODES
C     M1    = THE NUMBER OF INTERVALS ON THE FIRST SIDE OF THE RECTANGLE
C     IMAP  = CONTROLS THE DEFINITION OF THE VALUES OF THE COORDINATE
C             OF THE INTERIOR NODES.
C           = 1 MEANS SET ALL COORDINATES TO 0.
C           = 2 MEANS SET ALL COORDINATES TO THE CENTROID OF PERIMETER
C           = 3 MEANS APPLY THE UNIT SQUARE TRANSFORMATION FROM
C             W.A.COOK'S PAPER  (THIOKOL REPORT AFRDL - TR - 71 - 51)
C     ERR   = .TRUE. IF ERRORS WERE ENCOUNTERED
C     XN    = GLOBAL X VALUES OF NODES
C     YN    = GLOBAL Y VALUES OF NODES
C     NUID  = GLOBAL NODE UNIQUE IDENTIFIERS
C     LXK   = LINES PER ELEMENT
C     KXL   = ELEMENTS PER LINE
C     NXL   = NODES PER LINE
C     LXN   = LINES PER NODE
C  NOTE:
C     FOR *XN TABLES A NEGATIVE FLAG IN THE FOURTH COLUMN MEANS
C     GO TO THAT ROW FOR A CONTINUATION OF THE LIST.  IN THAT ROW
C     THE FIRST ELEMENT WILL BE NEGATED TO INDICATE THAT THIS IS
C     A CONTINUATION ROW.   (RMESH ITSELF GENERATES NO SUCH NEGATIVES.)
C     A NEGATIVE FLAG IN THE SECOND COLUMN OF THE LXN ARRAY MEANS
C     THAT THIS NODE IS A BOUNDARY NODE.

C***********************************************************************

      DIMENSION X (NPER), Y (NPER), NID (NPER)
      DIMENSION XN (MXND), YN (MXND), NUID (MXND)
      DIMENSION LXK (4, MXND), KXL (2, 3 * MXND)
      DIMENSION NXL (2, 3 * MXND), LXN (4, MXND)

      LOGICAL ERR, NOROOM

C  NOTE: NOROOM SHOULD NEVER BE TRUE WITH CROSS BEING CALLED IN RMESH.
C        THUS,  IT IS NEVER PASSED BACK TO QMESH.

      KKK = M1 * M2 + KKKOLD
      LLL =  (M1 *  (M2 + 1)) +  (M2 *  (M1 + 1)) + LLLOLD
      NNN =  (M1 + 1) *  (M2 + 1) + NNNOLD
      ERR = .TRUE.

C  CHECK INPUT

      IF (2 *  (NPER / 2) .NE. NPER) THEN
         CALL MESAGE ('IN RMESH,  NO. OF PERIMETER NODES IS ODD')
         RETURN
      ELSEIF ( (M1 .LT. 1) .OR. (M2 .LT. 1)) THEN
         WRITE ( * , 10000)NPER, M1
         RETURN
      ENDIF

C  COMPUTE CONSTANTS

      NLC = 2 * M1 + 1
      M1P1 = M1 + 1
      M2P1 = M2 + 1

C  PRODUCE LXK ARRAY

C  LINES FOR FIRST ELEMENT

      LXK (1, KKKOLD + 1) = 1 + LLLOLD
      LXK (2, KKKOLD + 1) = M1 + 1 + LLLOLD
      LXK (3, KKKOLD + 1) = M1 + 2 + LLLOLD
      LXK (4, KKKOLD + 1) = NLC + 1 + LLLOLD

C  FIRST ROW  (SHIFT FIRST ELEMENT TO SECOND,  ETC.)

      IF (M1 .GT. 1) THEN
         DO 110 K = 2, M1
            DO 100 I = 1, 4
               LXK (I, K + KKKOLD) = LXK (I, K + KKKOLD - 1) + 1
  100       CONTINUE
  110    CONTINUE
      ENDIF

C  SUCCEEDING ROWS  (SHIFT FIRST COLUMN TO SECOND, ETC.)

      IF (M2 .GT. 1) THEN
         K = M1 + KKKOLD
         DO 140 K2 = 2, M2
            DO 130 K1 = 1, M1
               K = K + 1
               KL = K - M1
               DO 120 I = 1, 4
                  LXK (I, K) = LXK (I, KL) + NLC
  120          CONTINUE
  130       CONTINUE
  140    CONTINUE
      ENDIF

C  PREPARE KXL TABLE BY USING SUBROUTINE CROSS ON THE LXK TABLE

      CALL CCROSS (4, KKK, 2, LLL, LXK, KXL, KKKOLD + 1, LLLOLD + 1,
     &   NOROOM, ERR)
      IF (ERR) RETURN
      ERR = .TRUE.

C  MAKE NXL TABLE

C  FIRST DO HORIZONTAL LINES

      DO 160 J = 1, M2P1
         NODE = 1 +  (J - 1) *  (M1 + 1) + NNNOLD
         L = 1 +  ( (J - 1) * NLC) + LLLOLD
         DO 150 I = 1, M1
            NXL (1, L) = NODE
            NXL (2, L) = NODE + 1
            NODE = NODE + 1
            L = L + 1
  150    CONTINUE
  160 CONTINUE

C  NEXT DO VERTICAL LINES

      DO 180 J = 1, M1P1
         NODE = J + NNNOLD
         L = J + M1 + LLLOLD
         DO 170 I = 1, M2
            NXL (1, L) = NODE
            NXL (2, L) = NODE + M1P1
            NODE = NODE + M1P1
            L = L + NLC
  170    CONTINUE
  180 CONTINUE

C  PREPARE LXN TABLE FROM NXL TABLE

      CALL CCROSS (2, LLL, 4, NNN, NXL, LXN, LLLOLD + 1, NNNOLD + 1,
     &   NOROOM, ERR)
      IF (ERR) RETURN
      ERR = .TRUE.

C  LINK - UP AVAILABLE LXN SPACE

      IAVAIL = NNN + 1
      NAVAIL = MXND - NNN
      DO 190 I = IAVAIL, MXND
         LXN (1, I) = 0
         LXN (2, I) = 0
         LXN (3, I) = 0
         LXN (4, I) = I + 1
  190 CONTINUE

C  LOGICAL CONNECTION TABLES ARE COMPLETE
C  FILL IN THE CO - ORDINATES OF THE INTERIOR POINTS
C  USE THE UNIT SQUARE TRANSFORMATION  (OF COOK / THIOKOL)

      IF ( (M1 .GT. 1) .AND. (M2 .GT. 1)) THEN

C  GET NODE NUMBERS FOR CORNERS

         I1Z = M1 + 1
         IZ1 = NPER -  (M2 - 1)
         I11 = M1 + M2 + 1

C  COLUMN LOOP

         DO 210 J = 2, M2
            KL = NPER + 2 - J
            KR = M1 + J
            ETA = DBLE(J - 1) / DBLE(M2)
            OMETA = 1.0 - ETA

C  ROW LOOP

            DO 200 I = 2, M1
               KB = I
               KT = IZ1 + 1 - I
               EPS = DBLE(I - 1) / DBLE(M1)
               OMEPS = 1.0 - EPS
               IM =  (J - 1) * M1P1 + I + NNNOLD
               XN (IM) =  (OMETA * X (KB)) + (ETA * X (KT)) +
     &            (OMEPS * X (KL)) + (EPS * X (KR)) -
     &            ((X (1) * OMETA * OMEPS) + (X (I1Z) * OMETA * EPS) +
     &            (X (IZ1) * ETA * OMEPS) + (X (I11) * ETA * EPS))
               YN (IM) =  (OMETA * Y (KB)) + (ETA * Y (KT)) +
     &            (OMEPS * Y (KL)) + (EPS * Y (KR)) -
     &            ((Y (1) * OMETA * OMEPS) + (Y (I1Z) * OMETA * EPS) +
     &            (Y (IZ1) * ETA * OMEPS) + (Y (I11) * ETA * EPS))
  200       CONTINUE
  210    CONTINUE
      ENDIF

C  DEFINE THE COORDINATES OF THE PERIMETER NODES.
C  ALSO FLAG SECOND ELEMENTS OF LXN ARRAY TO INDICATE
C  WHICH NODES ARE BOUNDARY NODES.
C  DEFINE UNIQUE NODE ID NUMBERS ALSO.

      DO 220 I = NNNOLD + 1, NNN
         NUID (I) = 0
  220 CONTINUE

C  BOTTOM

      IP = 0
      DO 230 I = 1, M1P1
         IM = I + NNNOLD
         IP = IP + 1
         LXN (2, IM) =  - LXN (2, IM)
         NUID (IM) = NID (IP)
         XN (IM) = X (IP)
         YN (IM) = Y (IP)
  230 CONTINUE

C  RIGHT

      IP = M1P1
      DO 240 I = 2, M2P1
         IP = IP + 1
         IM = IM + M1P1
         LXN (2, IM) =  - LXN (2, IM)
         NUID (IM) = NID (IP)
         XN (IM) = X (IP)
         YN (IM) = Y (IP)
  240 CONTINUE

C  TOP

      DO 250 I = 2, M1P1
         IP = IP + 1
         IM = IM - 1
         LXN (2, IM) =  - LXN (2, IM)
         NUID (IM) = NID (IP)
         XN (IM) = X (IP)
         YN (IM) = Y (IP)
  250 CONTINUE

C  LEFT

      DO 260 I = 2, M2
         IP = IP + 1
         IM = IM - M1P1
         LXN (2, IM) =  - LXN (2, IM)
         NUID (IM) = NID (IP)
         XN (IM) = X (IP)
         YN (IM) = Y (IP)
  260 CONTINUE

C  EXIT

      ERR = .FALSE.

      RETURN

10000 FORMAT (' IN RMESH, N = ', I5, ' AND M1 = ', I5,
     &   ' ARE INCOMPATIBLE')

      END

max_stars_repo_name: utkarshayachit/seacas

max_issues_repo_name: utkarshayachit/seacas

max_forks_repo_name: utkarshayachit/seacas


text:
C
C $Id: line3w.f,v 1.5 2008-07-27 00:17:34 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE LINE3W (XA,YA,XB,YB)
      SAVE
      COMMON /PRM31/  Q          ,L
      GO TO ( 10, 30, 40),L
   10 UA = Q
      UB = Q
      VA = XA
      VB = XB
   20 WA = YA
      WB = YB
      GO TO  50
   30 UA = XA
      UB = XB
      VA = Q
      VB = Q
      GO TO  20
   40 UA = XA
      UB = XB
      VA = YA
      VB = YB
      WA = Q
      WB = Q
   50 CALL LINE3 (UA,VA,WA,UB,VB,WB)
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
PROGRAM obsmake
!=======================================================================
!
! [PURPOSE:] For generating observations
!
! [HISTORY:]
!   01/26/2009 Takemasa Miyoshi  created
!
!=======================================================================
!$USE OMP_LIB
  USE common
  USE common_speedy
  USE common_obs_speedy

  IMPLICIT NONE

  LOGICAL,PARAMETER :: verbose=.FALSE.

  INTEGER,PARAMETER :: obkind=5
  REAL(r_size),ALLOCATABLE :: v3d(:,:,:,:)
  REAL(r_size),ALLOCATABLE :: v2d(:,:,:)
  REAL(r_size),ALLOCATABLE :: error(:)
  REAL(r_sngl) :: wk(6)
  REAL(r_size) :: oberr(obkind)
  INTEGER :: obelm(obkind)
  LOGICAL :: ex
  REAL(r_size) :: a
  INTEGER :: i,j,k,n,nn,ios
  INTEGER :: nobkind,nstation
  INTEGER :: nobs
  CHARACTER(100) :: cdummy

  CHARACTER(8) :: truef='true.grd'
  CHARACTER(7) :: obsf='obs.dat'
  !
  ! Initial settings
  !
  CALL set_common_speedy
  !
  ! Read true (nature run)
  !
  ALLOCATE(v3d(nlon,nlat,nlev,nv3d))
  ALLOCATE(v2d(nlon,nlat,nv2d))
  CALL read_grd(truef,v3d,v2d)
  !
  ! read table obserr.tbl
  !
  OPEN(11,FILE='obserr.tbl')
  READ(11,'(A)') cdummy
  READ(11,'(A)') cdummy
  nobkind = 0
  DO
    READ(11,'(L1,X,I5,X,ES9.2)',IOSTAT=ios) ex,i,a
    IF(ios /= 0) EXIT
    IF(ex) THEN
      nobkind = nobkind+1
      obelm(nobkind) = i
      oberr(nobkind) = a
    END IF
  END DO
  CLOSE(11)
  !
  ! Locations of station (station.tbl)
  !
  OPEN(10,FILE='station.tbl')
  READ(10,'(A)') cdummy
  READ(10,'(A)') cdummy
  nstation = 0
  DO
    READ(10,'(2I3)',IOSTAT=ios) i,j
    IF(ios /= 0) EXIT
    nstation = nstation + 1
  END DO
  CLOSE(10)
  IF(verbose) PRINT '(A,I)','nstation = ',nstation
  !
  ! Count number of obs
  !
  nobs = 0
  DO n=1,nobkind
    IF(obelm(n) == id_ps_obs) THEN
      nobs = nobs + nstation
    ELSE
      nobs = nobs + nstation*nlev
    END IF
  END DO
  IF(verbose) PRINT '(A,I)','nobs = ',nobs
  !
  ! Random number
  !
  ALLOCATE(error(nobs))
  CALL com_randn(nobs,error)
  !
  ! Output OBS
  !
  nn = 0
  OPEN(90,FILE=obsf,FORM='unformatted',ACCESS='sequential')
  OPEN(10,FILE='station.tbl')
  READ(10,'(A)') cdummy
  READ(10,'(A)') cdummy
  DO
    READ(10,'(2I3)',IOSTAT=ios) i,j
    IF(ios /= 0) EXIT
    DO n=1,nobkind
      IF(obelm(n) == id_ps_obs) THEN
        wk(1)=REAL(obelm(n),r_sngl)
        wk(2)=REAL(lon(i),r_sngl)
        wk(3)=REAL(lat(j),r_sngl)
        wk(4)=REAL(phi0(i,j),r_sngl)
        nn = nn+1
        wk(5)=REAL(v2d(i,j,iv2d_ps)+error(nn)*oberr(n),r_sngl)/100.0 ! hPa
        wk(6)=REAL(oberr(n),r_sngl)/100.0 ! hPa
        WRITE(90) wk
      ELSE
        wk(1)=REAL(obelm(n),r_sngl)
        wk(2)=REAL(lon(i),r_sngl)
        wk(3)=REAL(lat(j),r_sngl)
        wk(6)=REAL(oberr(n),r_sngl)
        DO k=1,nlev
          wk(4)=REAL(v2d(i,j,iv2d_ps)*sig(k),r_sngl)/100.0 ! hPa
          nn = nn+1
          SELECT CASE(obelm(n))
          CASE(id_u_obs)
            wk(5)=REAL(v3d(i,j,k,iv3d_u)+error(nn)*oberr(n),r_sngl)
          CASE(id_v_obs)
            wk(5)=REAL(v3d(i,j,k,iv3d_v)+error(nn)*oberr(n),r_sngl)
          CASE(id_t_obs)
            wk(5)=REAL(v3d(i,j,k,iv3d_t)+error(nn)*oberr(n),r_sngl)
          CASE(id_q_obs)
            wk(5)=REAL(v3d(i,j,k,iv3d_q)+error(nn)*oberr(n),r_sngl)
          END SELECT
          IF(obelm(n) == id_q_obs .AND. k >= 5) CYCLE !No upper-level Q obs
          WRITE(90) wk
        END DO
      END IF
    END DO
  END DO
  CLOSE(90)
  !
  ! Count number of obs
  !
  IF(verbose) THEN
    CALL get_nobs(obsf,6,nn)
  END IF

  DEALLOCATE(v3d,v2d,error)

  STOP
END PROGRAM obsmake

max_stars_repo_name: skotsuki/speedy-lpf

max_issues_repo_name: skotsuki/speedy-lpf

max_forks_repo_name: skotsuki/speedy-lpf


text:
      subroutine loop
      integer i, j, k

c     check points
      integer a1, a2, a3

      i = 1
      j = 123
 10   if (i .gt. 10) goto 40
      k = 456
      i = i + 1
      goto 10

 40   a1 = i     ! Bot
      a2 = j     ! Const (Int 123)
      a3 = k     ! Const (Int 456)
      end

      subroutine doloop
      integer i, j, k
      integer b(10)

c     check points
      integer b1, b2, b3, b4
      j = 123
      do 99  i = 2, 10
        b(i) = i
        k = 456
 99   continue

      b1 = b(1)   ! Bot
      b2 = b(2)   ! Bot
      b3 = j      ! Const (Int 123)
      b4 = k      ! Const (Int 456)
      end

max_stars_repo_name: camfort/fortran-vars

max_issues_repo_name: camfort/fortran-vars

max_forks_repo_name: camfort/fortran-vars


text:
      program time_shift

!     Time-shifts a SAC file.  One is prompted for the input & output
!     filenames and tshift (new-old).
!     The program calls SAC library routines for forward and inverse
!        Fourier transforms.

!   To compile:
!     gfortran -o time_shift time_shift.f ‘sac-config --cflags --libs libsac libsacio

!   Sample run format: (input.sac and a time-shift (new-old) of -0.05 are input.)
!     time_shift input.sac output.sac -0.05    -

      implicit none
      integer nmax
      parameter(nmax = 131072)

      real*4 signal(nmax),re(nmax),im(nmax)
      character filename_in*80, filename_out*80, name*20
      real*4 dt,tshift,secs,b
      integer npts, nmarg, nerr, nfn
      integer lenc
C
C
      nmarg = iargc()
      if (nmarg .eq. 0) then
         write(*,*)'Usage: time_shift filename_in filename_out tshift'
         write(*,*)'    tshift is new-old'
         stop
      end if
      call getarg(1,filename_in)
      call getarg(2,filename_out)
      call getarg(3,name)
      read(name,'(f10.0)') tshift
C
      nfn = lenc(filename_in)
      call rsac1(filename_in(1:nfn),signal,npts,secs,dt,nmax,nerr)
      if (nerr .ne. 0) then
         write(*,*) 'Error opening file ',filename_in(1:nfn)
         stop
      end if

      ! Take out mean and trend and then taper the data
      call getfhv('b',b,nerr)
      call remove_trend(signal, npts, dt, b)
      call taper_width(signal, npts, 2, 0.05)

      call timeshift(signal,re,im,npts,nmax,dt,tshift)

      write(*,'(a,f10.3)') 'Time shift (new-old) of',tshift

      call setfhv('user9',tshift,nerr)

      ! Write out timeshifted file
      nfn = lenc(filename_out)
      call wsac0(filename_out(1:nfn),signal,signal,nerr)
      stop
      end

      subroutine timeshift(signal,re,im,npts,nmax,dt,tshift)

!     time shifts signal by tshift (new-old)

      implicit none

      ! Input Parameters
      real*4 signal(*),re(*),im(*)
      integer npts, nf      ! Length of signal
      integer nmax
      real dt, tshift

      !  real*4 array signal is changed.

      integer k, ntran
      real ttot, df,fnyq

      ttot  = dt*(npts-1)
      ntran = min(4*npts,nmax)
      write(*,*) "dt, npts, ntran", dt, npts, ntran

      ! Pad Input with Zeros
      if (npts < ntran) then
         do k=npts+1,ntran
            signal(k) = 0.0
         end do
      end if

! nf must be a power of 2 that is >= ntran

      nf = 4
      do while (nf < ntran)
        nf = nf * 2
      enddo

      ! Forward FFT  signal, ntran, and nf input; re and im are output
      call fft(signal, ntran, re, im, nf)   !  SAC library routine
      df    = 1./(nf*dt)
      fnyq = (nf/2)*df
      write(*,*) "nf, df, fnyq", nf,df,fnyq

      ! Time shift
      if (abs(tshift) .ge. 1e-8) then
         call shiftt(nf,re,im,df,tshift)
      end if

      ! Inverse FFT  nf, re, im, nf input;  signal output
      call ifft(signal, ntran, re, im, nf)
      write(*,*) 'ntran nf',ntran, nf
      return
      end

      SUBROUTINE SHIFTT(nf,re,im,DF,TSHIFT)

      COMPLEX FSIG
!      complex  now
      real*4 re(*),im(*)
      TWOPI = 8.0*ATAN(1.0)
      N = nf/2 + 1
      DO J=1,N
         fsig = cmplx(re(j),im(j))
!         now = fsig
!        write(*,*) df*(j-1),real(fsig(j)),aimag(fsig(j))
         WT = -TWOPI*(J-1)*DF*TSHIFT    !SAC default sign convention
         FSIG = FSIG*CMPLX(COS(WT),SIN(WT))
!        write(*,*) j, df*(j-1), now, fsig
         re(j) = real(fsig)
         im(j) = aimag(fsig)
      END DO
!      if (ntran.gt.0) stop
      do j=N+1,nf
        re(j) = re(nf-j+2)
        im(j) = -im(nf-j+2)
      enddo
      RETURN
      END

      function lenc(string)

C     Returns length of character variable STRING excluding right-hand
C     most blanks or nulls

      character*(*) string
      length = len(string)      ! total length
      if (length .eq. 0) then
         lenc = 0
         return
      end if
      if(ichar(string(length:length)).eq.0)string(length:length) = ' '
      do j=length,1,-1
         lenc = j
         if (string(j:j).ne.' ' .and. ichar(string(j:j)).ne.0) return
      end do
      lenc = 0
      return
      end


text:
program MC_ex01_p01
    implicit none
    real,parameter :: pi=4.0*atan(1.0)
    integer::seed=13,iarg,i,j,k,n,hits
    real(8)::results(7,100),l=2,d=10,errors(7)=0
    character(len=80)::arg

    !RNG seed from arguments
    iarg=command_argument_count()
    if (iarg==1) then
        call get_command_argument(1,arg)
        read(arg,*) seed
    end if

    call srand(seed)

    !Simulation part
    n=1
    do i=1,7
        !Update the amount of throws for next simulations
        n=n*10
        !Do 100 simulations with the new throw count
        do j=1,100
            !Simulate a single set of throws
            hits=0
            do k=1,n
                hits=hits+throw()
            end do
            if(hits==0) then
                results(i,j)=0
            else
                results(i,j)=n*2*l/(d*hits)
            end if
            errors(i)=errors(i)+abs(results(i,j)-pi)
        end do
        errors(i)=errors(i)/100
    end do

    !Saving and plotting of the error data
    open(1,action="write",file="data.txt",status="replace")
    do i=1,7
        write(1,*) i, errors(i)
    end do
    close(1)
    call execute_command_line("python plot.py")

contains
    integer function throw()
        implicit none
        real(8)::y,alpha
        !Randomize the position and the angle of the needle.
        !It's only necessary to look at the other half of the throwing space, since the situation can always be mirrored
        y=rand()*d/2
        alpha=rand()*pi/2
        !Check if the end of the needle crosses a line
        if (l/2*sin(alpha)>=y) then
            throw=1
        else
            throw=0
        end if
    end function

end program MC_ex01_p01

max_stars_repo_name: toicca/basics-of-mc

max_issues_repo_name: toicca/basics-of-mc

max_forks_repo_name: toicca/basics-of-mc


text:
 MODULE GDSWZD_GAUSSIAN_MOD
!$$$  MODULE DOCUMENTATION BLOCK
!
! $Revision: 74917 $
!
! MODULE:  GDSWZD_GAUSSIAN_MOD  GDS WIZARD MODULE FOR GAUSSIAN CYLINDRICAL
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: - CONVERT FROM EARTH TO GRID COORDINATES OR VICE VERSA.
!           - COMPUTE VECTOR ROTATION SINES AND COSINES.
!           - COMPUTE MAP JACOBIANS.
!           - COMPUTE GRID BOX AREA.
!
! PROGRAM HISTORY LOG:
!   2015-01-21  GAYNO   INITIAL VERSION FROM A MERGER OF
!                       ROUTINES GDSWIZ04 AND GDSWZD04.
!   2015-09-17  GAYNO   RENAME MODULE AS "GDSWZD_GAUSSIAN_MOD"
!
! USAGE:  "USE GDSWZD_GAUSSIAN_MOD"  THEN CALL THE PUBLIC DRIVER
!         ROUTINE "GDSWZD_GAUSSIAN".
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
!
 IMPLICIT NONE

 PRIVATE

 PUBLIC                         :: GDSWZD_GAUSSIAN

 REAL,            PARAMETER     :: PI=3.14159265358979
 REAL,            PARAMETER     :: DPR=180./PI

 INTEGER                        :: J1, JH

 REAL,            ALLOCATABLE   :: BLAT(:)
 REAL                           :: DLON, RERTH
 REAL,            ALLOCATABLE   :: YLAT_ROW(:)

 CONTAINS

 SUBROUTINE GDSWZD_GAUSSIAN(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, &
                            XPTS,YPTS,RLON,RLAT,NRET, &
                            CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  GDSWZD_GAUSSIAN   GDS WIZARD FOR GAUSSIAN CYLINDRICAL
!   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
!
! ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB 2 GRID DEFINITION
!           TEMPLATE (PASSED IN INTEGER FORM AS DECODED BY THE
!           NCEP G2 LIBRARY) AND RETURNS ONE OF THE FOLLOWING:
!             (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES
!             (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES
!           WORKS FOR GAUSSIAN CYLINDRICAL PROJECTIONS.
!           IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT
!           BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT
!           OUTPUT ELEMENTS ARE SET TO FILL VALUES.
!           THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO.
!           OPTIONALLY, THE VECTOR ROTATIONS, THE MAP JACOBIANS AND
!           THE GRID BOX AREAS MAY BE RETURNED AS WELL.  TO COMPUTE
!           THE VECTOR ROTATIONS, THE OPTIONAL ARGUMENTS 'SROT' AND 'CROT'
!           MUST BE PRESENT.  TO COMPUTE THE MAP JACOBIANS, THE
!           OPTIONAL ARGUMENTS 'XLON', 'XLAT', 'YLON', 'YLAT' MUST BE PRESENT.
!           TO COMPUTE THE GRID BOX AREAS, THE OPTIONAL ARGUMENT
!           'AREA' MUST BE PRESENT.
!
! PROGRAM HISTORY LOG:
!   96-04-10  IREDELL
!   97-10-20  IREDELL  INCLUDE MAP OPTIONS
! 1999-04-08  IREDELL  USE SUBROUTINE SPLAT
! 2001-06-18  IREDELL  CORRECT AREA COMPUTATION
! 2012-08-01  GAYNO    CORRECT AREA COMPUTATION AT POLE.
!                      CORRECT YLAT COMPUTATION.
! 2015-01-21  GAYNO    MERGER OF GDSWIZ04 AND GDSWZD04.  MAKE
!                      CROT,SORT,XLON,XLAT,YLON,YLAT AND AREA
!                      OPTIONAL ARGUMENTS.  MAKE PART OF A MODULE.
!                      MOVE VECTOR ROTATION, MAP JACOBIAN AND GRID
!                      BOX AREA COMPUTATIONS TO SEPARATE SUBROUTINES.
! 2015-07-13  GAYNO    CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAY
!                      WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAY.
!                      RENAME AS "GDSWZD_GAUSSIAN".
!
! USAGE:    CALL GDSWZD_GAUSSIAN(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL,
!    &                           XPTS,YPTS,RLON,RLAT,NRET,
!    &                           CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
!
!   INPUT ARGUMENT LIST:
!     IGDTNUM  - INTEGER GRID DEFINITION TEMPLATE NUMBER.
!                CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                MUST BE "40" FOR GAUSSIAN GRIDS.
!     IGDTMPL  - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY.
!                CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE FOR SECTION
!                THREE:
!                 (1):  SHAPE OF EARTH, OCTET 15
!                 (2):  SCALE FACTOR OF SPHERICAL EARTH RADIUS,
!                       OCTET 16
!                 (3):  SCALED VALUE OF RADIUS OF SPHERICAL EARTH,
!                       OCTETS 17-20
!                 (4):  SCALE FACTOR OF MAJOR AXIS OF ELLIPTICAL EARTH,
!                       OCTET 21
!                 (5):  SCALED VALUE OF MAJOR AXIS OF ELLIPTICAL EARTH,
!                       OCTETS 22-25
!                 (6):  SCALE FACTOR OF MINOR AXIS OF ELLIPTICAL EARTH,
!                       OCTET 26
!                 (7):  SCALED VALUE OF MINOR AXIS OF ELLIPTICAL EARTH,
!                       OCTETS 27-30
!                 (8):  NUMBER OF POINTS ALONG A PARALLEL, OCTS 31-34
!                 (9):  NUMBER OF POINTS ALONG A MERIDIAN, OCTS 35-38
!                 (10): BASIC ANGLE OF INITIAL PRODUCTION DOMAIN,
!                       OCTETS 39-42
!                 (11): SUBDIVISIONS OF BASIC ANGLE, OCTETS 43-46
!                 (12): LATITUDE OF FIRST GRID POINT, OCTETS 47-50
!                 (13): LONGITUDE OF FIRST GRID POINT, OCTETS 51-54
!                 (14): RESOLUTION AND COMPONENT FLAGS, OCTET 55
!                 (15): LATITUDE OF LAST GRID POINT, OCTETS 56-59
!                 (16): LONGITUDE OF LAST GRID POINT, OCTETS 60-63
!                 (17): I-DIRECTION INCREMENT, OCTETS 64-67
!                 (18): NUMBER OF PARALLELS BETWEEN POLE AND EQUATOR,
!                       OCTETS 68-71
!                 (19): SCANNING MODE, OCTET 72
!     IGDTLEN  - INTEGER NUMBER OF ELEMENTS (19) OF THE GRID DEFINITION
!                TEMPLATE ARRAY.  CORRESPONDS TO THE GFLD%IGDTLEN
!                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!     IOPT     - INTEGER OPTION FLAG
!                (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
!                (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
!     NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
!     FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
!                (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
!     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
!     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
!     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
!                (ACCEPTABLE RANGE: -360. TO 360.)
!     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
!                (ACCEPTABLE RANGE: -90. TO 90.)
!
!   OUTPUT ARGUMENT LIST:
!     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
!     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
!     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
!     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
!     NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
!     CROT     - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES
!     SROT     - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION SINES
!                (UGRID=CROT*UEARTH-SROT*VEARTH;
!                 VGRID=SROT*UEARTH+CROT*VEARTH)
!     XLON     - REAL, OPTIONAL (NPTS) DX/DLON IN 1/DEGREES
!     XLAT     - REAL, OPTIONAL (NPTS) DX/DLAT IN 1/DEGREES
!     YLON     - REAL, OPTIONAL (NPTS) DY/DLON IN 1/DEGREES
!     YLAT     - REAL, OPTIONAL (NPTS) DY/DLAT IN 1/DEGREES
!     AREA     - REAL, OPTIONAL (NPTS) AREA WEIGHTS IN M**2
!
! EXTERNAL SUBPROGRAMS CALLED:
!   IP_SPLAT      COMPUTE LATITUDE FUNCTIONS
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE
!
 INTEGER,         INTENT(IN   ) :: IGDTNUM, IGDTLEN
 INTEGER,         INTENT(IN   ) :: IGDTMPL(IGDTLEN)
 INTEGER,         INTENT(IN   ) :: IOPT, NPTS
 INTEGER,         INTENT(  OUT) :: NRET
!
 REAL,            INTENT(IN   ) :: FILL
 REAL,            INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
 REAL,            INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
 REAL, OPTIONAL,  INTENT(  OUT) :: CROT(NPTS),SROT(NPTS)
 REAL, OPTIONAL,  INTENT(  OUT) :: XLON(NPTS),XLAT(NPTS)
 REAL, OPTIONAL,  INTENT(  OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
!
 INTEGER                        :: ISCAN, JSCAN, IM, JM
 INTEGER                        :: J, JA, JG
 INTEGER                        :: ISCALE, N
!
 LOGICAL                        :: LROT, LMAP, LAREA
!
 REAL,            ALLOCATABLE   :: ALAT(:), ALAT_JSCAN(:)
 REAL,            ALLOCATABLE   :: ALAT_TEMP(:),BLAT_TEMP(:)
 REAL                           :: HI, RLATA, RLATB, RLAT1, RLON1, RLON2
 REAL                           :: XMAX, XMIN, YMAX, YMIN, YPTSA, YPTSB
 REAL                           :: ECCEN_SQUARED, WB
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 IF(PRESENT(CROT)) CROT=FILL
 IF(PRESENT(SROT)) SROT=FILL
 IF(PRESENT(XLON)) XLON=FILL
 IF(PRESENT(XLAT)) XLAT=FILL
 IF(PRESENT(YLON)) YLON=FILL
 IF(PRESENT(YLAT)) YLAT=FILL
 IF(PRESENT(AREA)) AREA=FILL
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! IS THIS A GAUSSIAN GRID?
 IF(IGDTNUM/=40) THEN
   CALL GAUSSIAN_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
   RETURN
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 CALL EARTH_RADIUS(IGDTMPL,IGDTLEN,RERTH,ECCEN_SQUARED)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! ENSURE RADIUS OF EARTH IS DEFINED.
 IF(RERTH<0.) THEN
   CALL GAUSSIAN_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
   RETURN
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 IF(PRESENT(CROT).AND.PRESENT(SROT))THEN
   LROT=.TRUE.
 ELSE
   LROT=.FALSE.
 ENDIF
 IF(PRESENT(XLON).AND.PRESENT(XLAT).AND.PRESENT(YLON).AND.PRESENT(YLAT))THEN
   LMAP=.TRUE.
 ELSE
   LMAP=.FALSE.
 ENDIF
 IF(PRESENT(AREA))THEN
   LAREA=.TRUE.
 ELSE
   LAREA=.FALSE.
 ENDIF
 IM=IGDTMPL(8)
 JM=IGDTMPL(9)
 ISCALE=IGDTMPL(10)*IGDTMPL(11)
 IF(ISCALE==0) ISCALE=10**6
 RLAT1=FLOAT(IGDTMPL(12))/FLOAT(ISCALE)
 RLON1=FLOAT(IGDTMPL(13))/FLOAT(ISCALE)
 RLON2=FLOAT(IGDTMPL(16))/FLOAT(ISCALE)
 JG=IGDTMPL(18)*2
 ISCAN=MOD(IGDTMPL(19)/128,2)
 JSCAN=MOD(IGDTMPL(19)/64,2)
 HI=(-1.)**ISCAN
 JH=(-1)**JSCAN
 DLON=HI*(MOD(HI*(RLON2-RLON1)-1+3600,360.)+1)/(IM-1)
 ALLOCATE(ALAT_TEMP(JG))
 ALLOCATE(BLAT_TEMP(JG))
 CALL IP_SPLAT(4,JG,ALAT_TEMP,BLAT_TEMP)
 ALLOCATE(ALAT(0:JG+1))
 ALLOCATE(BLAT(0:JG+1))
!$OMP PARALLEL DO PRIVATE(JA) SCHEDULE(STATIC)
 DO JA=1,JG
   ALAT(JA)=DPR*ASIN(ALAT_TEMP(JA))
   BLAT(JA)=BLAT_TEMP(JA)
 ENDDO
!$OMP END PARALLEL DO
 DEALLOCATE(ALAT_TEMP,BLAT_TEMP)
 ALAT(0)=180.-ALAT(1)
 ALAT(JG+1)=-ALAT(0)
 BLAT(0)=-BLAT(1)
 BLAT(JG+1)=BLAT(0)
 J1=1
 DO WHILE(J1.LT.JG.AND.RLAT1.LT.(ALAT(J1)+ALAT(J1+1))/2)
   J1=J1+1
 ENDDO
 IF(LMAP)THEN
   ALLOCATE(ALAT_JSCAN(JG))
   DO JA=1,JG
     ALAT_JSCAN(J1+JH*(JA-1))=ALAT(JA)
   ENDDO
   ALLOCATE(YLAT_ROW(0:JG+1))
   DO JA=2,(JG-1)
     YLAT_ROW(JA)=2.0/(ALAT_JSCAN(JA+1)-ALAT_JSCAN(JA-1))
   ENDDO
   YLAT_ROW(1)=1.0/(ALAT_JSCAN(2)-ALAT_JSCAN(1))
   YLAT_ROW(0)=YLAT_ROW(1)
   YLAT_ROW(JG)=1.0/(ALAT_JSCAN(JG)-ALAT_JSCAN(JG-1))
   YLAT_ROW(JG+1)=YLAT_ROW(JG)
   DEALLOCATE(ALAT_JSCAN)
 ENDIF
 XMIN=0
 XMAX=IM+1
 IF(IM.EQ.NINT(360/ABS(DLON))) XMAX=IM+2
 YMIN=0.5
 YMAX=JM+0.5
 NRET=0
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
 IF(IOPT.EQ.0.OR.IOPT.EQ.1) THEN
!$OMP PARALLEL DO PRIVATE(N,J,WB,RLATA,RLATB) REDUCTION(+:NRET) SCHEDULE(STATIC)
   DO N=1,NPTS
     IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND. &
        YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
       RLON(N)=MOD(RLON1+DLON*(XPTS(N)-1)+3600,360.)
       J=YPTS(N)
       WB=YPTS(N)-J
       RLATA=ALAT(J1+JH*(J-1))
       RLATB=ALAT(J1+JH*J)
       RLAT(N)=RLATA+WB*(RLATB-RLATA)
       NRET=NRET+1
       IF(LROT) CALL GAUSSIAN_VECT_ROT(CROT(N),SROT(N))
       IF(LMAP) CALL GAUSSIAN_MAP_JACOB(YPTS(N),&
                                        XLON(N),XLAT(N),YLON(N),YLAT(N))
       IF(LAREA) CALL GAUSSIAN_GRID_AREA(YPTS(N),AREA(N))
     ELSE
       RLON(N)=FILL
       RLAT(N)=FILL
     ENDIF
   ENDDO
!$OMP END PARALLEL DO
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
 ELSEIF(IOPT.EQ.-1) THEN
!$OMP PARALLEL DO PRIVATE(N,JA,YPTSA, YPTSB, WB) REDUCTION(+:NRET) SCHEDULE(STATIC)
   DO N=1,NPTS
     XPTS(N)=FILL
     YPTS(N)=FILL
     IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90) THEN
       XPTS(N)=1+HI*MOD(HI*(RLON(N)-RLON1)+3600,360.)/DLON
       JA=MIN(INT((JG+1)/180.*(90-RLAT(N))),JG)
       IF(RLAT(N).GT.ALAT(JA)) JA=MAX(JA-2,0)
       IF(RLAT(N).LT.ALAT(JA+1)) JA=MIN(JA+2,JG)
       IF(RLAT(N).GT.ALAT(JA)) JA=JA-1
       IF(RLAT(N).LT.ALAT(JA+1)) JA=JA+1
       YPTSA=1+JH*(JA-J1)
       YPTSB=1+JH*(JA+1-J1)
       WB=(ALAT(JA)-RLAT(N))/(ALAT(JA)-ALAT(JA+1))
       YPTS(N)=YPTSA+WB*(YPTSB-YPTSA)
       IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND. &
          YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN
         NRET=NRET+1
         IF(LROT) CALL GAUSSIAN_VECT_ROT(CROT(N),SROT(N))
         IF(LMAP) CALL GAUSSIAN_MAP_JACOB(YPTS(N), &
                                          XLON(N),XLAT(N),YLON(N),YLAT(N))
         IF(LAREA) CALL GAUSSIAN_GRID_AREA(YPTS(N),AREA(N))
       ELSE
         XPTS(N)=FILL
         YPTS(N)=FILL
       ENDIF
     ENDIF
   ENDDO
!$OMP END PARALLEL DO
 ENDIF
 DEALLOCATE(ALAT, BLAT)
 IF (ALLOCATED(YLAT_ROW)) DEALLOCATE(YLAT_ROW)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE GDSWZD_GAUSSIAN
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 SUBROUTINE GAUSSIAN_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  GAUSSIAN_ERROR   ERROR HANDLER
!   PRGMMR: GAYNO       ORG: W/NMC23       DATE: 2015-07-13
!
! ABSTRACT: UPON AN ERROR, THIS SUBPROGRAM ASSIGNS
!           A "FILL" VALUE TO THE OUTPUT FIELDS.
!
! PROGRAM HISTORY LOG:
! 2015-07-13  GAYNO     INITIAL VERSION
!
! USAGE:    CALL GAUSSIAN_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
!
!   INPUT ARGUMENT LIST:
!     IOPT     - INTEGER OPTION FLAG
!                (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
!                (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
!     NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
!     FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
!                (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
!   OUTPUT ARGUMENT LIST:
!     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
!     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
!     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
!     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE
!
 INTEGER, INTENT(IN   ) :: IOPT, NPTS
!
 REAL,    INTENT(IN   ) :: FILL
 REAL,    INTENT(  OUT) :: RLAT(NPTS),RLON(NPTS)
 REAL,    INTENT(  OUT) :: XPTS(NPTS),YPTS(NPTS)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 IF(IOPT>=0) THEN
   RLON=FILL
   RLAT=FILL
 ENDIF
 IF(IOPT<=0) THEN
   XPTS=FILL
   YPTS=FILL
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE GAUSSIAN_ERROR
!
 SUBROUTINE GAUSSIAN_VECT_ROT(CROT,SROT)
!
! SUBPROGRAM:  GAUSSIAN_VECT_ROT   VECTOR ROTATION FIELDS FOR
!                                  GAUSSIAN CYLINDRICAL GRIDS
!
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: THIS SUBPROGRAM COMPUTES THE VECTOR ROTATION SINES AND
!           COSINES FOR A GAUSSIAN CYLINDRICAL GRID.
!
! PROGRAM HISTORY LOG:
! 2015-01-21  GAYNO    INITIAL VERSION
! 2015-09-17  GAYNO    RENAME AS "GAUSSIAN_VECT_ROT"
!
! USAGE:    CALL GAUSSIAN_VECT_ROT(CROT,SROT)
!
!   INPUT ARGUMENT LIST:
!     NONE
!
!   OUTPUT ARGUMENT LIST:
!     CROT     - CLOCKWISE VECTOR ROTATION COSINES (REAL)
!     SROT     - CLOCKWISE VECTOR ROTATION SINES (REAL)
!                (UGRID=CROT*UEARTH-SROT*VEARTH;
!                 VGRID=SROT*UEARTH+CROT*VEARTH)
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
!
 IMPLICIT NONE

 REAL,                INTENT(  OUT) :: CROT, SROT

 CROT=1.0
 SROT=0.0

 END SUBROUTINE GAUSSIAN_VECT_ROT
!
 SUBROUTINE GAUSSIAN_MAP_JACOB(YPTS, XLON, XLAT, YLON, YLAT)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  GAUSSIAN_MAP_JACOB  MAP JACOBIANS FOR
!                                  GAUSSIAN CYLINDRICAL GRIDS
!
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: THIS SUBPROGRAM COMPUTES THE MAP JACOBIANS FOR
!           A GAUSSIAN CYLINDRICAL GRID.
!
! PROGRAM HISTORY LOG:
! 2015-01-21  GAYNO    INITIAL VERSION
! 2015-09-17  GAYNO    RENAME AS "GAUSSIAN_MAP_JACOB"
!
! USAGE:  CALL GAUSSIAN_MAP_JACOB(YPTS,XLON,XLAT,YLON,YLAT)
!
!   INPUT ARGUMENT LIST:
!     YPTS     - Y-INDEX OF GRID POINT (REAL)
!
!   OUTPUT ARGUMENT LIST:
!     XLON     - DX/DLON IN 1/DEGREES (REAL)
!     XLAT     - DX/DLAT IN 1/DEGREES (REAL)
!     YLON     - DY/DLON IN 1/DEGREES (REAL)
!     YLAT     - DY/DLAT IN 1/DEGREES (REAL)
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$

 IMPLICIT NONE

 REAL,                INTENT(IN   ) :: YPTS
 REAL,                INTENT(  OUT) :: XLON, XLAT, YLON, YLAT

 XLON=1/DLON
 XLAT=0.
 YLON=0.
 YLAT=YLAT_ROW(NINT(YPTS))

 END SUBROUTINE GAUSSIAN_MAP_JACOB
!
 SUBROUTINE GAUSSIAN_GRID_AREA(YPTS,AREA)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  GAUSSIAN_GRID_AREA  GRID BOX AREA FOR
!                                  GAUSSIAN CYLINDRICAL GRIDS
!
!   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
!
! ABSTRACT: THIS SUBPROGRAM COMPUTES THE GRID BOX AREA FOR
!           A GAUSSIAN CYLINDRICAL GRID.
!
! PROGRAM HISTORY LOG:
! 2015-01-21  GAYNO    INITIAL VERSION
! 2015-09-17  GAYNO    RENAME AS "GAUSSIAN_GRID_AREA"
!
! USAGE:  CALL GAUSSIAN_GRID_AREA(YPTS,AREA)
!
!   INPUT ARGUMENT LIST:
!     YPTS     - Y-INDEX OF GRID POINT (REAL)
!
!   OUTPUT ARGUMENT LIST:
!     AREA     - AREA WEIGHTS IN M**2 (REAL)
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE

 REAL,            INTENT(IN   ) :: YPTS
 REAL,            INTENT(  OUT) :: AREA

 INTEGER                        :: J

 REAL                           :: WB, WLAT, WLATA, WLATB

 J = YPTS
 WB=YPTS-J
 WLATA=BLAT(J1+JH*(J-1))
 WLATB=BLAT(J1+JH*J)
 WLAT=WLATA+WB*(WLATB-WLATA)
 AREA=RERTH**2*WLAT*DLON/DPR

 END SUBROUTINE GAUSSIAN_GRID_AREA

 END MODULE GDSWZD_GAUSSIAN_MOD

max_stars_repo_name: noaa-ocs-modeling/PaHM

max_issues_repo_name: noaa-ocs-modeling/PaHM


text:
! (C) copyright 2020 Met Office
!
! this software is licensed under the terms of the apache licence version 2.0
! which can be obtained at http://www.apache.org/licenses/license-2.0.

!> Fortran module constants used throughout the rttovonedvarcheck filter

module ufo_rttovonedvarcheck_constants_mod

use kinds
use rttov_const, only: surftype_land, surftype_sea, surftype_seaice

implicit none
private

!-----------------------------------------------------------------------------
! 1. miscellaneous definitions
!-----------------------------------------------------------------------------
integer, parameter, public :: max_string = 200 !< maximum string length

! RTTOV values for surface land, sea and seaice
integer, parameter, public :: &
  RTLand = surftype_land, & !< integer for land surface type
  RTSea  = surftype_sea, &  !< integer for sea surface type
  RTIce  = surftype_seaice  !< integer for seaice surface type

!-----------------------------------------------------------------------------
! 2. Physical Constants
!-----------------------------------------------------------------------------

real(kind_real), parameter, public :: &
  MaxTemperature = 340.0_kind_real,   & !< Maximum temperature ( K )
  MinTemperature =  70.0_kind_real      !< Minimum temperature ( K )

end module ufo_rttovonedvarcheck_constants_mod

max_stars_repo_name: NOAA-EMC/ufo

max_issues_repo_name: NOAA-EMC/ufo

max_forks_repo_name: NOAA-EMC/ufo