!
!     ------------------------------------------------------------
!     I N I T M 2
!     ------------------------------------------------------------
!
!     Determine the Rydberg constant, RMASS constant, and
!     initialize the radial grid
!
      SUBROUTINE INITM2(IPROB)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE INOUT_C
      USE PARAM_C
      USE RADIAL_C
      USE RYDBERG_C
      USE COR_C
!...Translated by Pacific-Sierra Research 77to90  4.3E  15:58:56  11/18/01
!...Switches:
      IMPLICIT NONE
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: IPROB
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: NWD = 128
      INTEGER, PARAMETER :: NELEMENT = 50
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: J
      REAL , DIMENSION(NELEMENT) :: MN
      REAL, DIMENSION(NELEMENT - 9) :: FZ
      REAL(DOUBLE) :: ME, RY
      CHARACTER :: YES
!-----------------------------------------------
!
!
      DATA ME/ 548.579903D-6/
      DATA RY/ 109737.31534/
      DATA MN/ 1, 4, 7, 9, 11, 12, 14, 16, 19, 20, 23, 24, 27, 28, 31, 32, 35, &
         40, 39, 40, 45, 48, 51, 52, 55, 56, 59, 58, 63, 64, 69, 74, 75, 80, 79&
         , 84, 85, 88, 89, 90, 93, 98, 99, 102, 103, 108, 107, 114, 115, 120/
!
      DATA FZ/ 163.5, 198.4, 239.6, 283.9, 332.7, 386.1, 444.2, 507.8, 576.4, &
         650.7, 731.2, 817.9, 911.2, 1011.4, 1119.3, 1234.9, 1359.3, 1491.9, &
         1634.4, 1786.6, 1948.7, 2122.7, 2308.2, 2507.4, 2718.5, 2944.3, 3187.5&
         , 3442.2, 3720.3, 4015.4, 4326.6, 4664.8, 5022.4, 5405.9, 5811.9, &
         6245.0, 6706.8, 7201.3, 7728.7, 8296.5, 8889.3/
!

!
! ****   DETERMINE THE RYDBERG CONSTANT AND F(Z) PARAMETER
!
      WRITE (ISCW, '(/1X,A)') ' Default Rydberg constant (y/n)'
!      READ(iread,'(A1)') YES
      YES = 'y'
      IF (YES=='Y' .OR. YES=='y' .AND. INT(Z)<=NELEMENT) THEN
         ZMU = MN(INT(Z))
      ELSE
         WRITE (ISCW, '(1X,A)') 'Enter the mass of the atom'
!         READ(iread,*) ZMU
         ZMU = 200
      ENDIF
      RYDBRG = RY/(1. + ME/ZMU)
      WRITE (6, *) ' Rydberg constant used (cm-1) = ', RYDBRG
      WRITE (6, *)
      IF (IPROB /= 1) THEN
         IF (INT(Z) < 10) THEN
            FZE = 1.566432
         ELSE IF (INT(Z) <= NELEMENT) THEN
            FZE = FZ(INT(Z-9))
         ELSE
            WRITE (ISCW, '(1X,A)') 'Enter the f(Z) value in MHz'
            READ (IREAD, *) FZE
         ENDIF
      ENDIF
      RMASS = ME/ZMU
!
! ****  INITIALIZE THE RADIAL GRID
!
      RHO = -4.0
      DO J = 1, NO
         R(J) = DEXP(RHO)/Z
         RR(J) = R(J)*R(J)
         R2(J) = DSQRT(R(J))
         RHO = RHO + H
      END DO
      RETURN
      END SUBROUTINE INITM2


text:
module test_aomp

implicit none

integer, parameter :: rstd = 8
integer :: nsize
REAL(rstd), allocatable :: a_dev(:,:,:), b_dev(:,:,:), c_dev(:,:,:)
integer :: def_async_q = 0

contains
    subroutine dec_val_dev()

    end subroutine dec_val_dev

    subroutine _compute_dev()
        integer i,j,k,l,m
#define PRIVATE_ARRAY
#ifdef PRIVATE_ARRAY
        real(rstd), dimension(nsize,nsize,nsize) :: tmp
#else
        real(rstd) :: tmp
#endif

!$omp target teams
        do i=1,nsize
!$omp parallel do simd private(tmp)
            do j=1,nsize
                do k=1,nsize
#ifdef PRIVATE_ARRAY
                    tmp(i,j,k) = b_dev(i,j,k) * c_dev(i,j,k) * i * nsize*nsize + j * nsize + k
                    if (tmp(i,j,k) > 10) then
                        a_dev(i,j,k) = tmp(i,j,k)
                    else
                        a_dev(i,j,k) = 0
                    end if
#else
                    tmp = b_dev(i,j,k) * c_dev(i,j,k) * i * nsize*nsize + j * nsize + k
                    if (tmp > 10) then
                        a_dev(i,j,k) = tmp
                    else
                        a_dev(i,j,k) = 0
                    end if
#endif
                end do
            end do
        end do
!$omp end target teams
    end subroutine _compute_dev

    subroutine compute_dev()
!$omp target update to(b_dev,c_dev,nsize)
        CALL _compute_dev()
!$omp target update from(a_dev)
    end subroutine compute_dev

    subroutine init(n_size)
        integer n_size
        integer i,j,k

        nsize = n_size
        allocate(a_dev(nsize,nsize,nsize), b_dev(nsize,nsize,nsize), c_dev(nsize,nsize,nsize))
        do i=1,nsize
            do j=1,nsize
                do k=1,nsize
                    b_dev(i,j,k) = 1
                    c_dev(i,j,k) = 2
                end do
            end do
        end do
    end subroutine init

    subroutine deinit()
        deallocate(a_dev, b_dev, c_dev)
    end subroutine deinit

end module test_aomp

program test
   use test_aomp
   call init(100)
   call compute_dev()
   call deinit()

end program test

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text: Andrew Cholerton is a student at UC Davis, double majoring in Film Studies and Technocultural Studies and also double minoring in Theatre and Psychology.  He was studying in Santander, Spain but is now back in Davis.  He owns his own business called http://www.iepfilms.com/ Immortal Eclipse Productions.  Please feel free to visit them at http://www.iepfilms.com IEPfilms.com.  He has been working with film for six years.  He currently works as a videographer, photographer, and DJ.  In addition, he works for student housing television and other such groups.  He is also the current peer advisor for Film Studies and the president of Filmmakers Ambitions.  Furthermore, he is currently producing a short zombie film, documentary, and tutorial.  Please take time to check it out, and order a dvd of http://www.lostinzombieville.com Lost in Zombieville.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


Module charge_grid_module

  Use, Intrinsic :: iso_fortran_env, Only :  wp => real64

  Use lattice_module, Only : lattice

  Implicit None

  Public :: charge_grid_calculate, charge_grid_find_range, charge_grid_forces

  Private

  Real( wp ), Parameter, Private :: pi = 3.141592653589793238462643383279502884197_wp

Contains

  Subroutine charge_grid_calculate( l, alpha, q, r, range_gauss, q_grid, error )

    Use, Intrinsic :: iso_fortran_env, Only :  wp => real64

    !$ Use omp_lib

    Use lattice_module, Only : lattice

    Implicit none

    Type( lattice )                    , Intent( In    ) :: l
    Real( wp ),                          Intent( In    ) :: alpha
    Real( wp ), Dimension( 1: ),         Intent( In    ) :: q
    Real( wp ), Dimension( 1:, 1: ),     Intent( In    ) :: r
    Integer   , Dimension( 1:3        ), Intent( In    ) :: range_gauss
    Real( wp ), Dimension( 0:, 0:, 0: ), Intent(   Out ) :: q_grid
    Integer                            , Intent(   Out ) :: error

    Real( wp ), Parameter :: pi = 3.141592653589793238462643383279502884197_wp

    ! If set adjust the charge so it adds to zero. If we have zero charge there is
    ! gauranteed to be a solution of the linear equations, as the RHS vector
    ! has no component of the null space of the LHS matrix
    ! If we have to do this this really tells us that our representation of the charge density
    ! is not good enough, so print a warning
    Logical   , Parameter :: stabilise_q     = .True.

    Real( wp ), Parameter :: stabilise_q_tol = 1e-10_wp

    Real( wp ), Dimension( :, :, :, : ), Allocatable :: q_grid_red_hack

    Real( wp ), Dimension( 1:3 ) :: ri
    Real( wp ), Dimension( 1:3 ) :: fi
    Real( wp ), Dimension( 1:3 ) :: f_point
    Real( wp ), Dimension( 1:3 ) :: r_point
    Real( wp ), Dimension( 1:3 ) :: grid_vec

    Real( wp ) :: q_norm
    Real( wp ) :: qi_norm
    Real( wp ) :: q_val
    Real( wp ) :: q_tot, q_av, c, y, t

    Integer, Dimension( 1:3 ) :: n_grid
    Integer, Dimension( 1:3 ) :: i_atom_centre
    Integer, Dimension( 1:3 ) :: i_atom_grid
    Integer, Dimension( 1:3 ) :: i_point
    Integer, Dimension( 1:3 ) :: i_grid

    Integer :: n
    Integer :: n_th, iam
    Integer :: i1, i2, i3
    Integer :: i, i_th

    error = 0

    n = Size( q )
    n_grid = Ubound( q_grid ) + 1

    q_norm = ( ( ( alpha * alpha ) / pi ) ** 1.5_wp )

    ! What follows is a complete HACK to avoid gfortran stupidly putting
    ! arrays it is reducing on the stack and hence seg faulting once the
    ! array gets beyond a certain size - beyond belief, going to whinge about this one!
    !
    ! The proper fix is to //ise the loops over the grid points assoicated with a given
    ! atom. This is a bit complicated currently due to the Modulo meaning if we have
    ! a wide gaussian we can get a race condition. But we have to fix exactly this sort of
    ! thing when going to MPI, hence fix this properly then. At the moment
    ! just want correctness, hence simple HACK as grids really aren't that large

    n_th = 1
    ! Yes, omp_max_threads is a HACK as well ....
    !$ n_th = omp_get_max_threads()
    Allocate( q_grid_red_hack( Lbound( q_grid, Dim = 1 ):Ubound( q_grid, Dim = 1 ), &
         Lbound( q_grid, Dim = 2 ):Ubound( q_grid, Dim = 2 ), &
         Lbound( q_grid, Dim = 3 ):Ubound( q_grid, Dim = 3 ), 0:n_th - 1 ) )

    !$omp parallel default( none ) shared( n, l, alpha, r, q, q_norm, n_grid, q_grid, range_gauss, q_grid_red_hack, n_th ) &
    !$omp                          private( i, i1, i2, i3, qi_norm, ri, fi, i_atom_centre,   &
    !$omp                                   i_atom_grid, i_point, f_point, r_point, grid_vec, q_val, i_grid, iam, i_th )
    !$omp do collapse( 3 )
    Do i3 = 0, n_grid( 3 ) - 1
       Do i2 = 0, n_grid( 2 ) - 1
          Do i1 = 0, n_grid( 1 ) - 1
             q_grid( i1, i2, i3 ) = 0.0_wp
          End Do
       End Do
    End Do
    iam = 0
    !$ iam = omp_get_thread_num()
    Do i3 = 0, n_grid( 3 ) - 1
       Do i2 = 0, n_grid( 2 ) - 1
          Do i1 = 0, n_grid( 1 ) - 1
             q_grid_red_hack( i1, i2, i3, iam ) = 0.0_wp
          End Do
       End Do
    End Do
    ! No need to sync at this point as in the next loop nest we only write to the
    ! bit we have just initalised
    ! Loop over atoms
    !$omp do
    Do i = 1, n
       ! Loop over points associated with atoms
       ! Find point nearest to the atom, and call this the centre for the atom grid
       ! Assumes atom in fractional 0 < ri < 1
       ri = r( :, i )
       qi_norm = q_norm * q( i )
       Call l%to_fractional( ri, fi )
       i_atom_centre = Nint( fi * n_grid )
       Do i3 = - range_gauss( 3 ), range_gauss( 3 )
          Do i2 = - range_gauss( 2 ), range_gauss( 2 )
             Do i1 = - range_gauss( 1 ), range_gauss( 1 )
                i_atom_grid = [ i1, i2, i3 ]
                ! The indices of the point in space
                i_point = i_atom_centre + i_atom_grid
                ! Transform to fractional coordinates
                f_point = Real( i_point, wp ) / n_grid
                ! And fractional to real
                Call l%to_direct( f_point, r_point )
                ! Calculate the contribution to the total charge at the point
                ! R_POINT due to the charge distribution I
                grid_vec = r_point - ri
                q_val = qi_norm * Exp( - alpha * alpha * Dot_product( grid_vec, grid_vec ) )
                ! Reflect grid indices into reference Cell
                i_grid = Modulo( i_point, n_grid )
                ! And add in
                q_grid_red_hack( i_grid( 1 ), i_grid( 2 ), i_grid( 3 ), iam ) = &
                     q_grid_red_hack( i_grid( 1 ), i_grid( 2 ), i_grid( 3 ), iam ) + q_val
             End Do
          End Do
       End Do
    End Do
    !$omp end do
    ! Synced here so OK to add up
    ! Now do the reduction manually - not needed once hacky way is fixed
    Do i_th = 0, n_th - 1
       !$omp do collapse( 3 )
       Do i3 = 0, n_grid( 3 ) - 1
          Do i2 = 0, n_grid( 2 ) - 1
             Do i1 = 0, n_grid( 1 ) - 1
                q_grid( i1, i2, i3 ) = q_grid( i1, i2, i3 ) + q_grid_red_hack( i1, i2, i3, i_th )
             End Do
          End Do
       End Do
       !$omp end do
    End Do
    !$omp end parallel

    ! If required carefully make sure the charge on the grid adds to zero
    If( stabilise_q ) Then
       ! Use Kahan summation as adding lots of very small values
       q_tot = 0.0_wp
       c = 0.0_wp
       Do i3 = 0, n_grid( 3 ) - 1
          Do i2 = 0, n_grid( 2 ) - 1
             Do i1 = 0, n_grid( 1 ) - 1
                y = q_grid( i1, i2, i3 ) - c
                t = q_tot + y
                c = ( t - q_tot ) - y
                q_tot = t
             End Do
          End Do
       End Do
!!$       Write( *, * ) 'Total charge before stabilisation ', q_tot
       If( Abs( q_tot ) > stabilise_q_tol ) Then
!!$          Write( *, * ) 'WARNING: Sum of charge on grid greater than tolerance, sum = ', &
!!$               q_tot, ' tol = ', stabilise_q_tol
          error = -1
       End If
       q_av = q_tot / Size( q_grid )
       Do i3 = 0, n_grid( 3 ) - 1
          Do i2 = 0, n_grid( 2 ) - 1
             Do i1 = 0, n_grid( 1 ) - 1
                q_grid( i1, i2, i3 ) = q_grid( i1, i2, i3 ) - q_av
             End Do
          End Do
       End Do
       q_tot = 0.0_wp
       c = 0.0_wp
       Do i3 = 0, n_grid( 3 ) - 1
          Do i2 = 0, n_grid( 2 ) - 1
             Do i1 = 0, n_grid( 1 ) - 1
                y = q_grid( i1, i2, i3 ) - c
                t = q_tot + y
                c = ( t - q_tot ) - y
                q_tot = t
             End Do
          End Do
       End Do
!!$       Write( *, * ) 'Total charge after  stabilisation ', q_tot
       If( Abs( q_tot ) > stabilise_q_tol ) Then
          error = -2
       End If
    End If

  End Subroutine charge_grid_calculate

  Subroutine charge_grid_forces( l, alpha, q, r, range_gauss, q_grid, pot_grid, ei, f )

    Use, Intrinsic :: iso_fortran_env, Only :  wp => real64

    Use lattice_module, Only : lattice

    Implicit none

    Type( lattice )                    , Intent( In    ) :: l
    Real( wp ),                          Intent( In    ) :: alpha
    Real( wp ), Dimension( 1: ),         Intent( In    ) :: q
    Real( wp ), Dimension( 1:, 1: ),     Intent( In    ) :: r
    Integer   , Dimension( 1:3        ), Intent( In    ) :: range_gauss
    Real( wp ), Dimension( 0:, 0:, 0: ), Intent( In    ) :: q_grid
    Real( wp ), Dimension( 0:, 0:, 0: ), Intent( In    ) :: pot_grid
    Real( wp ), Dimension( 1: )        , Intent(   Out ) :: ei
    Real( wp ), Dimension( 1:, 1: )    , Intent(   Out ) :: f

    Real( wp ), Parameter :: pi = 3.141592653589793238462643383279502884197_wp

    Real( wp ), Dimension( 1:3, 1:3 ) :: stress

    Real( wp ), Dimension( 1:3 ) :: ri
    Real( wp ), Dimension( 1:3 ) :: fi
    Real( wp ), Dimension( 1:3 ) :: f_point
    Real( wp ), Dimension( 1:3 ) :: r_point
    Real( wp ), Dimension( 1:3 ) :: grid_vec

    Real( wp ) :: q_norm
    Real( wp ) :: qi_norm
    Real( wp ) :: g_val
    Real( wp ) :: dV
    Real( wp ) :: s

    Integer, Dimension( 1:3 ) :: n_grid
    Integer, Dimension( 1:3 ) :: i_atom_centre
    Integer, Dimension( 1:3 ) :: i_atom_grid
    Integer, Dimension( 1:3 ) :: i_point
    Integer, Dimension( 1:3 ) :: i_grid

    Integer :: n
    Integer :: i1, i2, i3
    Integer :: i, i_alpha, i_beta

    n      = Size( q )
    n_grid = Ubound( q_grid ) + 1
    dV     = l%get_volume() / Product( n_grid )

    q_norm = ( ( ( alpha * alpha ) / pi ) ** 1.5_wp )

    stress = 0.0_wp

    !$omp parallel default( none ) shared( n, l, alpha, r, q, q_norm, n_grid, q_grid, pot_grid, range_gauss, dV, ei, f, stress ) &
    !$omp                          private( i, i1, i2, i3, qi_norm, ri, fi, i_atom_centre,   &
    !$omp                                   i_atom_grid, i_point, f_point, r_point, grid_vec, g_val, i_grid, s )
    ! Loop over atoms
    ! Stress small array, size invariant, so reduction won't stack smash
    !$omp do reduction( +:stress )
    Do i = 1, n
       ! Loop over points associated with atoms
       ! Find point nearest to the atom, and call this the centre for the atom grid
       ! Assumes atom in fractional 0 < ri < 1
       ri = r( :, i )
       qi_norm = q_norm * q( i ) * dV
       Call l%to_fractional( ri, fi )
       i_atom_centre = Nint( fi * n_grid )
       !
       ei(    i ) = 0.0_wp
       f ( :, i ) = 0.0_wp
       Do i3 = - range_gauss( 3 ), range_gauss( 3 )
          Do i2 = - range_gauss( 2 ), range_gauss( 2 )
             Do i1 = - range_gauss( 1 ), range_gauss( 1 )
                i_atom_grid = [ i1, i2, i3 ]
                ! The indices of the point in space
                i_point = i_atom_centre + i_atom_grid
                ! Transform to fractional coordinates
                f_point = Real( i_point, wp ) / n_grid
                ! And fractional to real
                Call l%to_direct( f_point, r_point )
                ! Vector to the point of interest from the centre of the gaussin
                grid_vec = r_point - ri
                ! Gaussian at that point times normalisation times the volume element
                g_val = qi_norm * Exp( - alpha * alpha * Dot_product( grid_vec, grid_vec ) )
                ! Reflect grid indices into reference Cell
                i_grid = Modulo( i_point, n_grid )
                ! Include the potential term
                g_val = g_val * pot_grid( i_grid( 1 ), i_grid( 2 ), i_grid( 3 ) )
                ! Add into the per particle energy and the force
                ei(    i ) = ei(    i ) + 0.5_wp *                            g_val
                f ( :, i ) = f ( :, i ) - 2.0_wp * alpha * alpha * grid_vec * g_val
                ! Stress term - TOTALLY UNTESTED AND PROBABLY INCOMPLETE
                ! Need short range term due to differentiation of coulomb operator,
                ! and SIC term
                Do i_beta = 1, 3
                   Do i_alpha = i_beta, 3
                      s = - 2.0_wp * g_val * alpha * alpha * grid_vec( i_alpha ) * grid_vec( i_beta )
                      stress( i_alpha, i_beta ) = stress( i_alpha, i_beta ) + s
                   End Do
                End Do
             End Do
          End Do
       End Do
    End Do
    !$omp end do
    !$omp end parallel

    ! Stress - Untested!
    ! Symmetrize stress tensor
    Do i_beta = 1, 2
       Do i_alpha = i_beta + 1, 3
          stress( i_beta, i_alpha ) = stress( i_alpha, i_beta )
       End Do
    End Do

  End Subroutine charge_grid_forces

  Subroutine charge_grid_find_range( l, alpha, n_grid, range_gauss )

    Use, Intrinsic :: iso_fortran_env, Only :  wp => real64, li => int64

    Use lattice_module, Only : lattice

    Implicit None

    Type( lattice )             , Intent( In    ) :: l
    Real( wp ),                   Intent( In    ) :: alpha
    Integer   , Dimension( 1:3 ), Intent( In    ) :: n_grid
    Integer   , Dimension( 1:3 ), Intent(   Out ) :: range_gauss

    Real( wp ), Dimension( 1:3 ) :: G
    Real( wp ), Dimension( 1:3 ) :: grid_vec

    Real( wp ) :: q_norm
    Real( wp ) :: q_val

    Integer, Dimension( 1:3 ) :: n_vec

    Integer :: i_dir

    n_vec = [ 1, 0, 0 ] ! Will use cshift to move onto next vector
    q_norm = ( ( ( alpha * alpha ) / pi ) ** 1.5_wp )
    Do i_dir = 1, 3
       Call l%get_dir_vec( n_vec, G )
       range_gauss( i_dir ) = 1
       Do
          grid_vec = ( range_gauss( i_dir ) * G ) / ( n_grid( i_dir ) )
          q_val = q_norm * Exp( - alpha * alpha * Dot_product( grid_vec, grid_vec ) )
          If( q_val < 1e-15_wp ) Exit
          range_gauss( i_dir ) = range_gauss( i_dir ) + 1
       End Do
       n_vec = Cshift( n_vec, -1 )
    End Do

  End Subroutine charge_grid_find_range

End Module charge_grid_module


text:
c-----------------------------------------------------------------------
      subroutine rw_param(nelt,nelv,ndim,nelt_all,nelv_all,nfld,ifile)
c
c     Read .rea parameters etc. and write out result iff ifile=1
c

      character*80 string
      character*1  temp

      if (ifile.eq.1) then

         dnps  = 0
         call scanout(string,'PARAMETERS FOLLOW',17,10,11)
         call lineout(11,string,80)
         read(string,*) nparam
         do i=1,nparam
            call blank(string,80)
            read(10,'(a80)') string
            write(6,*) string
            if(i.eq.23) read(string,*) dnps
            call lineout(11,string,80)
         enddo
         nps = dnps

         call scanout(string,'LOGICAL',7,10,11)  ! output to 11
         call lineout(11,string,80)
         read(string,*) nlogic                    !number of logicals
         do i = 1,nlogic
            call blank(string,80)
            read(10,'(a80)') string
            if (indx1(string,'IFHEAT',6).ne.0) then
                read(string,*)temp
                if(temp.eq.'T'.or.temp.eq.'t') nfld = 2
            endif
            call lineout(11,string,80)
         enddo

         call scanout(string,'MESH DATA',9,10,11)  ! output to 11
         call lineout(11,string,80)

         nfld = nfld  + nps
      else
         call scanout(string,'MESH DATA',9,10, 0)  ! 0 = no output
      endif

      read (10,*)  nelt, ndim, nelv

      nelt_all = nelt_all + nelt
      nelv_all = nelv_all + nelv

      return
      end
c-----------------------------------------------------------------------
      subroutine rd_xyz (x,y,z,nel,ndim)
c
      real x(8,nel),y(8,nel),z(8,nel)
      integer e
      character*80 string

      if (ndim.eq.3) then
         do e=1,nel
            read (10,80) string
            read (10,*)   (x(k,e),k=1,4)
            read (10,*)   (y(k,e),k=1,4)
            read (10,*)   (z(k,e),k=1,4)
            read (10,*)   (x(k,e),k=5,8)
            read (10,*)   (y(k,e),k=5,8)
            read (10,*)   (z(k,e),k=5,8)
         enddo
      else
         do e=1,nel
            read (10,80) string
            read (10,*)   (x(k,e),k=1,4)
            read (10,*)   (y(k,e),k=1,4)
         enddo
      endif
   80 format(a80)

      return
      end
c-----------------------------------------------------------------------
      subroutine rd_curve(ncurve,ccurve,curve,nel,ndim)

      character*1 ccurve  (12,1),cc
      real         curve(6,12,1)
      integer e,f

      real buf(6)

      buf(6) = 0

      call blank(ccurve,12*nel)
      call rzero(curve ,72*nel)

      read(10,*)
      read(10,*) ncurve
      if (ncurve.gt.0) then
         do icurve=1,ncurve
            if (nel.lt.1000) then
               read(10,60) f,e,(buf(k),k=1,5),cc
            elseif (nel.lt.1 000 000) then
               read(10,61) f,e,(buf(k),k=1,5),cc
            else
               read(10,62) f,e,(buf(k),k=1,5),cc
            endif
            ccurve(f,e) = cc
            call copy(curve(1,f,e),buf,6)
         enddo
   60    format(i3,i3,5g14.6,1x,a1)
   61    format(i2,i6,5g14.6,1x,a1)
   62    format(i2,i12,5g14.6,1x,a1)
      endif

      call cleanr(curve ,72*nel)  ! clean up small zeros

      return
      end
c-----------------------------------------------------------------------
      subroutine rd_bdry(cbc,bc,ibc,string,nel,nelo,ndim,nfld,lelt)

      character*3  cbc (6,lelt,nfld) ! note: lelt reqd
      real         bc(5,6,lelt,nfld) ! for memory reasons
      character*80 string
      integer e,f
      integer ibc (6,lelt,nfld)

      real*8 bc8(5)

      call blank(string,80)
      read (10,80) string ! ***** BOUNDARY CONDITIONS *****

      nface = 2*ndim
      do j = 1,nfld
         read (10,80) string
         if (indx1(string,'BOUN',4).ne.0) then ! we might have bcs
            if (indx1(string,'NO ',3).eq.0)  then ! we have bcs, read and count
               do e=1,nel
               do f=1,nface

                  if (nel.lt.1 000) then
                     read(10,20) cbc(f,e,j),(bc(k,f,e,j),k=1,5)
                  elseif (nel.lt.1 000 000) then
                     read(10,21) cbc(f,e,j),(bc(k,f,e,j),k=1,5)
                  else
                     read(10,22) cbc(f,e,j),(bc8(k),k=1,5)
                     ibc(f,e,j) = bc8(1)
                     do ii = 1,5
                        bc(k,f,e,j) = bc8(k)
                     enddo
                  endif
   20             format(1x,a3,6x,5g14.7)
   21             format(1x,a3,6x,5g14.7)
   22             format(1x,a3,12x,5g18.11)

c                 Update E-E and P-P pointers:
c
c                    .Assumes that P-P pointer boundaries are not
c                     eliminated during mesh merge
c
                  if (cbc(f,e,j).eq.'P  '.or.cbc(f,e,j).eq.'E  ')
     $               bc(1,f,e,j) = bc(1,f,e,j) + nelo
                  if (cbc(f,e,j).eq.'P  '.or.cbc(f,e,j).eq.'E  '.and.
     $                 nel.ge.1000000)
     $               ibc(f,e,j) = ibc(f,e,j) + nelo

               enddo
               enddo

               call cleanr(bc(1,1,1,j),30*nel)  ! clean up small zeros

            endif
         else
            goto 51
         endif

      enddo
   51 continue
   80 format(a80)

      return
      end
c-----------------------------------------------------------------------

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


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

C=======================================================================
      FUNCTION PLTPGZ(N,X,Y,Z,XQ,YQ)
      DIMENSION X(3),Y(3),Z(3),D(2)

      X21 = X(2) - X(1)
      X31 = X(3) - X(1)
      Y21 = Y(2) - Y(1)
      Y31 = Y(3) - Y(1)
      Z21 = Z(2) - Z(1)
      Z31 = Z(3) - Z(1)
      A = Y31*Z21 - Z31*Y21
      B = Z31*X21 - X31*Z21
      C = X31*Y21 - Y31*X21
      IF (C.EQ.0.) THEN
         DO 2200 I = 1,N - 2
            D(1) = X(I+1) - X(I)
            D(2) = Y(I+1) - Y(I)
            DDD = D(1)*D(1) + D(2)*D(2)
            IF (DDD.NE.0.) THEN
               GO TO 2210

            END IF

 2200    CONTINUE
 2210    CONTINUE
         IF (DDD.EQ.0.) THEN
            PLTPGZ = 0.0
            RETURN

         END IF

         DDP = D(1)*XQ + D(2)*YQ
         DDP1 = D(1)*X(1) + D(2)*Y(1)
         ALPHA = (DDP-DDP1)/DDD
         PLTPGZ = Z(1) + ALPHA* (Z(2)-Z(1))
         RETURN

      END IF

      PLTPGZ = ((A* (X(1)-XQ)+B* (Y(1)-YQ))/C) + Z(1)
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


program test_string

use, intrinsic:: iso_fortran_env, only:  stderr=>error_unit

use h5fortran, only : hdf5_file

implicit none (type, external)

character(*), parameter :: path='test_string.h5'

call test_rw(path)
print *, "OK: string read/write"
call test_overwrite(path)
print *, "OK: string overwrite"

print *,'PASSED: HDF5 string write/read'

contains

subroutine test_rw(fn)

character(*), intent(in) :: fn

type(hdf5_file) :: h
character(2) :: value
character(80) :: val1k

call h%open(fn, action='w')
call h%write('/little', '42')
call h%write('/MySentence', 'this is a little sentence.')
call h%close()

call h%open(fn, action='r')
call h%read('/little', value)

if(len_trim(value) /= 2) error stop "test_string: read length /= 2"
if (value /= '42') error stop 'test_string:  read/write verification failure. Value: '// value

call h%read('/little', val1k)
if (len_trim(val1k) /= 2) then
  write(stderr, '(a,i0,/,a)') 'expected character len_trim 2 but got len_trim() = ', len_trim(val1k), val1k
  error stop
endif

call h%close()

end subroutine test_rw

subroutine test_overwrite(fn)

character(*), intent(in) :: fn

type(hdf5_file) :: h
character(2) :: v

call h%open(fn, action='rw')
call h%write('/little', '73')
call h%close()

call h%open(fn, action='r')
call h%read('/little', v)
call h%close()

if (v /= '73') error stop 'test_string:  overwrite string failure. Value: '// v // " /= 73"

end subroutine test_overwrite

end program


! { dg-do compile }
subroutine foo1
   implicit none
   integer i
   open(1, iomsg=666)      ! { dg-error "IOMSG must be" }
   open(1, iomsg='sgk')    ! { dg-error "IOMSG must be" }
   open(1, iomsg=i)        ! { dg-error "IOMSG must be" }
   close(1, iomsg=666)     ! { dg-error "IOMSG must be" }
   close(1, iomsg='sgk')   ! { dg-error "IOMSG must be" }
   close(1, iomsg=i)       ! { dg-error "IOMSG must be" }
end subroutine foo1

subroutine foo
   implicit none
   integer i
   real :: x = 1
   write(1, *, iomsg='sgk') x   ! { dg-error "IOMSG must be" }
   write(1, *, iomsg=i)     x   ! { dg-error "IOMSG must be" }
   read(1,  *, iomsg='sgk') x   ! { dg-error "IOMSG must be" }
   read(1,  *, iomsg=i)     x   ! { dg-error "IOMSG must be" }
   flush(1,    iomsg='sgk')     ! { dg-error "IOMSG must be" }
   flush(1,    iomsg=i)         ! { dg-error "IOMSG must be" }
   rewind(1,   iomsg='sgk')     ! { dg-error "IOMSG must be" }
   rewind(1,   iomsg=i)         ! { dg-error "IOMSG must be" }
   backspace(1,iomsg='sgk')     ! { dg-error "IOMSG must be" }
   backspace(1,iomsg=i)         ! { dg-error "IOMSG must be" }
   wait(1,     iomsg='sgk')     ! { dg-error "IOMSG must be" }
   wait(1,     iomsg=i)         ! { dg-error "IOMSG must be" }
end subroutine foo

subroutine bar
   implicit none
   integer i
   real :: x = 1
   character(len=20) s(2)
   open(1, iomsg=s)         ! { dg-error "must be scalar" }
   close(1, iomsg=s)        ! { dg-error "must be scalar" }
   write(1, *, iomsg=s) x   ! { dg-error "must be scalar" }
   read(1,  *, iomsg=s) x   ! { dg-error "must be scalar" }
   flush(1,    iomsg=s)     ! { dg-error "must be scalar" }
   rewind(1,   iomsg=s)     ! { dg-error "must be scalar" }
   backspace(1,iomsg=s)     ! { dg-error "must be scalar" }
   wait(1,     iomsg=s)     ! { dg-error "must be scalar" }
end subroutine bar


text:
      subroutine run_iiiii(j,i1,i2,i3,i4,i5,DetGr,Xtwiddle0,Gtwiddle,
     . Shat6,N0)
      implicit none
C---Fixes Diiiii according to extension of Denner-Dittmaier
c---knowing D00iiii with a correction of order Delta*Diiiiii
      include 'pvDnames.f'
      include 'pvDv.f'
      include 'Darraydef.f'
      include 'Darrays.f'
      integer ep,N0,j,i1,i2,i3,i4,i5,n
      double precision DetGr,Xtwiddle0(3),Gtwiddle(3,3)
      double complex Shat6s(3,z5max,-2:0),Shat6(3,z5max,-2:0)

      do ep=-2,0
      do n=1,3
      Shat6s(n,z5(i1,i2,i3,i4,i5),ep)=Shat6(n,z5(i1,i2,i3,i4,i5),ep)
     .   -2d0*(delta(n,i1)*Dv(N0+dzziiii(z4(i2,i3,i4,i5)),ep)
     .        +delta(n,i2)*Dv(N0+dzziiii(z4(i1,i3,i4,i5)),ep)
     .        +delta(n,i3)*Dv(N0+dzziiii(z4(i1,i2,i4,i5)),ep)
     .        +delta(n,i4)*Dv(N0+dzziiii(z4(i1,i2,i3,i5)),ep)
     .        +delta(n,i5)*Dv(N0+dzziiii(z4(i1,i2,i3,i4)),ep))
      enddo

      Dv(diiiii(z5(i1,i2,i3,i4,i5))+N0,ep)=-(
     . +Gtwiddle(j,1)*Shat6s(1,z5(i1,i2,i3,i4,i5),ep)
     . +Gtwiddle(j,2)*Shat6s(2,z5(i1,i2,i3,i4,i5),ep)
     . +Gtwiddle(j,3)*Shat6s(3,z5(i1,i2,i3,i4,i5),ep)
     . -DetGr*Dv(diiiiii(z6(j,i1,i2,i3,i4,i5))+N0,ep))
     . /Xtwiddle0(j)

      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
   !check maximum and miminum thickness ratio
   ! uses
   !
   !   INTEGER (kind=4), INTENT(OUT) :: ierr    !flag for error detection
   !   REAL (kind=8) min_tr,    & !Minimum Allowable Thickness ratio
   !                 max_tr       !Maximum Allowable Thickness ratio
   !
   !   CHARACTER(len=1 ) letra  !for console input
   !
   !   TYPE (ele14), POINTER :: e    !pointer to an element data
   !   TYPE (section), POINTER :: sec    !pointer to a section data

   IF( e%lb < min_tr .OR. e%lb > max_tr)THEN      !check
     WRITE(*,"(' thickness ratio limit exceeded, TH_RATIO: ',f7.3, &
           &  /,' Do you want to continue [Y/N]')")e%lb
     DO                            !pause until answer
       READ(*,'(a1)')letra         !read the answer
       SELECT CASE (letra)         !process answer
       CASE ( 'y', 'Y')         !continue
         IF( e%lb < min_tr )THEN    !minimum exceeded
           WRITE(*,"(' New value for minimum thickness ratio')")
           READ(*,*) min_tr        !read value
           sec%rprop(2) = min_tr   !Store Minimum Thickness ratio
         ELSE                      !maximum exceeded
           WRITE(*,"(' New value for maximum thickness ratio')")
           READ(*,*) max_tr        !read value
           sec%rprop(3) = max_tr   !Store Maximum Thickness ratio
         END IF
       CASE ( 'n', 'N')         !stop process
         ierr = 1                  !modify flag
       CASE DEFAULT
         CYCLE                     !answer not valid, ask again
       END SELECT
       EXIT
     END DO
   END IF

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
! a module for storing the geometric information so we don't have to pass it

module sponge_module

  use bl_types
  use multifab_module
  use ml_layout_module

  implicit none

  real(dp_t), save :: topsponge_lo_r, topsponge_hi_r, botsponge_lo_r, botsponge_hi_r

  ! the sponge_start_density should be the density below which the
  ! sponge first turns on.  Different problems may compute this in
  ! different ways (i.e. not using sponge_center_density and
  ! sponge_start_factor), so we provide this public module variable to
  ! ensure that the rest of the code always knows at what density the
  ! sponge begins.
  real(dp_t), save, public :: sponge_start_density

  private

  public :: init_sponge, make_sponge

contains

  subroutine init_sponge(rho0)

    use geometry, only: dr, r_end_coord, spherical
    use bl_constants_module
    use probin_module, only: anelastic_cutoff, prob_hi, sponge_start_factor, dm_in, prob_lo

    real(kind=dp_t), intent(in   ) :: rho0(0:)

    real (kind = dp_t) :: rloc
    real (kind = dp_t) :: r_top
    integer            :: r

    real (kind=dp_t) :: prob_lo_r

    if (spherical .eq. 1) then
       prob_lo_r = 0.d0
    else
       prob_lo_r = prob_lo(dm_in)
    end if

    r_top = prob_lo_r + dble(r_end_coord(1,1)+1) * dr(1)
    topsponge_lo_r = r_top

    ! we do things differently than what is in Source/sponge.f90
    sponge_start_density = sponge_start_factor*anelastic_cutoff

    do r=0,r_end_coord(1,1)
       rloc = prob_lo_r + (dble(r)+HALF) * dr(1)
       if (rho0(r) < sponge_start_density) then
          topsponge_lo_r = rloc
          exit
       endif
    enddo

    do r=0,r_end_coord(1,1)
       rloc = prob_lo_r + (dble(r)+HALF) * dr(1)
       if (rho0(r) < anelastic_cutoff) then
          topsponge_hi_r = rloc
          exit
       endif
    enddo

    do r=0,r_end_coord(1,1)
       rloc = prob_lo_r + (dble(r)+HALF) * dr(1)
       if (rho0(r) < 6.d7) then
          botsponge_lo_r = rloc
          exit
       endif
    enddo

    do r=0,r_end_coord(1,1)
       rloc = prob_lo_r + (dble(r)+HALF) * dr(1)
       if (rho0(r) < 5.d7) then
          botsponge_hi_r = rloc
          exit
       endif
    enddo

    if ( parallel_IOProcessor() ) write(6,1000) topsponge_lo_r, topsponge_hi_r

1000 format('inner sponge: topsponge_lo_r      , topsponge_hi_r      : ',e20.12,2x,e20.12)

  end subroutine init_sponge

  subroutine make_sponge(sponge,dx,dt,mla)

    use bl_constants_module
    use ml_cc_restriction_module, only: ml_cc_restriction

    type(multifab) , intent(inout) :: sponge(:)
    real(kind=dp_t), intent(in   ) :: dx(:,:),dt
    type(ml_layout), intent(in   ) :: mla

    ! Local variables
    real(kind=dp_t), pointer :: sp(:,:,:,:)
    integer :: i,n,ng_sp
    integer :: lo(get_dim(sponge(1))),hi(get_dim(sponge(1))),dm,nlevs

    dm = get_dim(sponge(1))
    nlevs = mla%nlevel

    ng_sp = sponge(1)%ng

    do n=1,nlevs

       do i = 1, nfabs(sponge(n))
          sp => dataptr(sponge(n), i)
          lo =  lwb(get_box(sponge(n), i))
          hi =  upb(get_box(sponge(n), i))
          select case (dm)
          case (2)
             call mk_sponge_2d(sp(:,:,1,1),ng_sp,lo,hi,dx(n,:),dt)
          case (3)
             call mk_sponge_3d(sp(:,:,:,1),ng_sp,lo,hi,dx(n,:),dt)
          end select
       end do

    end do

    ! the loop over nlevs must count backwards to make sure the finer grids are done first
    do n=nlevs,2,-1
       ! set level n-1 data to be the average of the level n data covering it
       call ml_cc_restriction(sponge(n-1),sponge(n),mla%mba%rr(n-1,:))
    end do

  end subroutine make_sponge

  subroutine mk_sponge_2d(sponge,ng_sp,lo,hi,dx,dt)

    use bl_constants_module
    use probin_module, only: prob_lo, xrb_use_bottom_sponge, sponge_min

    integer        , intent(in   ) :: lo(:),hi(:),ng_sp
    real(kind=dp_t), intent(inout) :: sponge(lo(1)-ng_sp:,lo(2)-ng_sp:)
    real(kind=dp_t), intent(in   ) :: dx(:),dt

    integer         :: j
    real(kind=dp_t) :: y

    if (xrb_use_bottom_sponge) then

       do j = lo(2),hi(2)
          y = prob_lo(2) + (dble(j)+HALF)*dx(2)
          if(y .le. botsponge_lo_r) then
             sponge(:,j) = sponge_min
          else if(y .le. botsponge_hi_r) then
             sponge(:,j) = -HALF*(ONE-sponge_min) &
                  * cos(M_PI*(y-botsponge_lo_r)/(botsponge_hi_r-botsponge_lo_r)) &
                  + HALF*(ONE+sponge_min)
          else if(y .le. topsponge_lo_r) then
             sponge(:,j) = ONE
          else if (y .le. topsponge_hi_r) then
             sponge(:,j) = HALF*(ONE-sponge_min) &
                  * cos(M_PI*(y-topsponge_lo_r)/(topsponge_hi_r-topsponge_lo_r)) &
                  + HALF*(ONE+sponge_min)
          else
             sponge(:,j) = sponge_min
          end if
       end do

    else

       do j = lo(2),hi(2)
          y = prob_lo(2) + (dble(j)+HALF)*dx(2)
          if(y .le. topsponge_lo_r) then
             sponge(:,j) = ONE
          else if (y .le. topsponge_hi_r) then
             sponge(:,j) = HALF*(ONE-sponge_min) &
                  * cos(M_PI*(y-topsponge_lo_r)/(topsponge_hi_r-topsponge_lo_r)) &
                  + HALF*(ONE+sponge_min)
          else
             sponge(:,j) = sponge_min
          end if
       end do

    end if

  end subroutine mk_sponge_2d

  subroutine mk_sponge_3d(sponge,ng_sp,lo,hi,dx,dt)

    use bl_constants_module
    use probin_module, only: prob_lo, xrb_use_bottom_sponge, sponge_min

    integer        , intent(in   ) :: lo(:),hi(:),ng_sp
    real(kind=dp_t), intent(inout) :: sponge(lo(1)-ng_sp:,lo(2)-ng_sp:,lo(3)-ng_sp:)
    real(kind=dp_t), intent(in   ) :: dx(:),dt

    integer         :: k
    real(kind=dp_t) :: z

    if (xrb_use_bottom_sponge) then

       do k = lo(3),hi(3)
          z = prob_lo(3) + (dble(k)+HALF)*dx(3)
          if(z .le. botsponge_lo_r) then
             sponge(:,:,k) = sponge_min
          else if(z .le. botsponge_hi_r) then
             sponge(:,:,k) = -HALF*(ONE-sponge_min) &
                  * cos(M_PI*(z-botsponge_lo_r)/(botsponge_hi_r-botsponge_lo_r)) &
                  + HALF*(ONE+sponge_min)
          else if(z .le. topsponge_lo_r) then
             sponge(:,:,k) = ONE
          else if (z .le. topsponge_hi_r) then
             sponge(:,:,k) = HALF*(ONE-sponge_min) &
                  * cos(M_PI*(z-topsponge_lo_r)/(topsponge_hi_r-topsponge_lo_r)) &
                  + HALF*(ONE+sponge_min)
          else
             sponge(:,:,k) = sponge_min
          end if
       end do

    else

       do k = lo(3),hi(3)
          z = prob_lo(3) + (dble(k)+HALF)*dx(3)
          if(z .le. topsponge_lo_r) then
             sponge(:,:,k) = ONE
          else if (z .le. topsponge_hi_r) then
             sponge(:,:,k) = HALF*(ONE-sponge_min) &
                  * cos(M_PI*(z-topsponge_lo_r)/(topsponge_hi_r-topsponge_lo_r)) &
                  + HALF*(ONE+sponge_min)
          else
             sponge(:,:,k) = sponge_min
          end if
       end do

    end if

   end subroutine mk_sponge_3d

end module sponge_module

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


!-----------------------------------------------------------------------
! Surface exchange coefficients
!-----------------------------------------------------------------------
subroutine SURF_EXCH(snowdepth,Tsurf,z0,CH)

use CONSTANTS, only : &
  g,                 &! Acceleration due to gravity (m/s^2)
  vkman               ! Von Karman constant

use DRIVING, only : &
  Ta,                &! Air temperature (K)
  Ua,                &! Wind speed (m/s)
  zT,                &! Temperature and humidity measurement height (m)
  zU,                &! Wind measurement height (m)
  zvar                ! Subtract snow depth from measurement height

use MODELS, only: &
  em                  ! Surface exchange model   0 - fixed
                      !                          1 - stability correction
use PARAMETERS, only : &
  bstb                ! Stability slope parameter

implicit none

real, intent(in) :: &
  snowdepth,         &! Snow depth (m)
  Tsurf,             &! Surface temperature (K)
  z0                  ! Roughness length for momentum (m)

real, intent(out) :: &
  CH                  ! Transfer coefficient for heat and moisture

real :: &
  CD,                &! Drag coefficient
  fh,                &! Stability correction
  RiB,               &! Bulk Richardson number
  z0h,               &! Roughness length for heat and moisture (m)
  zTs                 ! Adjusted temperature measurement height (m)

zTs = zT
if (zvar) then
  zTs = zT - snowdepth
  zTs = max(zTs, 1.)
end if

! Neutral exchange coefficients
z0h = 0.1*z0
CD = (vkman / log(zU/z0))**2
CH = vkman**2 / (log(zU/z0)*log(zTs/z0h))

! Stability correction (Louis et al. 1982, quoted by Beljaars 1992)
if (em == 1) then
  RiB = g*(Ta - Tsurf)*zU**2 / (zTs*Ta*Ua**2)
  if (RiB > 0) then
    fh = 1/(1 + 3*bstb*RiB*sqrt(1 + bstb*RiB))
  else
    fh = 1 - 3*bstb*RiB / (1 + 3*bstb**2*CD*sqrt(-RiB*zU/z0))
  end if
  CH = fh*CH
end if

end subroutine SURF_EXCH


subroutine anbkerror(gradx,grady)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    anbkerror apply anisotropic background error covariance
!   prgmmr: parrish          org: np22                date: 2005-02-03
!
! abstract: apply regional anisotropic background error.
!
! program history log:
!   2005-02-08  parrish
!   2005-04-29  parrish - replace coarse2fine with fgrid2agrid;
!                         remove ansmoothrf_reg_d
!   2006-11-30  todling - add fpsproj as arg to (t)balance routine(s)
!   2008-10-10  derber - add strong constraint to background error
!   2008-12-29  todling - update interface to strong_bk/bk_ad
!   2009-04-13  derber - move strong_bk into balance
!   2009-07-01  sato - update for global mode
!   2010-05-05  todling - update to use gsi_bundle
!   2010-06-22  todling - update to better handle bundle pointers
!   2010-06-29  lueken - replaced tv with t in call to gsi_bundlegetpointer
!   2010-08-19  lueken - add only to module use
!   2012-10-09  Gu - add fut2ps as arg to (t)balance routine(s)
!   2013-05-23  zhu    - add ntclen for aircraft temperature bias correction
!   2014-02-07  pondeca - update to handle motley variables. rename p_st to p_sf
!   2014-02-14  pondeca - update to handle optional separation of sf and vp control variables
!                         into land-only and water-only parts
!   2015-05-02  parrish - add subroutine ansmoothrf_reg_sub2slab_option, and
!                         parameter rtma_bkerr_sub2slab.
!                         rtma_bkerr_sub2slab = F, then use ansmoothrf_reg_subdomain_option
!                         rtma_bkerr_sub2slab = T, then use ansmoothrf_reg_sub2slab_option
!                         This allows dual resolution for the anisotropic recursive filter, which
!                         currently can only be used with full horizontal domain (slab) storage.
!   2015-07-02  pondeca - update slab mode option to work with any number of control variables
!
!   input argument list:
!     gradx    - input field
!
!   output
!     grady    - background structure * gradx
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2
  use jfunc, only: nsclen,npclen,ntclen
  use balmod, only: balance,tbalance,ke_vp,bvk
  use berror, only: varprd,fpsproj,fut2ps
  use constants, only: zero
  use control_vectors, only: control_vector,assignment(=)
  use control_vectors, only: mvars,nrf,nrf_var,nrf_3d
  use gsi_4dvar, only: nsubwin
  use timermod, only: timer_ini,timer_fnl
  use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundlemerge,gsi_bundle,gsi_bundledup,gsi_bundledestroy
  implicit none

! Declare passed variables
  type(control_vector),intent(inout) :: gradx
  type(control_vector),intent(inout) :: grady

! Declare local variables
  integer(i_kind) i,j,k,ii,istatus
  real(r_kind),dimension(:,:,:),pointer::p_t  =>NULL()
  real(r_kind),dimension(:,:,:),pointer::p_sf =>NULL()
  real(r_kind),dimension(:,:,:),pointer::p_vp =>NULL()
  real(r_kind),dimension(:,:  ),pointer::p_ps =>NULL()
  real(r_kind),dimension(:,:,:),pointer::p_sfwter =>NULL()
  real(r_kind),dimension(:,:,:),pointer::p_vpwter =>NULL()
  real(r_kind),pointer::rank2a(:,:)  =>NULL()
  real(r_kind),pointer::rank2b(:,:)  =>NULL()
  real(r_kind),pointer::rank3a(:,:,:)=>NULL()
  real(r_kind),pointer::rank3b(:,:,:)=>NULL()
  logical do_balance, do_balancewter
  integer(i_kind), parameter :: myvars = 6
  integer(i_kind) :: ipnts(myvars)
  type(gsi_bundle) :: mbundle
  character(len=6), parameter :: myvnames(myvars) = (/  &
       'sf    ', 'vp    ', 'ps    ', 't     ', 'sfwter', 'vpwter'/)

! Initialize timer
  call timer_ini('anbkerror')

! Put things in grady first since operations change input variables
  grady=gradx

! Since each internal vector [step(jj)] of grad has the same structure, pointers
! are the same independent of the subwindow jj
call gsi_bundlegetpointer (grady%step(1),myvnames,ipnts,istatus)

! Define what to do depending on what's in CV and SV
do_balance=ipnts(1)>0.and.ipnts(2)>0.and.ipnts(3)>0.and.ipnts(4)>0
do_balancewter=ipnts(5)>0.and.ipnts(6)>0

! Loop on control steps
  do ii=1,nsubwin

!    Create temporary bundle which merges grady%step(ii) with grady%motley(ii)
     if(mvars>0) then
        call gsi_bundlemerge(mbundle,grady%step(ii),grady%motley(ii),' add motley to step',istatus)
     else
        call gsi_bundledup(grady%step(ii),mbundle,' copy of step ',istatus)
     end if

!    Transpose of balance equation
     if(do_balancewter) then
        call gsi_bundlegetpointer (mbundle,'sfwter',p_sfwter,  istatus)
        call gsi_bundlegetpointer (mbundle,'vpwter',p_vpwter,  istatus)

       !Adjoint of contribution to velocity potential from streamfunction.
        do k=1,ke_vp
           do j=1,lon2
              do i=1,lat2
                 p_sfwter(i,j,k)=p_sfwter(i,j,k)+bvk(i,j,k)*p_vpwter(i,j,k)
              end do
           end do
        end do
     endif

     if(do_balance) then
        call gsi_bundlegetpointer (mbundle,'sf',p_sf,  istatus)
        call gsi_bundlegetpointer (mbundle,'vp',p_vp,  istatus)
        call gsi_bundlegetpointer (mbundle,'ps',p_ps,  istatus)
        call gsi_bundlegetpointer (mbundle,'t ',p_t,   istatus)
        call gsi_bundlegetpointer (mbundle,'sfwter',p_sfwter,  istatus)
        call gsi_bundlegetpointer (mbundle,'vpwter',p_vpwter,  istatus)
        call tbalance(p_t,p_ps,p_sf,p_vp,fpsproj,fut2ps)
     endif

!    Apply variances, as well as vertical & horizontal parts of background error
     call anbkgcov(mbundle)

!    Balance equation
     if(do_balance) call balance(p_t,p_ps,p_sf,p_vp,fpsproj,fut2ps)

     if(do_balancewter) then
       do k=1,ke_vp
          do j=1,lon2
             do i=1,lat2
                p_vpwter(i,j,k)=p_vpwter(i,j,k)+bvk(i,j,k)*p_sfwter(i,j,k)
             end do
          end do
       end do
     endif

!    Transfer step part of mbundle back to grady%step(ii)
     do i=1,nrf
        if(nrf_3d(i)) then
           call gsi_bundlegetpointer(mbundle,trim(nrf_var(i)),rank3a,istatus)
           call gsi_bundlegetpointer(grady%step(ii),trim(nrf_var(i)),rank3b,istatus)
           rank3b=rank3a
        else
           call gsi_bundlegetpointer(mbundle,trim(nrf_var(i)),rank2a,istatus)
           call gsi_bundlegetpointer(grady%step(ii),trim(nrf_var(i)),rank2b,istatus)
           rank2b=rank2a
        end if
     end do

!    Transfer motley part of mbundle back to grady%motley(ii)
     do i=1,mvars
        call gsi_bundlegetpointer(mbundle,trim(nrf_var(nrf+i)),rank2a,istatus)
        call gsi_bundlegetpointer(grady%motley(ii),trim(nrf_var(nrf+i)),rank2b,istatus)
        rank2b=rank2a
     end do
  end do

! clean work space
  call gsi_bundledestroy(mbundle,istatus)
  if(istatus/=0) then
     write(6,*) ' in anbkerror: trouble destroying work mbundle'
     call stop2(999)
  endif

! Take care of background error for bias correction terms
  do i=1,nsclen
     grady%predr(i)=grady%predr(i)*varprd(i)
  end do
  do i=1,npclen
     grady%predp(i)=grady%predp(i)*varprd(nsclen+i)
  end do
  if(ntclen>0)then
     do i=1,ntclen
        grady%predt(i)=grady%predt(i)*varprd(nsclen+npclen+i)
     end do
  end if

! Finalize timer
  call timer_fnl('anbkerror')

end subroutine anbkerror

subroutine anbkgcov(bundle)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    anbkgcov    apply anisotropic background error covar
!   prgmmr: parrish        org: np22                date: 2005-02-14
!
! abstract: apply regional anisotropic background error covariance
!
! program history log:
!   2005-02-14  parrish
!   2009-07-01  sato - update for global mode
!   2010-05-20  todling - update fit interface to sug2grid/grid2sub (bundle)
!   2010-06-22  todling - update interface (remove cwmr since it's in bunlde)
!   2010-06-29  lueken - added if(ipnts(2)>0) to second call of anbkgvar
!   2011-02-22  zhu - replace the argument list of ansmoothrf_reg_subdomain_option by a bundle
!   2012-06-25  parrish - replace sub2grid and grid2sub calls with general_sub2grid, general_grid2sub.
!                 NOTE:  This will not work with sst and the motley variables stl,sti.  However
!                        this is not currently used in this version of RTMA.
!   2014-03-19  pondeca - add wspd10m
!   2014-04-10  pondeca - add td2m,mxtm,mitm,pmsl
!   2014-05-07  pondeca - add howv
!   2016-05-03  pondeca - add uwnd10m, vwnd10m
!
!   input argument list:
!     t        - t on subdomain
!     p        - p surface pressure on subdomain
!     q        - q on subdomain
!     oz       - ozone on subdomain
!     skint    - skin temperature on subdomain
!     sst      - sea surface temperature on subdomain
!     stl      - land surface temperature on subdomain
!     sti      - ice surface temperature on subdomain
!     gust     - 10-m gust on subdomain
!     vis      - surface visibility on subdomain
!     wspd10m  - 10m-wind speed on subdomain
!     td2m     - td on subdomain
!     mxtm     - daily maximum temperature
!     mitm     - daily  minimum temperature
!     pmsl     - pressure at mean sea level
!     howv     - significant wave height
!     tcamt    - total cloud amount
!     lcbas    - lowest cloud base height
!     cldch    - cloud ceiling height
!     uwnd10m  - 10m-uwind on subdomain
!     vwnd10m  - 10m-vwind on subdomain
!     pswter   - water surface pressure on subdomain
!     twter    - water 2m-temperature on subdomain
!     qwter    - water 2m-specific humidity on subdomain
!     gustwter - water 10m-gust on subdomain
!     wspd10mwter  - water 10m-wind speed on subdomain
!     td2mwter - water td on subdomain
!     mxtmwter - water daily maximum temperature
!     mitmwter - water daily minimum temperature
!     uwnd10mwter  - water 10m-uwind on subdomain
!     vwnd10mwter  - water 10m-vwind on subdomain
!
!   output argument list:
!                 all after smoothing, combining scales
!     t        - t on subdomain
!     p        - p surface pressure on subdomain
!     q        - q on subdomain
!     oz       - ozone on subdomain
!     skint    - skin temperature on subdomain
!     sst      - sea surface temperature on subdomain
!     stl      - land surface temperature on subdomain
!     sti      - ice surface temperature on subdomain
!     gust     - 10-m gust on subdomain
!     vis      - surface visibility on subdomain
!     wspd10m  - 10m-wind speed on subdomain
!     td2m     - td on subdomain
!     mxtm     - daily maximum temperature
!     mitm     - daily minimum temperature
!     pmsl     - pressure at mean sea level
!     howv     - significant wave height
!     tcamt    - total cloud amount
!     lcbas    - lowest cloud base height
!     cldch    - cloud ceiling height
!     uwnd10m  - 10m-uwind on subdomain
!     vwnd10m  - 10m-vwind on subdomain
!     pswter   - water surface pressure on subdomain
!     twter    - water 2m-temperature on subdomain
!     qwter    - water 2m-specific humidity on subdomain
!     gustwter - water 10m-gust on subdomain
!     wspd10mwter  - water 10m-wind speed on subdomain
!     td2mwter - water td on subdomain
!     mxtmwter - water daily maximum temperature
!     mitmwter - water daily minimum temperature
!     uwnd10mwter  - water 10m-uwind on subdomain
!     vwnd10mwter  - water 10m-vwind on subdomain
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2,twodvar_regional
  use anberror, only: rtma_subdomain_option,nsmooth, nsmooth_shapiro,rtma_bkerr_sub2slab
  use constants, only: zero
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  implicit none

! Passed Variables
  type(gsi_bundle),                 intent(inout) :: bundle

! Local Variables

  integer(i_kind) n,istatus
  integer(i_kind) i_sst,i_stl,i_sti,i_ps,i_t,i_q,i_gust,i_wspd10m, &
                  i_td2m,i_mxtm,i_mitm,i_uwnd10m,i_vwnd10m, &
                  i_pswter,i_twter,i_qwter,i_gustwter,i_wspd10mwter, &
                  i_td2mwter,i_mxtmwter,i_mitmwter,i_uwnd10mwter,i_vwnd10mwter

  real(r_kind),dimension(lat2,lon2) :: skint,sst,stl,sti
  real(r_kind),dimension(lat2,lon2) :: field,fld,fldwter

  real(r_kind),dimension(:,:),   pointer :: ptrsst=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrstl=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrsti=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrps=>NULL()
  real(r_kind),dimension(:,:,:), pointer :: ptrt=>NULL()
  real(r_kind),dimension(:,:,:), pointer :: ptrq=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrgust=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrwspd10m=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrtd2m=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrmxtm=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrmitm=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptruwnd10m=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrvwnd10m=>NULL()
  real(r_kind),dimension(:,:,:), pointer :: ptr3d=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrpswter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrtwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrqwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrgustwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrwspd10mwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrtd2mwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrmxtmwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrmitmwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptruwnd10mwter=>NULL()
  real(r_kind),dimension(:,:),   pointer :: ptrvwnd10mwter=>NULL()

! Perform simple vertical smoothing while fields are in sudomain mode.
! The accompanying smoothing in the horizontal is performed inside the
! recursive filter. Motivation: Reduce possible high frequency noise in
! the analysis that would arise from the use of a "non-blending" RF algorithm.

  do n=1,bundle%n3d
     call gsi_bundlegetpointer ( bundle,bundle%r3(n)%shortname,ptr3d,istatus )
     call vert_smther(ptr3d,nsmooth,nsmooth_shapiro)
  end do

! Break up skin temp into components

! First get pointers
  !SST
  call gsi_bundlegetpointer (bundle, 'sst', i_sst,  istatus)
  call gsi_bundlegetpointer (bundle, 'stl', i_stl,  istatus)
  call gsi_bundlegetpointer (bundle, 'sti', i_sti,  istatus)
  if (i_sst>0) call gsi_bundlegetpointer (bundle, 'sst', ptrsst, istatus)
  if (i_stl>0) call gsi_bundlegetpointer (bundle, 'stl', ptrstl, istatus)
  if (i_sti>0) call gsi_bundlegetpointer (bundle, 'sti', ptrsti, istatus)

  if(twodvar_regional) then
     !SURFACE PRESSURE
     call gsi_bundlegetpointer (bundle, 'ps',     i_ps,      istatus)
     call gsi_bundlegetpointer (bundle, 'pswter', i_pswter,  istatus)
     if (i_ps>0)     call gsi_bundlegetpointer (bundle, 'ps',     ptrps,     istatus)
     if (i_pswter>0) call gsi_bundlegetpointer (bundle, 'pswter', ptrpswter, istatus)

     !TEMPERATURE
     call gsi_bundlegetpointer (bundle, 't',     i_t,      istatus)
     call gsi_bundlegetpointer (bundle, 'twter', i_twter,  istatus)
     if (i_t>0)     call gsi_bundlegetpointer (bundle, 't',     ptrt,     istatus)
     if (i_twter>0) call gsi_bundlegetpointer (bundle, 'twter', ptrtwter, istatus)

     !SPECIFIC HUMIDITY
     call gsi_bundlegetpointer (bundle, 'q',     i_q,      istatus)
     call gsi_bundlegetpointer (bundle, 'qwter', i_qwter,  istatus)
     if (i_q>0)     call gsi_bundlegetpointer (bundle, 'q',     ptrq,     istatus)
     if (i_qwter>0) call gsi_bundlegetpointer (bundle, 'qwter', ptrqwter, istatus)

     !WIND GUST
     call gsi_bundlegetpointer (bundle, 'gust',     i_gust,      istatus)
     call gsi_bundlegetpointer (bundle, 'gustwter', i_gustwter,  istatus)
     if (i_gust>0)     call gsi_bundlegetpointer (bundle, 'gust',     ptrgust,     istatus)
     if (i_gustwter>0) call gsi_bundlegetpointer (bundle, 'gustwter', ptrgustwter, istatus)

     !10-m WIND SPEED
     call gsi_bundlegetpointer (bundle, 'wspd10m',     i_wspd10m,      istatus)
     call gsi_bundlegetpointer (bundle, 'wspd10mwter', i_wspd10mwter,  istatus)
     if (i_wspd10m>0)     call gsi_bundlegetpointer (bundle, 'wspd10m',     ptrwspd10m,     istatus)
     if (i_wspd10mwter>0) call gsi_bundlegetpointer (bundle, 'wspd10mwter', ptrwspd10mwter, istatus)

     !2-m DEW POINT
     call gsi_bundlegetpointer (bundle, 'td2m',     i_td2m,      istatus)
     call gsi_bundlegetpointer (bundle, 'td2mwter', i_td2mwter,  istatus)
     if (i_td2m>0)     call gsi_bundlegetpointer (bundle, 'td2m',     ptrtd2m,     istatus)
     if (i_td2mwter>0) call gsi_bundlegetpointer (bundle, 'td2mwter', ptrtd2mwter, istatus)

     !MAXIMUM TEMPERATURE
     call gsi_bundlegetpointer (bundle, 'mxtm',     i_mxtm,      istatus)
     call gsi_bundlegetpointer (bundle, 'mxtmwter', i_mxtmwter,  istatus)
     if (i_mxtm>0)     call gsi_bundlegetpointer (bundle, 'mxtm',     ptrmxtm,     istatus)
     if (i_mxtmwter>0) call gsi_bundlegetpointer (bundle, 'mxtmwter', ptrmxtmwter, istatus)

     !MINIMUM TEMPERATURE
     call gsi_bundlegetpointer (bundle, 'mitm',     i_mitm,      istatus)
     call gsi_bundlegetpointer (bundle, 'mitmwter', i_mitmwter,  istatus)
     if (i_mitm>0)     call gsi_bundlegetpointer (bundle, 'mitm',     ptrmitm,     istatus)
     if (i_mitmwter>0) call gsi_bundlegetpointer (bundle, 'mitmwter', ptrmitmwter, istatus)

     !10-m UWIND
     call gsi_bundlegetpointer (bundle, 'uwnd10m',     i_uwnd10m,      istatus)
     call gsi_bundlegetpointer (bundle, 'uwnd10mwter', i_uwnd10mwter,  istatus)
     if (i_uwnd10m>0)     call gsi_bundlegetpointer (bundle, 'uwnd10m',     ptruwnd10m,     istatus)
     if (i_uwnd10mwter>0) call gsi_bundlegetpointer (bundle, 'uwnd10mwter', ptruwnd10mwter, istatus)

     !10-m VWIND
     call gsi_bundlegetpointer (bundle, 'vwnd10m',     i_vwnd10m,      istatus)
     call gsi_bundlegetpointer (bundle, 'vwnd10mwter', i_vwnd10mwter,  istatus)
     if (i_vwnd10m>0)     call gsi_bundlegetpointer (bundle, 'vwnd10m',     ptrvwnd10m,     istatus)
     if (i_vwnd10mwter>0) call gsi_bundlegetpointer (bundle, 'vwnd10mwter', ptrvwnd10mwter, istatus)
  endif

! Break up skin temp

  if(i_sst>0 .and. i_stl>0 .and. i_sti>0) then
    stl=zero
    sst=zero
    sti=zero
    skint=ptrsst

    call anbkgvar(skint,sst,stl,sti,0)

    ptrsst=sst
    ptrstl=stl
    ptrsti=sti
  endif

  if(twodvar_regional) then
    if(i_ps>0.and.i_pswter>0) then
      fld=zero
      fldwter=zero
      field=ptrps

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrps=fld
      ptrpswter=fldwter
    endif

    if(i_t>0.and.i_twter>0) then
      fld=zero
      fldwter=zero
      field(:,:)=ptrt(:,:,1)

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrt(:,:,1)=fld(:,:)
      ptrtwter=fldwter
    endif

    if(i_q>0.and.i_qwter>0) then
      fld=zero
      fldwter=zero
      field(:,:)=ptrq(:,:,1)

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrq(:,:,1)=fld(:,:)
      ptrqwter=fldwter
    endif

    if(i_gust>0.and.i_gustwter>0) then
      fld=zero
      fldwter=zero
      field=ptrgust

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrgust=fld
      ptrgustwter=fldwter
    endif

    if(i_wspd10m>0.and.i_wspd10mwter>0) then
      fld=zero
      fldwter=zero
      field=ptrwspd10m

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrwspd10m=fld
      ptrwspd10mwter=fldwter
    endif

    if(i_td2m>0.and.i_td2mwter>0) then
      fld=zero
      fldwter=zero
      field=ptrtd2m

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrtd2m=fld
      ptrtd2mwter=fldwter
    endif

    if(i_mxtm>0.and.i_mxtmwter>0) then
      fld=zero
      fldwter=zero
      field=ptrmxtm

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrmxtm=fld
      ptrmxtmwter=fldwter
    endif

    if(i_mitm>0.and.i_mitmwter>0) then
      fld=zero
      fldwter=zero
      field=ptrmitm

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrmitm=fld
      ptrmitmwter=fldwter
    endif

    if(i_uwnd10m>0.and.i_uwnd10mwter>0) then
      fld=zero
      fldwter=zero
      field=ptruwnd10m

      call anbkgvar_lw(field,fld,fldwter,0)

      ptruwnd10m=fld
      ptruwnd10mwter=fldwter
    endif

    if(i_vwnd10m>0.and.i_vwnd10mwter>0) then
      fld=zero
      fldwter=zero
      field=ptrvwnd10m

      call anbkgvar_lw(field,fld,fldwter,0)

      ptrvwnd10m=fld
      ptrvwnd10mwter=fldwter
    endif
  endif

! Apply auto-covariance
  if(rtma_subdomain_option) then
      if(rtma_bkerr_sub2slab) then
         call ansmoothrf_reg_sub2slab_option(bundle)
      else
         call ansmoothrf_reg_subdomain_option(bundle)
      end if
  else
     call ansmoothrf(bundle)
  end if

!Adjoint of simple vertical smoothing
  do n=bundle%n3d,1,-1
     call gsi_bundlegetpointer ( bundle,bundle%r3(n)%shortname,ptr3d,istatus )
     call tvert_smther(ptr3d,nsmooth,nsmooth_shapiro)
  end do

  if(twodvar_regional) then
    if(i_vwnd10m>0.and.i_vwnd10mwter>0) then
      fld=ptrvwnd10m
      fldwter=ptrvwnd10mwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrvwnd10m=field
!     ignore content of remaining arrays
    endif

    if(i_uwnd10m>0.and.i_uwnd10mwter>0) then
      fld=ptruwnd10m
      fldwter=ptruwnd10mwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptruwnd10m=field
!     ignore content of remaining arrays
    endif

    if(i_mitm>0.and.i_mitmwter>0) then
      fld=ptrmitm
      fldwter=ptrmitmwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrmitm=field
!     ignore content of remaining arrays
    endif

    if(i_mxtm>0.and.i_mxtmwter>0) then
      fld=ptrmxtm
      fldwter=ptrmxtmwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrmxtm=field
!     ignore content of remaining arrays
    endif

    if(i_td2m>0.and.i_td2mwter>0) then
      fld=ptrtd2m
      fldwter=ptrtd2mwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrtd2m=field
!     ignore content of remaining arrays
    endif

    if(i_wspd10m>0.and.i_wspd10mwter>0) then
      fld=ptrwspd10m
      fldwter=ptrwspd10mwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrwspd10m=field
!     ignore content of remaining arrays
    endif

    if(i_gust>0.and.i_gustwter>0) then
      fld=ptrgust
      fldwter=ptrgustwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrgust=field
!     ignore content of remaining arrays
    endif

    if(i_q>0.and.i_qwter>0) then
      fld(:,:)=ptrq(:,:,1)
      fldwter=ptrqwter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrq(:,:,1)=field(:,:)
!     ignore content of remaining arrays
    endif

    if(i_t>0.and.i_twter>0) then
      fld(:,:)=ptrt(:,:,1)
      fldwter=ptrtwter
      field(:,:)=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrt(:,:,1)=field(:,:)
!     ignore content of remaining arrays
    endif

    if(i_ps>0.and.i_pswter>0) then
      fld=ptrps
      fldwter=ptrpswter
      field=zero

      call anbkgvar_lw(field,fld,fldwter,1)

      ptrps=field
!     ignore content of remaining arrays
    endif
  endif

!==> combine sst,stl, and sti into skin temperature field
  if(i_sst>0 .and. i_stl>0 .and. i_sti>0) then
    sst=ptrsst
    stl=ptrstl
    sti=ptrsti
    skint=zero

    call anbkgvar(skint,sst,stl,sti,1)

    ptrsst=skint
!   ignore contents of remaining arrays
  endif

end subroutine anbkgcov

subroutine anbkgvar(skint,sst,stl,sti,iflg)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    anbkgvar_reg      manipulate skin temp
!   prgmmr: parrish          org: np22                date: 1990-10-06
!
! abstract: manipulate skin temp <--> sst,sfc temp, and ice temp fields
!
! program history log:
!   2005-01-22  parrish
!   2008-06-05  safford - rm unused uses
!   2012-06-25  parrish - remove _i_kind from integer constants
!
!   input argument list:
!     skint    - skin temperature grid values
!     sst      - sst grid values
!     stl      - land surface temperature grid values
!     sti      - snow/ice covered surface temperature grid values
!     iflg     - flag for skin temperature manipulation
!                0: skint --> sst,stl,sti
!                1: sst,stl,sti --> skint
!
!   output argument list:
!     skint    - skin temperature grid values
!     sst      - sst grid values
!     stl      - land surface temperature grid values
!     sti      - snow/ice covered surface temperature grid values
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2
  use guess_grids, only: isli2
  implicit none

! Declare passed variables
  integer(i_kind)                  ,intent(in   ) :: iflg
  real(r_kind),dimension(lat2,lon2),intent(inout) :: skint,sst,stl,sti

! Declare local variables
  integer(i_kind) i,j

       do j=1,lon2
          do i=1,lat2
             if(iflg == 0) then
! Break skin temperature into components
!          If land point
                if(isli2(i,j) == 1) then
                   stl(i,j)=skint(i,j)
!          If ice
                else if(isli2(i,j) == 2) then
                   sti(i,j)=skint(i,j)
!          Else treat as a water point
                else
                   sst(i,j)=skint(i,j)
                end if

             else if (iflg==1) then
! Combine sst,stl, and sti into skin temperature field
!          Land point, load land sfc t into skint
                if(isli2(i,j) == 1) then
                   skint(i,j)=stl(i,j)
!          Ice, load ice temp into skint
                else if(isli2(i,j) == 2) then
                   skint(i,j)=sti(i,j)
!          Treat as a water point, load sst into skint
                else
                   skint(i,j)=sst(i,j)
                end if
             end if
          end do
       end do

  return
end subroutine anbkgvar

subroutine anbkgvar_lw_original(field,fld,fldwter,iflg)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    anbkgvar_lw land/water field manipulation
!   prgmmr: pondeca          org: np22                date: 2013-04-02
!
! abstract: manipulate field  <--> fld,fldwter
!           (based on anbkgvar)
!
! program history log:
!   2013-04-02 pondeca
!
!   input argument list:
!     field   - field grid values
!     fld     - field grid values over land
!     fldwter - field grid values over water
!     iflg    - flag for field manipulation
!               0: field --> fld,fldwter
!               1: fld,fldwter --> field
!
!   output argument list:
!     field   - field grid values
!     fld     - field grid values over land
!     fldwter - field grid values over water
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2
  use guess_grids, only: isli2
  implicit none

! Declare passed variables
  integer(i_kind)                  ,intent(in   ) :: iflg
  real(r_kind),dimension(lat2,lon2),intent(inout) :: field, &
                                                     fld,fldwter

! Declare local variables
  integer(i_kind) i,j

       do j=1,lon2
          do i=1,lat2
             if(iflg == 0) then
! Break field into components
!          If land point
                if(isli2(i,j) == 1) then
                   fld(i,j)=field(i,j)

!          Else treat as a water point
                else
                   fldwter(i,j)=field(i,j)
                end if

             else if (iflg==1) then
! Combine fld,fldwter into field
!          Land point
                if(isli2(i,j) == 1) then
                   field(i,j)=fld(i,j)
!          Treat as a water point,
                else
                   field(i,j)=fldwter(i,j)
                end if
             end if
          end do
       end do

  return
end subroutine anbkgvar_lw_original

subroutine anbkgvar_lw(field,fld,fldwter,iflg)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    anbkgvar_lw land/water field manipulation
!   prgmmr: pondeca          org: np22                date: 2013-04-02
!
! abstract: manipulate field  <--> fld,fldwter
!           (based on anbkgvar)
!
! program history log:
!   2013-04-02 pondeca
!
!   input argument list:
!     field   - field grid values
!     fld     - field grid values over land
!     fldwter - field grid values over water
!     iflg    - flag for field manipulation
!               0: field --> fld,fldwter
!               1: fld,fldwter --> field
!
!   output argument list:
!     field   - field grid values
!     fld     - field grid values over land
!     fldwter - field grid values over water
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2,region_lat,region_lon, &
                     nlon_regional,nlat_regional,istart,jstart
  use guess_grids, only: isli2
  use mpimod, only: mype
  use constants, only: rad2deg
  implicit none

! Declare passed variables
  integer(i_kind)                  ,intent(in   ) :: iflg
  real(r_kind),dimension(lat2,lon2),intent(inout) :: field, &
                                                     fld,fldwter

! Declare local variables
  integer(i_kind) i,j
  integer(i_kind) iglob,jglob
  integer(i_kind) mm1
  real(r_kind) flon,flat
  logical glerlarea

! Declare local parameters
!    Great Lakes
  real(r_kind),parameter::flon1=-93._r_kind
  real(r_kind),parameter::flon2=-75._r_kind
  real(r_kind),parameter::flat1=40.5_r_kind
  real(r_kind),parameter::flat2=49.5_r_kind

!    Great Salt Lake
  real(r_kind),parameter::slon1=-113._r_kind
  real(r_kind),parameter::slon2=-112._r_kind
  real(r_kind),parameter::slat1=40.6_r_kind
  real(r_kind),parameter::slat2=41.7_r_kind

  mm1=mype+1

       do j=1,lon2
          jglob=max(1,min(j+jstart(mm1)-2,nlon_regional))
          do i=1,lat2
             iglob=max(1,min(i+istart(mm1)-2,nlat_regional))
             flat=region_lat(iglob,jglob)*rad2deg
             flon=region_lon(iglob,jglob)*rad2deg
             glerlarea=(flat>=flat1.and.flat<=flat2).and.(flon>=flon1.and.flon<=flon2)
             glerlarea=glerlarea.or.((flat>=slat1.and.flat<=slat2).and.(flon>=slon1.and.flon<=slon2))

             if(iflg == 0) then
! Break field into components
!          If land point
                if(isli2(i,j) == 1 .or. .not.glerlarea) then
                   fld(i,j)=field(i,j)

!          Else treat as a water point
                else
                   fldwter(i,j)=field(i,j)
                end if

             else if (iflg==1) then
! Combine fld,fldwter into field
!          Land point
                if(isli2(i,j) == 1 .or. .not.glerlarea) then
                   field(i,j)=fld(i,j)
!          Treat as a water point,
                else
                   field(i,j)=fldwter(i,j)
                end if
             end if
          end do
       end do

  return
end subroutine anbkgvar_lw

subroutine ansmoothrf(cstate)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ansmoothrf  anisotropic rf for regional mode
!   prgmmr: parrish          org: np22                date: 2005-02-14
!
! abstract: apply anisotropic rf for regional mode
!
! program history log:
!   2005-02-14  parrish
!   2008-12-04  sato - update for global mode
!   2009-01-02  todling - get mype from mpimod directly
!   2015-07-01  pondeca - rewrite to handle motley variables
!
!   input argument list:
!     cstate     - bundle containing horizontal fields to be smoothed
!
!   output argument list:
!     cstate     - bundle containing smoothed horizontal fields
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: indices,indices_p,ngauss,pf2aP1, &
                      filter_all,filter_p2,filter_p3
  use patch2grid_mod, only: patch2grid, tpatch2grid
  use mpimod, only:  npe
  use constants, only: zero
  use gridmod, only: nlat,nlon,lat2,lon2,nsig,regional
  use control_vectors, only: nvars,nrf,nrf_var,nrf_3d
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use general_commvars_mod, only: s2g_raf
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  use fgrid2agrid_mod, only: fgrid2agrid,tfgrid2agrid
  use raflib, only: raf4_ad,raf4
  use raflib, only: raf4_ad_wrap,raf4_wrap
  implicit none

! Declare passed variables
  type(gsi_bundle),intent(inout) :: cstate

! Declare local variables
  integer(i_kind):: ips,ipe,jps,jpe,kps,kpe,kds,kde
  integer(i_kind):: p_ips,p_ipe,p_jps,p_jpe,p_kps,p_kpe
  integer(i_kind) i,igauss,j,k,kk,n,istatus
  real(r_kind),allocatable:: fields(:,:,:,:)
  real(r_kind),allocatable:: work(:,:,:,:)
  real(r_kind),allocatable:: worka(:,:,:)
  real(r_single),allocatable:: workb(:,:,:,:)
  real(r_kind),allocatable,dimension(:,:,:)  :: workanp,workasp
  real(r_single),allocatable,dimension(:,:,:,:):: workbnp,workbsp
  real(r_kind),pointer::rank2(:,:)
  real(r_kind),pointer::rank3(:,:,:)

! Convert from subdomain to full horizontal field distributed among processors
  kds=1
  kde=s2g_raf%num_fields
  kps=s2g_raf%kbegin_loc
  kpe=s2g_raf%kend_loc
  ips=indices%ips
  ipe=indices%ipe
  jps=indices%jps
  jpe=indices%jpe

  p_kps=s2g_raf%kbegin_loc  !indices_p%kps
  p_kpe=s2g_raf%kend_loc    !indices_p%kpe
  p_ips=indices_p%ips
  p_ipe=indices_p%ipe
  p_jps=indices_p%jps
  p_jpe=indices_p%jpe

               !WHY THE "1" IN THE ALLOCATE STATEMENTS? / MPondeca
  allocate(fields(1,lat2,lon2,kds:kde))             !MUST DEALLOCATE / MPondeca
  allocate(work(1,nlat,nlon,kps:max(kps,kpe)))   !MUST DEALLOCATE / MPondeca

  allocate(worka(ips:ipe,jps:jpe,kps:max(kps,kpe)))
  allocate(workb(ngauss,ips:ipe,jps:jpe,kps:max(kps,kpe)))

  if(.not.regional) then
     allocate(workanp(p_ips:p_ipe, p_jps:p_jpe, p_kps:p_kpe))
     allocate(workasp(p_ips:p_ipe ,p_jps:p_jpe ,p_kps:p_kpe))

     allocate(workbnp(ngauss, p_ips:p_ipe, p_jps:p_jpe, p_kps:p_kpe))
     allocate(workbsp(ngauss, p_ips:p_ipe, p_jps:p_jpe, p_kps:p_kpe ))
  end if

  fields=zero
  kk=0
  do n=1,nvars
     if (n<=nrf .and. nrf_3d(min(n,nrf))) then
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
        if(istatus==0) then
           do k=1,nsig
              kk=kk+1
              do j=1,lon2
                 do i=1,lat2
                    fields(1,i,j,kk)=rank3(i,j,k)
                 end do
              end do
           end do
        endif
     else        !2d flds including motley flds
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
        if(istatus==0) then
           kk=kk+1
           do j=1,lon2
              do i=1,lat2
                 fields(1,i,j,kk)=rank2(i,j)
              end do
           end do
        endif
     endif
  end do

! Convert from subdomain to full horizontal fields distributed among processors

  call general_sub2grid(s2g_raf,fields,work)

!  adjoint of coarse to fine grid
  do k=kps,kpe
     if(regional) then
        call tfgrid2agrid(pf2aP1,work(1,1,1,k),worka(ips,jps,k))
     else
        call tpatch2grid(work(1,1,1,k), &
                         worka  (ips,jps,k),   &
                         workanp(p_ips,p_jps,k), &
                         workasp(p_ips,p_jps,k))
     end if
  end do

!  transfer coarse grid fields to ngauss copies
  do k=kps,kpe
     do j=jps,jpe
        do i=ips,ipe
           do igauss=1,ngauss
              workb(igauss,i,j,k)=worka(i,j,k)
           end do
        end do
     end do
  end do

  if(.not.regional) then
     do k=p_kps,p_kpe
        do j=p_jps,p_jpe
           do i=p_ips,p_ipe
              do igauss=1,ngauss
                 workbnp(igauss,i,j,k)=workanp(i,j,k)
                 workbsp(igauss,i,j,k)=workasp(i,j,k)
              end do
           end do
        end do
     end do
  end if

!   apply recursive filter

!  call raf4_wrap(   workb,filter_all,ngauss,indices,npe)  /MPondeca
!  call raf4_ad_wrap(workb,filter_all,ngauss,indices,npe)  /MPondeca

  call raf4 (workb,filter_all,ngauss,ips,ipe,jps,jpe,ips,ipe,jps,jpe,kps,kpe,npe)
  call raf4_ad(workb,filter_all,ngauss,ips,ipe,jps,jpe,ips,ipe,jps,jpe,kps,kpe,npe)

  if(.not.regional) then
     call raf4_wrap(   workbnp,filter_p2,ngauss,indices_p,npe)   !WORK ON THIS / MPondeca
     call raf4_ad_wrap(workbnp,filter_p2,ngauss,indices_p,npe)   !WORK ON THIS / MPondeca
     call raf4_wrap(   workbsp,filter_p3,ngauss,indices_p,npe)   !WORK ON THIS / MPondeca
     call raf4_ad_wrap(workbsp,filter_p3,ngauss,indices_p,npe)   !WORK ON THIS / MPondeca
  end if

!  add together ngauss copies
  worka=zero
  do k=kps,kpe
     do j=jps,jpe
        do i=ips,ipe
           do igauss=1,ngauss
              worka(i,j,k)=worka(i,j,k)+workb(igauss,i,j,k)
           end do
        end do
     end do
  end do
  deallocate(workb)

  if(.not.regional) then
     workanp=zero
     workasp=zero
     do k=p_kps,p_kpe
        do j=p_jps,p_jpe
           do i=p_ips,p_ipe
              do igauss=1,ngauss
                 workanp(i,j,k)=workanp(i,j,k)+workbnp(igauss,i,j,k)
                 workasp(i,j,k)=workasp(i,j,k)+workbsp(igauss,i,j,k)
              end do
           end do
        end do
     end do
  end if

!  coarse to fine grid
  do k=kps,kpe
     if(regional) then
        call fgrid2agrid(pf2aP1,worka(ips,jps,k),work(1,1,1,k))
     else
        call patch2grid(work(1,1,1,k), &
                        worka  (ips,jps,k),&
                        workanp(p_ips,p_jps,k),&
                        workasp(p_ips,p_jps,k))
     end if
  end do
  deallocate(worka)

  call general_grid2sub(s2g_raf,work,fields)
  deallocate(work)

  if(.not.regional) then
     deallocate(workanp)
     deallocate(workbnp)
     deallocate(workasp)
     deallocate(workbsp)
  end if

!   transfer from work back to bundle

  kk=0
  do n=1,nvars
     if (n<=nrf .and. nrf_3d(min(n,nrf))) then
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
        if(istatus==0) then
           do k=1,nsig
              kk=kk+1
              do j=1,lon2
                 do i=1,lat2
                    rank3(i,j,k)=fields(1,i,j,kk)
                 end do
              end do
           end do
        endif
     else        !2d flds including motley flds
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
        if(istatus==0) then
           kk=kk+1
           do j=1,lon2
              do i=1,lat2
                 rank2(i,j)=fields(1,i,j,kk)
              end do
           end do
        endif
     endif
  end do
  deallocate(fields)

end subroutine ansmoothrf

subroutine vert_smther(g,nsmooth,nsmooth_shapiro)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    vert_smther
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-15  lueken - added subprogram doc block
!
!   input argument list:
!    nsmooth
!    nsmooth_shapiro
!    g
!
!   output argument list:
!    g
!
! Notes:  nsmooth > 0 ==> apply 1-2-1 smoother
!         nsmooth_shapiro 0 ==> apply second moment preserving
!                               "shapiro smoother"
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use constants, only: quarter,half
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nsig
  implicit none

! Declare passed variables
  integer(i_kind)                             ,intent(in   ) :: nsmooth,nsmooth_shapiro
  real(r_kind),dimension(1:lat2,1:lon2,1:nsig),intent(inout) :: g

! Declare local variables
  integer(i_kind) i,j,l,k,kp,km,kp3,km3
  real(r_kind), allocatable:: gaux(:)

  if (nsig==1)return

  allocate(gaux(1:nsig))

  if (nsmooth > 0 ) then
     do i=1,lat2
        do j=1,lon2
           do l=1,nsmooth
              gaux(1:nsig)=g(i,j,1:nsig)
              do k=1,nsig
                 kp=min(k+1,nsig) ; km=max(1,k-1)
                 g(i,j,k)=quarter*(gaux(kp)+gaux(km))+half*gaux(k)
              enddo
           enddo
        enddo
     enddo
  endif

  if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then
     do i=1,lat2
        do j=1,lon2
           do l=1,nsmooth_shapiro
              gaux(1:nsig)=g(i,j,1:nsig)
              do k=1,nsig
                 kp=min(k+1,nsig) ; km=max(1,k-1)
                 kp3=min(k+3,nsig) ; km3=max(1,k-3)
                 g(i,j,k)=.28125_r_kind*(gaux(kp)+gaux(km))+half*gaux(k)-.03125_r_kind*(gaux(kp3)+gaux(km3))
              enddo
           enddo
        enddo
     enddo
  endif
  deallocate(gaux)

  return
end subroutine vert_smther

subroutine tvert_smther(g,nsmooth,nsmooth_shapiro)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    tvert_smther
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-09-15  lueken - added subprogram doc block
!
!   input argument list:
!    nsmooth
!    nsmooth_shapiro
!    g
!
!   output argument list:
!    g
!
! Notes:  nsmooth > 0 ==>  1-2-1 smoother
!         nsmooth_shapiro 0 ==>  second moment preserving
!                               "shapiro smoother"
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use constants, only: zero,quarter,half
  use kinds, only: r_kind,i_kind
  use gridmod, only: lat2,lon2,nsig
  implicit none

! Declare passed variables
  integer(i_kind)                             ,intent(in   ) :: nsmooth,nsmooth_shapiro
  real(r_kind),dimension(1:lat2,1:lon2,1:nsig),intent(inout) :: g

! Declare local variables
  integer(i_kind) i,j,l,k,kp,km,kp3,km3
  real(r_kind), allocatable:: gaux(:)

  if (nsig==1)return

  allocate(gaux(1:nsig))

  if (nsmooth > 0 ) then
     do i=1,lat2
        do j=1,lon2
           do l=1,nsmooth
              gaux(1:nsig)=zero
              do k=1,nsig
                 kp=min(k+1,nsig) ; km=max(1,k-1)
                 gaux(k)=gaux(k)+half*g(i,j,k)
                 gaux(km)=gaux(km)+quarter*g(i,j,k)
                 gaux(kp)=gaux(kp)+quarter*g(i,j,k)
              enddo
              g(i,j,1:nsig)=gaux(1:nsig)
           enddo
        enddo
     enddo
  endif

  if (nsmooth_shapiro > 0 .and. nsmooth <= 0) then
     do i=1,lat2
        do j=1,lon2
           do l=1,nsmooth_shapiro
              gaux(1:nsig)=zero
              do k=1,nsig
                 kp=min(k+1,nsig) ; km=max(1,k-1)
                 kp3=min(k+3,nsig) ; km3=max(1,k-3)
                 gaux(km3)=gaux(km3)-.03125_r_kind*g(i,j,k)
                 gaux(kp3)=gaux(kp3)-.03125_r_kind*g(i,j,k)
                 gaux(k)=gaux(k)+half*g(i,j,k)
                 gaux(km)=gaux(km)+.28125_r_kind*g(i,j,k)
                 gaux(kp)=gaux(kp)+.28125_r_kind*g(i,j,k)
              enddo
              g(i,j,1:nsig)=gaux(1:nsig)
           enddo
        enddo
     enddo
  endif
  deallocate(gaux)

  return
end subroutine tvert_smther

subroutine ansmoothrf_reg_sub2slab_option(cstate)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ansmoothrf_reg_subdomain_option  anisotropic rf for regional mode
!   prgmmr: parrish          org: np22                date: 2015-05-01
!
! abstract: apply anisotropic rf for regional mode (input subdomains but slab for filter)
!
! program history log:
!   2005-02-14  parrish
!   2011-02-22  zhu - use cstate to replace argument list such as p,t,q,vp,st
!   2014-02-07  pondeca - update to include motley variables as well in 2d set of filtered variables
!   2015-05-02  parrish - make copy of ansmoothrf_reg_subdomain_option and
!                            modify to allow conversion from subdomains to slabs.
!
!   input argument list:
!     t,p,q,oz,st,stl,sti,cwmr,st,vp   -  fields to be smoothed
!
!   output argument list:
!     t,p,q,oz,st,stl,sti,cwmr,st,vp   -  smoothed fields
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: filter_all,ngauss
  use anberror, only: pf2aP1
  use mpimod, only: npe
  use constants, only: zero,zero_single
  use gridmod, only: lat2,lon2,nsig
  use raflib, only: raf4_ad,raf4
  use control_vectors, only: nvars,nrf,nrf_var,nrf_3d
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use general_commvars_mod, only: s2g_raf
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  use fgrid2agrid_mod, only: fgrid2agrid,tfgrid2agrid
  implicit none

! Declare passed variables
  type(gsi_bundle),intent(inout) :: cstate

! Declare local variables
  integer(i_kind) i,igauss,j,k,kk,n,istatus
  real(r_kind),allocatable:: worka(:,:,:,:)
  real(r_single),allocatable:: workb(:,:,:,:)
  real(r_kind),pointer::rank2(:,:)
  real(r_kind),pointer::rank3(:,:,:)
  real(r_kind),allocatable::slaba(:,:,:,:)
  real(r_kind),allocatable::slabb(:,:,:)

  integer(i_kind):: kds,kde,kps,kpe,nlatf,nlonf,nlat,nlon

  kds=1
  kde=s2g_raf%num_fields
  kps=s2g_raf%kbegin_loc
  kpe=s2g_raf%kend_loc
  nlat=pf2ap1%nlata
  nlon=pf2ap1%nlona
  nlatf=pf2ap1%nlatf
  nlonf=pf2ap1%nlonf

  allocate(worka(1,lat2,lon2,kds:kde))
  allocate(slaba(1,nlat,nlon,kps:max(kps,kpe)))
  allocate(slabb(nlatf,nlonf,kps:max(kps,kpe)))
  allocate(workb(ngauss,nlatf,nlonf,kps:max(kps,kpe)))

!   transfer from bundle to worka

  worka=zero
  kk=0
  do n=1,nvars
     if (n<=nrf .and. nrf_3d(min(n,nrf))) then
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
        if(istatus==0) then
           do k=1,nsig
              kk=kk+1
              do j=1,lon2
                 do i=1,lat2
                    worka(1,i,j,kk)=rank3(i,j,k)
                 end do
              end do
           end do
        endif
     else        !2d flds including motley flds
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
        if(istatus==0) then
           kk=kk+1
           do j=1,lon2
              do i=1,lat2
                 worka(1,i,j,kk)=rank2(i,j)
              end do
           end do
        endif
     endif
  end do

!   transfer from subdomains to slabs

  call general_sub2grid(s2g_raf,worka,slaba)

!  adjoint of coarse to fine grid interpolation:

  do k=kps,kpe
     call tfgrid2agrid(pf2aP1,slaba(1,1,1,k),slabb(1,1,k))
  end do

!  transfer to ngauss copies

  do k=kps,kpe
     do j=1,nlonf
        do i=1,nlatf
           do igauss=1,ngauss
              workb(igauss,i,j,k)=slabb(i,j,k)
           end do
        end do
     end do
  end do

!   apply recursive filter

  call raf4   (workb,filter_all,ngauss,1,nlatf,1,nlonf,1,nlatf,1,nlonf,kps,kpe,npe)
  call raf4_ad(workb,filter_all,ngauss,1,nlatf,1,nlonf,1,nlatf,1,nlonf,kps,kpe,npe)

!  add together ngauss copies

  slabb=zero
  do k=kps,kpe
     do j=1,nlonf
        do i=1,nlatf
           do igauss=1,ngauss
              slabb(i,j,k)=workb(igauss,i,j,k)+slabb(i,j,k)
           end do
        end do
     end do
  end do
  deallocate(workb)

!  coarse to fine grid interpolation:

  do k=kps,kpe
     call fgrid2agrid(pf2aP1,slabb(1,1,k),slaba(1,1,1,k))
  end do
  deallocate(slabb)

  call general_grid2sub(s2g_raf,slaba,worka)
  deallocate(slaba)

!   transfer from worka back to bundle

  kk=0
  do n=1,nvars
     if (n<=nrf .and. nrf_3d(min(n,nrf))) then
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
        if(istatus==0) then
           do k=1,nsig
              kk=kk+1
              do j=1,lon2
                 do i=1,lat2
                    rank3(i,j,k)=worka(1,i,j,kk)
                 end do
              end do
           end do
        endif
     else        !2d flds including motley flds
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
        if(istatus==0) then
           kk=kk+1
           do j=1,lon2
              do i=1,lat2
                 rank2(i,j)=worka(1,i,j,kk)
              end do
           end do
        endif
     endif
  end do
  deallocate(worka)

end subroutine ansmoothrf_reg_sub2slab_option

subroutine ansmoothrf_reg_subdomain_option(cstate)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    ansmoothrf_reg_subdomain_option  anisotropic rf for regional mode
!   prgmmr: parrish          org: np22                date: 2005-02-14
!
! abstract: apply anisotropic rf for regional mode (using subdomains instead of slabs)
!              NOTE: only works if filter grid is same as analysis grid
!
! program history log:
!   2005-02-14  parrish
!   2011-02-22  zhu - use cstate to replace argument list such as p,t,q,vp,st
!   2014-02-07  pondeca - update to include motley variables as well in 2d set of filtered variables
!
!   input argument list:
!     t,p,q,oz,st,stl,sti,cwmr,st,vp   -  fields to be smoothed
!
!   output argument list:
!     t,p,q,oz,st,stl,sti,cwmr,st,vp   -  smoothed fields
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!$$$
  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: indices, filter_all,ngauss,halo_update_reg
  use mpimod, only: mype,npe
  use constants, only: zero,zero_single
  use gridmod, only: istart,jstart,nsig
  use raflib, only: raf4_ad_wrap,raf4_wrap
  use control_vectors, only: nvars,nrf,nrf_var,nrf_3d
  use gsi_bundlemod, only: gsi_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  implicit none

! Declare passed variables
  type(gsi_bundle),intent(inout) :: cstate

! Declare local variables
  integer(i_kind) i,igauss,iloc,j,jloc,k,kk,mm1,n,istatus
  real(r_single),dimension(ngauss, &
                           indices%ips:indices%ipe,&
                           indices%jps:indices%jpe,&
                           indices%kps:indices%kpe):: workb
  real(r_kind),pointer::rank2(:,:)
  real(r_kind),pointer::rank3(:,:,:)

  integer(i_kind):: ids,ide,jds,jde,kds,kde,ips,ipe,jps,jpe,kps,kpe

  ids=indices%ids; ide=indices%ide
  jds=indices%jds; jde=indices%jde
  kds=indices%kds; kde=indices%kde
  ips=indices%ips; ipe=indices%ipe
  jps=indices%jps; jpe=indices%jpe
  kps=indices%kps; kpe=indices%kpe

  mm1=mype+1

  workb=zero_single

!  transfer variables to ngauss copies
  kk=0
  do n=1,nvars
     if (n<=nrf .and. nrf_3d(min(n,nrf))) then
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
        if(istatus==0) then
           do k=1,nsig
              kk=kk+1
              do j=jps,jpe
                 jloc=j-jstart(mm1)+2
                 do i=ips,ipe
                    iloc=i-istart(mm1)+2
                    do igauss=1,ngauss
                       workb(igauss,i,j,kk)=rank3(iloc,jloc,k)
                    end do
                 end do
              end do
           end do
        endif
     else        !2d flds including motley flds
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
        if(istatus==0) then
           kk=kk+1
           do j=jps,jpe
              jloc=j-jstart(mm1)+2
              do i=ips,ipe
                 iloc=i-istart(mm1)+2
                 do igauss=1,ngauss
                    workb(igauss,i,j,kk)=rank2(iloc,jloc)
                 end do
              end do
           end do
        endif
     endif
  end do

!   apply recursive filter

  call raf4_wrap(workb,filter_all,ngauss,indices,npe)
  call raf4_ad_wrap(workb,filter_all,ngauss,indices,npe)

!  add together ngauss copies
  kk=0
  do n=1,nvars
     if (n<=nrf .and. nrf_3d(min(n,nrf))) then
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank3,istatus)
        if(istatus==0) then
           do k=1,nsig
              kk=kk+1
              do j=jps,jpe
                 jloc=j-jstart(mm1)+2
                 do i=ips,ipe
                    iloc=i-istart(mm1)+2
                    rank3(iloc,jloc,k)=zero
                    do igauss=1,ngauss
                       rank3(iloc,jloc,k)=rank3(iloc,jloc,k)+workb(igauss,i,j,kk)
                    end do
                 end do
              end do
           end do
           call halo_update_reg(rank3,nsig)
        endif
     else
        call gsi_bundlegetpointer (cstate,trim(nrf_var(n)),rank2,istatus)
        if(istatus==0) then
           kk=kk+1
           do j=jps,jpe
              jloc=j-jstart(mm1)+2
              do i=ips,ipe
                 iloc=i-istart(mm1)+2
                 rank2(iloc,jloc)=zero
                 do igauss=1,ngauss
                    rank2(iloc,jloc)=rank2(iloc,jloc)+workb(igauss,i,j,kk)
                 end do
              end do
           end do
           call halo_update_reg(rank2,1)
        endif
     endif
  end do

end subroutine ansmoothrf_reg_subdomain_option


text:
      subroutine fill_old_items(old_items)

!     Deltares Software Centre

!>/File
!>      add the old items to the collection

      use timers         !< performance timers
      use processet      !< use processet definitions
      implicit none

      ! decalaration of arguments

      type(old_item_coll)                :: old_items         !< the old_items table to be filled

      ! common declarations

      include 'data.inc'                           ! tables read from proces definition file

      ! local declaration

      type(old_item)            :: a_old_item        ! single process
      integer                   :: i                 ! loop counter old items
      integer                   :: i2                ! index in collection
      integer(4)                :: ithndl = 0        ! handle for performance timer
      if (timon) call timstrt( "fill_old_items", ithndl )

      do i = 1, n_old_items

         a_old_item%old_name      = old_items_old_name(i)
         a_old_item%new_name      = old_items_new_name(i)
         a_old_item%old_default   = old_items_old_default(i)
         a_old_item%configuration = old_items_configuration(i)
         a_old_item%serial        = old_items_serial(i)
         a_old_item%action_type   = old_items_action_type(i)

         i2 = old_item_coll_add( old_items, a_old_item)

      enddo

      if (timon) call timstop( ithndl )
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
The http://nrs.ucop.edu/ UC Natural Reserve System protects over 135,000 acres of land in California for research and educational purposes. Most (not all) reserves have restricted access, meaning you either need to go with an organized group or http://nrs.ucop.edu/reserves/UseApps/UseApps.html apply for entrance. Each and every one of them is a very special place and well worth a visit for some Outdoor Activities outdoor fun.

 UC Davis Reserves

 Other notable reserves
  Sagehen Creek Field Station
  Russell Ranch Wildlife Area
  See also: other http://nrs.ucop.edu/ReserveMapCampus.htm UC Reserves listed by campus