From 6dee3922c34a1a9a1e4ebc5fa5f266ffa07afd7c Mon Sep 17 00:00:00 2001
From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com>
Date: Tue, 30 Jan 2024 08:32:33 -0500
Subject: [PATCH] orog.fd - Reduce dependency on SPLIB and some cleanup (#889)

- Remove routines MAKEMT, MAKEPC and MAKEOA, which were used for the old spectral GFS
gaussian grid (and have dependencies on SPLIB) and are now obsolete.
- Remove write of the GRIB1 version of the 'orog' file. That was only used for the old spectral
GFS and had dependencies on SPLIB.
- Remove write of old binary version of the 'orog' file, which was used for the old spectral GFS.
- Remove BACIO library from the build.
- Remove unused variables and perform some general cleanup.

Fixes #886.
---
 reg_tests/grid_gen/driver.hercules.sh         |    2 +-
 reg_tests/weight_gen/driver.hercules.sh       |    2 +-
 .../orog_mask_tools.fd/orog.fd/CMakeLists.txt |    1 -
 .../orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F | 1451 ++---------------
 sorc/orog_mask_tools.fd/orog.fd/netcdf_io.F90 |    4 +-
 ush/fv3gfs_make_orog.sh                       |    8 +-
 6 files changed, 112 insertions(+), 1356 deletions(-)

diff --git a/reg_tests/grid_gen/driver.hercules.sh b/reg_tests/grid_gen/driver.hercules.sh
index 3c46a2d45..5a323249c 100755
--- a/reg_tests/grid_gen/driver.hercules.sh
+++ b/reg_tests/grid_gen/driver.hercules.sh
@@ -30,7 +30,7 @@ module list
 set -x
 ulimit -s unlimited
 
-export WORK_DIR="${WORK_DIR:-/work/noaa/stmp/$LOGNAME}"
+export WORK_DIR="${WORK_DIR:-/work2/noaa/stmp/$LOGNAME}"
 export WORK_DIR="${WORK_DIR}/reg-tests/grid-gen"
 QUEUE="${QUEUE:-batch}"
 PROJECT_CODE=${PROJECT_CODE:-fv3-cpu}
diff --git a/reg_tests/weight_gen/driver.hercules.sh b/reg_tests/weight_gen/driver.hercules.sh
index 369796d61..5fba04abe 100755
--- a/reg_tests/weight_gen/driver.hercules.sh
+++ b/reg_tests/weight_gen/driver.hercules.sh
@@ -36,7 +36,7 @@ module use ../../modulefiles
 module load build.$target.$compiler
 module list
 
-export DATA="${WORK_DIR:-/work/noaa/stmp/$LOGNAME}"
+export DATA="${WORK_DIR:-/work2/noaa/stmp/$LOGNAME}"
 export DATA="${DATA}/reg-tests/weight_gen"
 
 #-----------------------------------------------------------------------------
diff --git a/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt b/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
index 04ab86742..bc0e5d5b1 100644
--- a/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
+++ b/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
@@ -25,7 +25,6 @@ target_include_directories(orog_lib INTERFACE ${mod_dir})
 target_link_libraries(
   orog_lib
   PUBLIC
-  bacio::bacio_4
   w3emc::w3emc_d
   ip::ip_d
   sp::sp_d
diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
index 96e3b38d8..add54f00c 100644
--- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
+++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
@@ -1,25 +1,9 @@
 C> @file
-C> Terrain maker for global spectral model.
+C> Terrain maker for the ufs weather model.
 C> @author Mark Iredell @date 92-04-16
   
-C> This program creates 7 terrain-related files computed from the navy
-C> 10-minute terrain dataset. The model physics grid parameters and
-C> spectral truncation and filter parameters are read by this program as
-C> input.
-C>      
-C> The 7 files produced are:
-C> 1. sea-land mask on model physics grid
-C> 2. gridded orography on model physics grid
-C> 3. mountain std dev on model physics grid
-C> 4. spectral orography in spectral domain
-C> 5. unfiltered gridded orography on model physics grid
-C> 6. grib sea-land mask on model physics grid
-C> 7. grib gridded orography on model physics grid
-C>
-C> The orography is only filtered for wavenumbers greater than nf0. For
-C> wavenumbers n between nf0 and nf1, the orography is filtered by the
-C> factor 1-((n-nf0)/(nf1-nf0))**2. The filtered orography will not have
-C> information beyond wavenumber nf1.
+C> This program creates land mask and fraction, terrain and gravity wave
+C> drag fields for a single tile of the cubed sphere grid. 
 C>
 C> PROGRAM HISTORY LOG:
 C> -  92-04-16  IREDELL
@@ -37,90 +21,76 @@
 C>  - 05-09-05   if test on HK and HLPRIM for GAMMA SQRT
 C>  - 07-08-07   replace 8' with 30" incl GICE, conintue w/ S-Y. lake slm
 C>  - 08-08-07  All input 30", UMD option, and filter as described below
-C> Quadratic filter applied by default.
-C> NF0 is normally set to an even value beyond the previous truncation, 
-C> for example, for jcap=382, NF0=254+2
-C> NF1 is set as jcap+2 (and/or nearest even), eg., for t382, NF1=382+2=384
-C> if no filter is desired then NF1=NF0=0 and ORF=ORO
-C> but if no filter but spectral to grid (with gibbs) then NF1=jcap+2, and NF1=jcap+1
 C>       
 C>   INPUT FILES:
-C>  -   UNIT5      - PHYSICS LONGITUDES (IM), PHYSICS LATITUDES (JM),
-C>                  SPECTRAL TRUNCATION (NM), RHOMBOIDAL FLAG (NR),
-C>                  AND FIRST AND SECOND FILTER PARAMETERS (NF0,NF1).
-C>                  RESPECTIVELY READ IN FREE FORMAT.
+C>  -   UNIT5      - RESOLUTION OF INPUT OROG DATA (MTNRES), I-DIMENSION
+C>                   OF CUBED SPHERE TILE (IM), J-DIMENSION OF TILE
+C>                   (JM), FACTOR TO ADJUST OROGRAPHY BY ITS VARIABLE
+C>                   (EFAC), LATITUDE REVERSAL FLAG (BLAT), 
+C>                   RESPECTIVELY READ IN FREE FORMAT.
 C>  -   UNIT235    - GTOPO 30" AVR for ZAVG elevation
 C>  -   UNIT10     - 30" UMD land (lake) cover mask  see MSKSRC switch
-C>  -  XUNIT11     - GTOPO AVR
-C>  -  XUNIT12     - GTOPO STD DEV
-C>  -  XUNIT13     - GTOPO MAX
 C>  -   UNIT14     - GTOPO SLM (10' NAVY if switched to get lakes 
 C>  -   UNIT15     - GICE Grumbine 30" RAMP Antarctica orog IMNx3616
 C>  -   UNIT25     - Ocean land-sea mask on gaussian grid         
+C>  -   NCID       - 'GRID' FILE FOR THE CUBED SPHERE TILE (OUTGRID).
+C>  -   NCID       - (INPUTOROG) Input orography/GWD file on gaussian
+C>                   grid. When specified, will be interpolated to model tile.
+C>  -   NCID       - Mask/land fraction on the cubed sphere tile.   
+C>                   (optional).
 C>
 C>   OUTPUT FILES:
-C>  -   UNIT51     - SEA-LAND MASK (IM,JM)
-C>  -   UNIT52     - GRIDDED OROGRAPHY (IM,JM)
-C>  -   UNIT54     - SPECTRAL OROGRAPHY ((NM+1)*((NR+1)*NM+2))
-C>  -   UNIT55     - UNFILTERED GRIDDED OROGRAPHY (IM,JM)
-C>  -   UNIT57     - GRIB GRIDDED OROGRAPHY (IM,JM)
+C>  -   NCID       - Mask/land fraction on the cubed sphere tile.   
+C>                   (optional).
+C>  -   NCID       - The final orography file on the cubed sphere tile.
+C>                   Contains mask, terrain and GWD fields.
 C>
 C>   SUBPROGRAMS CALLED:
-C>  -   UNIQUE:
 C>  -   TERSUB     - MAIN SUBPROGRAM
-C>  -   SPLAT      - COMPUTE GAUSSIAN LATITUDES OR EQUALLY-SPACED LATITUDES
-C>  -   LIBRARY:
-C>  -   SPTEZ      - SPHERICAL TRANSFORM
-C>  -   GBYTES     - UNPACK BITS
 C>
 C> @return 0 for success, error code otherwise.
+      implicit none
       include 'netcdf.inc'
       logical fexist, opened
-      integer fsize, ncid, error, id_dim, nx, ny
+      integer fsize, ncid, error, id_dim, nx, ny, imn, jmn
       character(len=256) :: OUTGRID = "none"
       character(len=256) :: INPUTOROG = "none"
       character(len=256) :: merge_file = "none"
       logical :: mask_only = .false.
-      integer :: MTNRES,IM,JM,NM,NR,NF0,NF1,EFAC,BLAT,NW
+      integer :: MTNRES,IM,JM,EFAC,BLAT
       fsize=65536
-      READ(5,*) MTNRES,IM,JM,NM,NR,NF0,NF1,EFAC,BLAT
+      READ(5,*) MTNRES,IM,JM,EFAC,BLAT
       READ(5,*) OUTGRID
       READ(5,*) INPUTOROG
       READ(5,*) mask_only
       READ(5,*) merge_file
-!      MTNRES=1
-!      IM=48
-!      JM=48
-!      NM=46
-!      NF0=0             
-!      NF1=0              
-!      efac=0
-!      blat=0
-!      NR=0
-!      OUTGRID = "C48_grid.tile1.nc"
-!      INPUTOROG = "oro.288x144.nc"
-      print*, "INPUTOROG=", trim(INPUTOROG)
-      print*, "IM,JM=", IM, JM
-      print*, "MASK_ONLY", mask_only
-      print*, "MERGE_FILE", trim(merge_file)
+
+      print*, "INPUTOROG: ", trim(INPUTOROG)
+      print*, "OUTGRID: ", trim(INPUTOROG)
+      print*, "IM,JM: ", IM, JM
+      print*, "MASK_ONLY: ", mask_only
+      print*, "MERGE_FILE: ", trim(merge_file)
+
 ! --- MTNRES defines the input (highest) elev resolution
 ! --- =1 is topo30 30" in units of 1/2 minute.
 !     so MTNRES for old values must be *2.
 !     =16 is now Song Yu's 8' orog the old ops standard
 ! --- other possibilities are =8 for 4' and =4 for 2' see
 !     HJ for T1000 test. Must set to 1 for now.
-      MTNRES=1
-      print*, MTNRES,IM,JM,NM,NR,NF0,NF1,EFAC,BLAT
-      NW=(NM+1)*((NR+1)*NM+2)
+
+      print*, "MTNRES,EFAC,BLAT: ",MTNRES,EFAC,BLAT
+
       IMN = 360*120/MTNRES
       JMN = 180*120/MTNRES
-      print *, ' Starting terr12 mtnlm7_slm30.f  IMN,JMN:',IMN,JMN
 
-! --- read the grid resolution if the OUTGRID exists.
-      if( trim(OUTGRID) .NE. "none" ) then
+! --- Read the model grid resolution from 'grid' file OUTGRID.
+
+      READ_GRID_FILE : if( trim(OUTGRID) .NE. "none" ) then
+
          inquire(file=trim(OUTGRID), exist=fexist)
          if(.not. fexist) then
-            print*, "file "//trim(OUTGRID)//" does not exist"
+            print*, "FATAL ERROR: file "
+     &             //trim(OUTGRID)//" does not exist"
             CALL ERREXIT(4)
          endif
          do ncid = 103, 512
@@ -128,7 +98,7 @@
            if( .NOT.opened )exit
          end do
 
-         print*, "outgrid=", trim(outgrid)
+         print*, "READ outgrid=", trim(outgrid)
          error=NF__OPEN(trim(OUTGRID),NF_NOWRITE,fsize,ncid)
          call netcdf_err(error, 'Open file '//trim(OUTGRID) )
          error=nf_inq_dimid(ncid, 'nx', id_dim)
@@ -157,10 +127,15 @@
          endif
          error=nf_close(ncid)
          call netcdf_err(error, 'close file '//trim(OUTGRID) )
+
+      else
+
+        print*, "FATAL ERROR: SPECIFY OUTGRID"
+        CALL ERREXIT(8)
          
-      endif         
+      endif READ_GRID_FILE
          
-      CALL TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
+      CALL TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
      &            OUTGRID,INPUTOROG,MASK_ONLY,MERGE_FILE)
       STOP
       END
@@ -171,11 +146,6 @@
 !! @param[in] JMN "j" dimension of the input terrain dataset.
 !! @param[in] IM "i" dimension of the model grid tile.
 !! @param[in] JM "j" dimension of the model grid tile.
-!! @param[in] NM Spectral truncation.
-!! @param[in] NR Rhomboidal flag.
-!! @param[in] NF0 First order spectral filter parameters.
-!! @param[in] NF1 Second order spectral filter parameters.
-!! @param[in] NW Number of waves.
 !! @param[in] EFAC Factor to adjust orography by its variance.
 !! @param[in] BLAT When less than zero, reverse latitude/
 !! longitude for output.
@@ -187,12 +157,12 @@
 !! @param[in] MASK_ONLY Flag to generate the Land Mask only
 !! @param[in] MERGE_FILE Ocean merge file 
 !! @author Jordan Alpert NOAA/EMC
-      SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
+      SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
      &     OUTGRID,INPUTOROG,MASK_ONLY,MERGE_FILE)
       implicit none
       include 'netcdf.inc'
 C
-      integer                      :: IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW
+      integer, intent(in)          :: IMN,JMN,IM,JM
       character(len=*), intent(in) :: OUTGRID
       character(len=*), intent(in) :: INPUTOROG
       character(len=*), intent(in) :: MERGE_FILE
@@ -200,20 +170,19 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
       logical, intent(in) :: mask_only
 
       real, parameter :: MISSING_VALUE=-9999.
-      real, PARAMETER :: PI=3.1415926535897931
       integer, PARAMETER :: NMT=14
 
       integer :: efac,blat,zsave1,zsave2,itopo,kount
       integer :: kount2,islmx,jslmx,oldslm,msksrc,mskocn,notocn
       integer :: i,j,nx,ny,ncid,js,jn,iw,ie,k,it,jt,i1,error,id_dim
       integer :: id_var,nx_in,ny_in,fsize,wgta,IN,INW,INE,IS,ISW,ISE
-      integer :: M,N,IMT,IRET,ios,iosg,latg2,istat,itest,jtest
+      integer :: IMT,ios,iosg,latg2,istat,itest,jtest
       integer :: i_south_pole,j_south_pole,i_north_pole,j_north_pole
       integer :: maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
       integer(1) :: i3save
       integer(2) :: i2save
 
-      integer, allocatable :: JST(:),JEN(:),KPDS(:),KGDS(:),numi(:)
+      integer, allocatable :: numi(:)
       integer, allocatable :: lonsperlat(:)
 
       integer, allocatable :: IST(:,:),IEN(:,:),ZSLMX(:,:)
@@ -223,12 +192,11 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
 
       integer, allocatable :: IWORK(:,:,:)
 
-      real :: DEGRAD,maxlat, minlat,timef,tbeg,tend,tbeg1
-      real :: PHI,DELXN,RS,RN,slma,oroa,vara,var4a,xn,XS,FFF,WWW
-      real :: sumdif,avedif
+      real :: maxlat, minlat,timef,tbeg,tend,tbeg1
+      real :: DELXN,RS,RN,slma,oroa,vara,var4a,xn,XS
 
-      real, allocatable :: COSCLT(:),WGTCLT(:),RCLT(:),XLAT(:),DIFFX(:)
-      real, allocatable :: XLON(:),ORS(:),oaa(:),ola(:),GLAT(:)
+      real, allocatable :: WGTCLT(:),XLAT(:)
+      real, allocatable :: XLON(:),oaa(:),ola(:),GLAT(:)
 
       real, allocatable :: GEOLON(:,:),GEOLON_C(:,:),DX(:,:)
       real, allocatable :: GEOLAT(:,:),GEOLAT_C(:,:),DY(:,:)
@@ -245,28 +213,23 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
       real, allocatable :: OA(:,:,:),OL(:,:,:),HPRIME(:,:,:)
       real, allocatable :: oa_in(:,:,:), ol_in(:,:,:)
 
-
-      complex :: ffj(im/2+1)
-
-      logical :: grid_from_file,output_binary,fexist,opened
-      logical :: SPECTR, REVLAT, FILTER
+      logical :: fexist,opened
+      logical :: REVLAT
 
       logical :: is_south_pole(IM,JM), is_north_pole(IM,JM)
-      logical :: LB(IM*JM)
 
-      output_binary = .false.
       tbeg1=timef()
       tbeg=timef()
       fsize = 65536
 ! integers
-      allocate (JST(JM),JEN(JM),KPDS(200),KGDS(200),numi(jm))
+      allocate (numi(jm))
       allocate (lonsperlat(jm/2))
       allocate (IST(IM,jm),IEN(IM,jm),ZSLMX(2700,1350))
       allocate (glob(IMN,JMN))
 
 ! reals
-      allocate (COSCLT(JM),WGTCLT(JM),RCLT(JM),XLAT(JM),DIFFX(JM/2))
-      allocate (XLON(IM),ORS(NW),oaa(4),ola(4),GLAT(JMN))
+      allocate (WGTCLT(JM),XLAT(JM))
+      allocate (XLON(IM),oaa(4),ola(4),GLAT(JMN))
 
       allocate (ZAVG(IMN,JMN))
       allocate (ZSLM(IMN,JMN))
@@ -275,9 +238,6 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
 !
 !  SET CONSTANTS AND ZERO FIELDS
 !
-      DEGRAD = 180./PI
-      SPECTR = NM .GT. 0     ! if NM <=0 grid is assumed lat/lon
-      FILTER = .TRUE.        ! Spectr Filter defaults true and set by NF1 & NF0
 !     MSKSRC = 0             ! MSKSRC=0 navy 10 lake msk, 1 UMD 30, -1 no lakes
       MSKSRC = 1             ! MSKSRC=0 navy 10 lake msk, 1 UMD 30, -1 no lakes
       REVLAT = BLAT .LT. 0   ! Reverse latitude/longitude for output
@@ -360,7 +320,7 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
 C- READ_G for global 30" terrain 
 C
        print *,' About to call read_g, ITOPO=',ITOPO
-          if ( ITOPO .ne. 0 ) call read_g(glob,ITOPO)
+          if ( ITOPO .ne. 0 ) call read_g(glob)
 ! --- transpose even though glob 30" is from S to N and NCEP std is N to S
        do j=1,jmn/2
        do I=1,imn
@@ -383,8 +343,8 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
 !
 
 !  --- IMN,JMN
-      print*, ' IM, JM, NM, NR, NF0, NF1, EFAC, BLAT'
-      print*, IM,JM,NM,NR,NF0,NF1,EFAC,BLAT
+      print*, ' IM, JM, EFAC, BLAT'
+      print*, IM,JM,EFAC,BLAT
        print *,'  imn,jmn,glob(imn,jmn)=',imn,jmn,glob(imn,jmn)
        print *,' UBOUND ZAVG=',UBOUND(ZAVG)
        print *,' UBOUND glob=',UBOUND(glob)
@@ -507,48 +467,14 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
         enddo
         print *,ios,latg2,'COMPUTE TERRAIN ON A REDUCED GAUSSIAN GRID',
      &          numi
-C       print *,ios,latg2,'COMPUTE TERRAIN ON A REDUCED GAUSSIAN GRID'
       endif
-!       print *,ios,latg2,'TERRAIN ON GAUSSIAN GRID',numi
-      
 !
 !    This code assumes that lat runs from north to south for gg!
 !
-
-      print *,' SPECTR=',SPECTR,' REVLAT=',REVLAT,' ** with GICE-07 **'
-      IF (SPECTR) THEN
-        CALL SPLAT(4,JM,COSCLT,WGTCLT)
-        DO J=1,JM/2
-          RCLT(J)      = ACOS(COSCLT(J))
-        ENDDO
-        DO J = 1,JM/2
-           PHI = RCLT(J) * DEGRAD
-           XLAT(J) = 90. - PHI
-           XLAT(JM-J+1) =  PHI - 90.
-        ENDDO
-      ELSE
-        CALL SPLAT(0,JM,COSCLT,WGTCLT)
-        DO J=1,JM
-          RCLT(J) = ACOS(COSCLT(J))
-          XLAT(J) = 90.0 - RCLT(J) * DEGRAD
-        ENDDO
-      ENDIF
+      print *,' REVLAT=',REVLAT,' ** with GICE-07 **'
 
       allocate (GICE(IMN+1,3601))
 !
-       sumdif = 0.
-       DO J = JM/2,2,-1
-       DIFFX(J) = xlat(J) - XLAT(j-1)
-       sumdif = sumdif + DIFFX(J)
-       ENDDO
-       avedif=sumdif/(float(JM/2))
-!      print *,' XLAT= avedif: ',avedif
-!      write (6,107) (xlat(J)-xlat(j-1),J=JM,2,-1)
-       print *,' XLAT='
-       write (6,106) (xlat(J),J=JM,1,-1)
-  106 format( 10(f7.3,1x))   
-  107 format( 10(f9.5,1x))   
-C
       DELXN = 360./IMN      ! MOUNTAIN DATA RESOLUTION
 C
       DO J=1,JMN
@@ -666,15 +592,15 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
       allocate (land_frac(IM,JM),lake_frac(IM,JM))
 
 !--- reading grid file.
-      grid_from_file = .false.
       is_south_pole = .false.
       is_north_pole = .false.
       i_south_pole = 0
       j_south_pole = 0
       i_north_pole = 0
       j_north_pole = 0
-      if( trim(OUTGRID) .NE. "none" ) then
-         grid_from_file = .true.
+
+! Read 'grid' file
+
          inquire(file=trim(OUTGRID), exist=fexist)
          if(.not. fexist) then
             print*, "file "//trim(OUTGRID)//" does not exist"
@@ -693,25 +619,7 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
      &                   trim(OUTGRID) )
          nx = 2*IM
          ny = 2*JM
-!         error=nf_inq_dimlen(ncid,id_dim,nx)
-!         print*, "nx = ", nx, id_dim
-!         call netcdf_err(error, 'inquire dimension nx length '//
-!     &       'from file '//trim(OUTGRID) )
-!         error=nf_inq_dimid(ncid, 'ny', id_dim)
-!         call netcdf_err(error, 'inquire dimension ny from file '//
-!     &                   trim(OUTGRID) )
-!         error=nf_inq_dimlen(ncid,id_dim,ny)
-!         call netcdf_err(error, 'inquire dimension ny length '//
-!     &       'from file '//trim(OUTGRID) )
-!        IM should equal nx/2 and JM should equal ny/2
-!         if(IM .ne. nx/2) then
-!            print*, "IM=",IM, " /= grid file nx/2=",nx/2
-!            CALL ERREXIT(4)
-!         endif
-!         if(JM .ne. ny/2) then
-!            print*, "JM=",JM, " /= grid file ny/2=",ny/2
-!            CALL ERREXIT(4)
-!         endif
+
          print*, "Read the grid from file "//trim(OUTGRID)
 
          allocate(tmpvar(nx+1,ny+1))
@@ -818,7 +726,6 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
             endif            
          endif
          
-         
          allocate(tmpvar(nx,ny))
          error=nf_inq_varid(ncid, 'area', id_var)
          call netcdf_err(error, 'inquire varid of area from file '
@@ -834,31 +741,10 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
                dy(i,j) = dx(i,j)
             enddo
          enddo
-!         allocate(tmpvar(nx,ny+1))
-
-!         error=nf_inq_varid(ncid, 'dx', id_var)
-!         call netcdf_err(error, 'inquire varid of dx from file '
-!     &                   //trim(OUTGRID) )
-!         error=nf_get_var_double(ncid, id_var, tmpvar)
-!         call netcdf_err(error, 'inquire data of dx from file '
-!     &                   //trim(OUTGRID) )
-!         dx(1:IM,1:JM+1) = tmpvar(2:nx:2,1:ny+1:2)
-!         deallocate(tmpvar)
-
-!          allocate(tmpvar(nx+1,ny))
-!         error=nf_inq_varid(ncid, 'dy', id_var)
-!         call netcdf_err(error, 'inquire varid of dy from file '
-!     &                   //trim(OUTGRID) )
-!         error=nf_get_var_double(ncid, id_var, tmpvar)
-!         call netcdf_err(error, 'inquire data of dy from file '
-!     &                   //trim(OUTGRID) )
-!         dy(1:IM+1,1:JM) = tmpvar(1:nx+1:2,2:ny:2)
-         deallocate(tmpvar)         
-      endif
+
       tend=timef()
       write(6,*)' Timer 1 time= ',tend-tbeg
                                 !
-      if(grid_from_file) then
        tbeg=timef()
 
          IF (MERGE_FILE == 'none') then
@@ -884,10 +770,6 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
 
       tend=timef()
       write(6,*)' MAKEMT2 time= ',tend-tbeg
-      else
-        CALL MAKEMT(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,GLAT,
-     &   IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
-      endif
 
        call minmxj(IM,JM,ORO,'     ORO')
        call minmxj(IM,JM,SLM,'     SLM')
@@ -914,16 +796,12 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
 C ===  Compute mtn principal coord HTENSR: THETA,GAMMA,SIGMA
 C
       allocate (THETA(IM,JM),GAMMA(IM,JM),SIGMA(IM,JM),ELVMAX(IM,JM))
-       if(grid_from_file) then       
+
       tbeg=timef()
          CALL MAKEPC2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,GLAT,
      1            IM,JM,IMN,JMN,geolon_c,geolat_c,SLM)
       tend=timef()
       write(6,*)' MAKEPC2 time= ',tend-tbeg
-       else
-         CALL MAKEPC(ZAVG,ZSLM,THETA,GAMMA,SIGMA,GLAT,
-     1            IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
-       endif
 
        call minmxj(IM,JM,THETA,'   THETA')
        call minmxj(IM,JM,GAMMA,'   GAMMA')
@@ -943,7 +821,7 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
 
        call minmxj(IM,JM,ORO,'     ORO')
        print*, "inputorog=", trim(INPUTOROG)
-       if(grid_from_file) then
+
          if(trim(INPUTOROG) == "none") then
            print*, "calling MAKEOA2 to compute OA, OL"
       tbeg=timef()
@@ -1057,21 +935,15 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
      &                   //trim(INPUTOROG) )
            
            print*, "calling MAKEOA3 to compute OA, OL"
-           CALL MAKEOA3(ZAVG,zslm,VAR,GLAT,OA,OL,IWORK,ELVMAX,ORO,SLM,
+           CALL MAKEOA3(ZAVG,VAR,GLAT,OA,OL,IWORK,ELVMAX,ORO,SLM,
      1            WORK1,WORK2,WORK3,WORK4,WORK5,WORK6,
      2            IM,JM,IMN,JMN,geolon_c,geolat_c,
-     3            geolon,geolat,is_south_pole,is_north_pole,nx_in,ny_in,
+     3            geolon,geolat,nx_in,ny_in,
      4            oa_in,ol_in,slm_in,lon_in,lat_in)
 
            deallocate(oa_in,ol_in,slm_in,lon_in,lat_in)
 
          endif
-       else
-         CALL MAKEOA(ZAVG,VAR,GLAT,OA,OL,IWORK,ELVMAX,ORO,
-     1            WORK1,WORK2,WORK3,WORK4,
-     2            WORK5,WORK6,
-     3            IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
-       endif
 
 !      Deallocate 2d vars
        deallocate(IST,IEN)
@@ -1366,67 +1238,23 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
       ENDDO
 !
        call mnmxja(IM,JM,ELVMAX,itest,jtest,'  ELVMAX')
-! --- Quadratic filter applied by default.
-! --- NF0 is normally set to an even value beyond the previous truncation, 
-! --- for example, for jcap=382, NF0=254+2
-! --- NF1 is set as jcap+2 (and/or nearest even), eg., for t382, NF1=382+2=384
-! --- if no filter is desired then NF1=NF0=0 and ORF=ORO
-! --- if no filter but spectral to grid (with gibbs) then NF1=jcap+2, and NF1=jcap+1
 !
       deallocate(VAR4)
       allocate (ORF(IM,JM))
-      IF ( NF1 - NF0 .eq. 0 ) FILTER=.FALSE.
-      print *,' NF1, NF0, FILTER=',NF1,NF0,FILTER
-      IF (FILTER) THEN
-C       SPECTRALLY TRUNCATE AND FILTER OROGRAPHY
-        do j=1,jm
-          if(numi(j).lt.im) then
-            ffj=cmplx(0.,0.)
-            call spfft1(numi(j),im/2+1,numi(j),1,ffj,oro(1,j),-1)
-            call spfft1(im,im/2+1,im,1,ffj,oro(1,j),+1)
-          endif
-        enddo
-        CALL SPTEZ(NR,NM,4,IM,JM,ORS,ORO,-1)
-!
-      print *,' about to apply spectral filter '
-        FFF=1./(NF1-NF0)**2
-        I=0
-        DO M=0,NM
-        DO N=M,NM+NR*M
-          IF(N.GT.NF0) THEN
-            WWW=MAX(1.-FFF*(N-NF0)**2,0.)
-            ORS(I+1)=ORS(I+1)*WWW
-            ORS(I+2)=ORS(I+2)*WWW
-          ENDIF
-        I=I+2
-        ENDDO
-        ENDDO
-!
-        CALL SPTEZ(NR,NM,4,IM,JM,ORS,ORF,+1)
-        do j=1,jm
-          if(numi(j).lt.im) then
-            call spfft1(im,im/2+1,im,1,ffj,orf(1,j),-1)
-            call spfft1(numi(j),im/2+1,numi(j),1,ffj,orf(1,j),+1)
-          endif
-        enddo
 
-      ELSE
-        IF (REVLAT) THEN
-          CALL REVERS(IM, JM, numi, SLM, WORK1)
-          CALL REVERS(IM, JM, numi, ORO, WORK1)
-          DO IMT=1,NMT
-            CALL REVERS(IM, JM, numi, HPRIME(1,1,IMT), WORK1)
-          ENDDO
-        ENDIF
-        ORS=0.
-        ORF=ORO
+      IF (REVLAT) THEN
+        CALL REVERS(IM, JM, SLM, WORK1)
+        CALL REVERS(IM, JM, ORO, WORK1)
+        DO IMT=1,NMT
+          CALL REVERS(IM, JM, HPRIME(1,1,IMT), WORK1)
+        ENDDO
       ENDIF
+      ORF=ORO
 
       deallocate (WORK1)
 
        call mnmxja(IM,JM,ELVMAX,itest,jtest,'  ELVMAX')
       print *,' ELVMAX(',itest,jtest,')=',ELVMAX(itest,jtest)
-      print *,' after spectral filter is applied'
        call minmxj(IM,JM,ORO,'     ORO')
        call minmxj(IM,JM,ORF,'     ORF')
 C
@@ -1459,152 +1287,14 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
       ENDDO
       tend=timef()
       write(6,*)' Timer 5 time= ',tend-tbeg
-      if (output_binary) then
-      tbeg=timef()
-C       OUTPUT BINARY FIELDS
-        print *,' OUTPUT BINARY FIELDS'
-        WRITE(51) REAL(SLM,4)
-        WRITE(52) REAL(ORF,4)
-        WRITE(53) REAL(HPRIME,4)
-        WRITE(54) REAL(ORS,4)
-        WRITE(55) REAL(ORO,4)
-        WRITE(66) REAL(THETA,4)
-        WRITE(67) REAL(GAMMA,4)
-        WRITE(68) REAL(SIGMA,4)
-! --- OCLSM is real(4) write only if ocean mask is present
-            if ( mskocn .eq. 1 ) then
-        ios=0
-         WRITE(27,iostat=ios) OCLSM
-         print *,' write OCLSM input:',ios
-!      print *,' LSM:',OCLSM(1,1),OCLSM(50,50),OCLSM(75,75),OCLSM(IM,JM)
-            endif
-C
-       call minmxj(IM,JM,ORO,'     ORO')
-      print *,' IM=',IM,' JM=',JM,' SPECTR=',SPECTR
-C---    Test binary file output:
-      WRITE(71) REAL(SLM,4)
-      DO IMT=1,NMT
-        WRITE(71) REAL(HPRIME(:,:,IMT),4)
-        print *,' HPRIME(',itest,jtest,imt,')=',HPRIME(itest,jtest,imt)
-      ENDDO
-      WRITE(71) REAL(ORO,4)
-      IF (SPECTR) THEN
-        WRITE(71) REAL(ORF,4)   ! smoothed spectral orography!
-      ENDIF
-C  OUTPUT GRIB FIELDS
-      KPDS=0
-      KPDS(1)=7
-      KPDS(2)=78
-      KPDS(3)=255
-      KPDS(4)=128
-      KPDS(5)=81
-      KPDS(6)=1
-      kpds(8)=2004
-      KPDS(9)=1
-      KPDS(10)=1
-      KPDS(13)=4
-      KPDS(15)=1
-      KPDS(16)=51
-      KPDS(17)=1
-      KPDS(18)=1
-      KPDS(19)=1
-      KPDS(21)=20
-      KPDS(22)=0
-      KGDS=0
-      KGDS(1)=4
-      KGDS(2)=IM
-      KGDS(3)=JM
-      KGDS(4)=90000-180000/PI*RCLT(1)
-      KGDS(6)=128
-      KGDS(7)=180000/PI*RCLT(1)-90000
-      KGDS(8)=-NINT(360000./IM)
-      KGDS(9)=NINT(360000./IM)
-      KGDS(10)=JM/2
-      KGDS(20)=255
-! --- SLM
-      CALL BAOPEN(56,'fort.56',IRET)
-      if (iret .ne. 0) print *,' BAOPEN ERROR UNIT 56: IRET=',IRET
-      CALL PUTGB(56,IM*JM,KPDS,KGDS,LB,SLM,IRET)
-      print *,' SLM: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-      if (iret .ne. 0) print *,' SLM PUTGB ERROR:  UNIT 56: IRET=',IRET
-      print *,' SLM: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-! --- OCLSM if present
-!           if ( mskocn .eq. 1 ) then
-!     CALL BAOPEN(27,'fort.27',IRET)
-!     if (iret .ne. 0) print *,' OCLSM BAOPEN ERROR UNIT 27:IRET=',IRET
-!     CALL PUTGB(27,IM*JM,KPDS,KGDS,LB,OCLSM,IRET)
-!     if (iret .ne. 0) print *,' OCLSM PUTGB ERROR: UNIT 27:IRET=',IRET
-!     print *,' OCLSM: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-!           endif
-
-      KPDS(5)=8
-      IF (SPECTR) THEN
-        CALL BAOPEN(57,'fort.57',IRET)
-        CALL PUTGB(57,IM*JM,KPDS,KGDS,LB,ORF,IRET)
-      print *,' ORF (ORO): putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-      ENDIF
-C
-C ===  write out theta (angle of land to East) using #101 (wave dir)
-C ===  [radians] and since < 1 scale adjust kpds(22)
-C
-      KPDS(5)=101
-        CALL BAOPEN(58,'fort.58',IRET)
-        CALL PUTGB(58,IM*JM,KPDS,KGDS,LB,THETA,IRET)
-      print *,' THETA: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-C
-C ===  write out (land aspect ratio or anisotropy)  using #102 
-C ===  (as in wind wave hgt)
-C
-      KPDS(22)=2
-      KPDS(5)=102
-        CALL BAOPEN(60,'fort.60',IRET)
-        CALL PUTGB(60,IM*JM,KPDS,KGDS,LB,SIGMA,IRET)
-      print *,' SIGMA: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-C
-C ===  write out (slope parameter sigma)  using #9 
-C ===  (as in std hgt)
-C
-      KPDS(22)=1
-      KPDS(5)=103
-        CALL BAOPEN(59,'fort.59',IRET)
-        CALL PUTGB(59,IM*JM,KPDS,KGDS,LB,GAMMA,IRET)
-      print *,' GAMMA: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-C
-      KPDS(22)=1
-      KPDS(5)=9
-        CALL BAOPEN(61,'fort.61',IRET)
-        CALL PUTGB(61,IM*JM,KPDS,KGDS,LB,HPRIME,IRET)
-      print *,' HPRIME: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-C
-C
-      KPDS(22)=0
-      KPDS(5)=8
-        CALL BAOPEN(62,'fort.62',IRET)
-        CALL PUTGB(62,IM*JM,KPDS,KGDS,LB,ELVMAX,IRET)
-      print *,' ELVMAX: putgb-KPDS(22,5),iret:',KPDS(22),KPDS(5),IRET
-      endif ! output_binary
-C
-      DELXN = 360./IM
-      do i=1,im
-        xlon(i) = DELXN*(i-1)
+
+      do j = 1, jm
+         xlat(j) = geolat(1,j)
       enddo
-      IF(trim(OUTGRID) == "none") THEN
-        do j=1,jm
-         do i=1,im
-           geolon(i,j) = xlon(i)
-           geolat(i,j) = xlat(j)
-         enddo
-        enddo
-      else
-        do j = 1, jm
-           xlat(j) = geolat(1,j)
-        enddo
-        do i = 1, im
-           xlon(i) = geolon(i,1)
-        enddo
-      endif
-      tend=timef()
-      write(6,*)' Binary output time= ',tend-tbeg
+      do i = 1, im
+         xlon(i) = geolon(i,1)
+      enddo
+
       tbeg=timef()
       CALL WRITE_NETCDF(IM,JM,SLM,land_frac,ORO,ORF,HPRIME,1,1,
      1                  GEOLON(1:IM,1:JM),GEOLAT(1:IM,1:JM), XLON,XLAT)
@@ -1615,8 +1305,8 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
       print *,' ===== Deallocate Arrays and ENDING MTN VAR OROG program'
 
 !     Deallocate 1d vars
-      deallocate(JST,JEN,KPDS,KGDS,numi,lonsperlat)
-      deallocate(COSCLT,WGTCLT,RCLT,XLAT,DIFFX,XLON,ORS,oaa,ola,GLAT)
+      deallocate(numi,lonsperlat)
+      deallocate(WGTCLT,XLAT,XLON,oaa,ola,GLAT)
 
 !     Deallocate 2d vars
       deallocate (OCLSM)
@@ -1628,183 +1318,7 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
       tend=timef()
       write(6,*)' Total runtime time= ',tend-tbeg1
       RETURN
-      END
-
-!> Create the orography, land-mask, standard deviation of
-!! orography and the convexity on a model gaussian grid.
-!! This routine was used for the spectral GFS model.
-!!
-!! @param[in] zavg The high-resolution input orography dataset.
-!! @param[in] zslm The high-resolution input land-mask dataset.
-!! @param[out] oro Orography on the model grid.
-!! @param[out] slm Land-mask on the model grid.
-!! @param[out] var Standard deviation of orography on the model grid.
-!! @param[out] var4 Convexity on the model grid.
-!! @param[out] glat Latitude of each row of the high-resolution 
-!! orography and land-mask datasets.
-!! @param[out] ist This is the 'i' index of high-resolution data set
-!! at the east edge of the model grid cell.
-!! the high-resolution dataset with respect to the 'east' edge
-!! @param[out] ien This is the 'i' index of high-resolution data set
-!! at the west edge of the model grid cell.
-!! @param[out] jst  This is the 'j' index of high-resolution data set
-!! at the south edge of the model grid cell.
-!! @param[out] jen This is the 'j' index of high-resolution data set
-!! at the north edge of the model grid cell.
-!! @param[in] im "i" dimension of the model grid.
-!! @param[in] jm "j" dimension of the model grid.
-!! @param[in] imn "i" dimension of the hi-res input orog/mask dataset.
-!! @param[in] jmn "j" dimension of the hi-res input orog/mask dataset.
-!! @param[in] xlat The latitude of each row of the model grid.
-!! @param[in] numi For reduced gaussian grids, the number of 'i' points
-!! for each 'j' row.
-!! @author Jordan Alpert NOAA/EMC
-      SUBROUTINE MAKEMT(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
-     1 GLAT,IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
-      DIMENSION GLAT(JMN),XLAT(JM)
-!     REAL*4 OCLSM
-      INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
-      DIMENSION ORO(IM,JM),SLM(IM,JM),VAR(IM,JM),VAR4(IM,JM)
-      DIMENSION IST(IM,jm),IEN(IM,jm),JST(JM),JEN(JM),numi(jm)
-      INTEGER mskocn,isave
-      LOGICAL FLAG, DEBUG
-C==== DATA DEBUG/.TRUE./
-      DATA DEBUG/.FALSE./
-C
-! ---- OCLSM holds the ocean (im,jm) grid
-! ---  mskocn=1 Use ocean model sea land mask, OK and present,
-! ---  mskocn=0 dont use Ocean model sea land mask, not OK, not present
-      print *,' _____ SUBROUTINE MAKEMT '
-C---- GLOBAL XLAT AND XLON ( DEGREE )
-C
-      JM1 = JM - 1
-      DELXN = 360./IMN      ! MOUNTAIN DATA RESOLUTION
-C
-      DO J=1,JMN
-         GLAT(J) = -90. + (J-1) * DELXN + DELXN * 0.5
-      ENDDO
-C
-C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
-C
-C  (*j*)  for hard wired zero offset (lambda s =0) for terr05 
-      DO J=1,JM
-       DO I=1,numi(j)
-         IM1 = numi(j) - 1
-         DELX  = 360./numi(j)       ! GAUSSIAN GRID RESOLUTION
-         FACLON  = DELX / DELXN
-         IST(I,j) = FACLON * FLOAT(I-1) - FACLON * 0.5 + 1
-         IEN(I,j) = FACLON * FLOAT(I) - FACLON * 0.5 + 1
-!        IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
-!        IEN(I,j) = FACLON * FLOAT(I)   + 0.0001
-C
-         IF (IST(I,j) .LE. 0)      IST(I,j) = IST(I,j) + IMN
-         IF (IEN(I,j) .LT. IST(I,j)) IEN(I,j) = IEN(I,j) + IMN
-!
-!          if ( I .lt. 10  .and. J .ge. JM-1 )
-!    1   PRINT*,' MAKEMT: I j IST IEN ',I,j,IST(I,j),IEN(I,j)
-       ENDDO
-!          if ( J .ge. JM-1 ) then
-!     print *,' *** FACLON=',FACLON, 'numi(j=',j,')=',numi(j)
-!          endif
-      ENDDO
-      print *,' DELX=',DELX,' DELXN=',DELXN
-      DO J=1,JM-1
-         FLAG=.TRUE.
-         DO J1=1,JMN
-            XXLAT = (XLAT(J)+XLAT(J+1))/2.
-            IF(FLAG.AND.GLAT(J1).GT.XXLAT) THEN
-              JST(J) = J1
-              JEN(J+1) = J1 - 1
-              FLAG = .FALSE.
-            ENDIF
-         ENDDO
-CX     PRINT*, ' J JST JEN ',J,JST(J),JEN(J),XLAT(J),GLAT(J1)
-      ENDDO
-      JST(JM) = MAX(JST(JM-1) - (JEN(JM-1)-JST(JM-1)),1)
-      JEN(1)  = MIN(JEN(2) + (JEN(2)-JST(2)),JMN)      
-!      PRINT*, ' JM JST JEN=',JST(JM),JEN(JM),XLAT(JM),GLAT(JMN)
-C
-C...FIRST, AVERAGED HEIGHT
-C
-      DO J=1,JM
-        DO I=1,numi(j)
-            ORO(I,J)  = 0.0
-            VAR(I,J)  = 0.0
-            VAR4(I,J) = 0.0
-            XNSUM = 0.0
-            XLAND = 0.0
-            XWATR = 0.0
-            XL1 = 0.0
-            XS1 = 0.0
-            XW1 = 0.0
-            XW2 = 0.0
-            XW4 = 0.0
-            DO II1 = 1, IEN(I,J) - IST(I,J) + 1
-               I1 = IST(I,J) + II1 - 1
-               IF(I1.LE.0.)  I1 = I1 + IMN
-               IF(I1.GT.IMN) I1 = I1 - IMN
-!        if ( i .le. 10 .and. i .ge. 1 ) then 
-!             if (j .eq. JM )
-!    &print *,' J,JST,JEN,IST,IEN,I1=',
-!    &J,JST(j),JEN(J),IST(I,j),IEN(I,j),I1
-!        endif
-               DO J1=JST(J),JEN(J)
-                  XLAND = XLAND + FLOAT(ZSLM(I1,J1))
-                  XWATR = XWATR + FLOAT(1-ZSLM(I1,J1))
-                  XNSUM = XNSUM + 1.
-                  HEIGHT = FLOAT(ZAVG(I1,J1)) 
-C.........
-                  IF(HEIGHT.LT.-990.) HEIGHT = 0.0
-                  XL1 = XL1 + HEIGHT * FLOAT(ZSLM(I1,J1))
-                  XS1 = XS1 + HEIGHT * FLOAT(1-ZSLM(I1,J1))
-                  XW1 = XW1 + HEIGHT
-                  XW2 = XW2 + HEIGHT ** 2
-C   check antarctic pole 
-!        if ( i .le. 10 .and. i .ge. 1 )then
-!           if (j .ge. JM-1 )then
-C=== degub testing
-!     print *," I,J,I1,J1,XL1,XS1,XW1,XW2:",I,J,I1,J1,XL1,XS1,XW1,XW2
-! 153 format(1x,' ORO,ELVMAX(i=',i4,' j=',i4,')=',2E14.5,3f5.1)
-!          endif
-!        endif
-               ENDDO
-            ENDDO
-            IF(XNSUM.GT.1.) THEN
-! --- SLM initialized with OCLSM calc from all land points except ....
-! ---  0 is ocean and 1 is land for slm
-! ---  Step 1 is to only change SLM after GFS SLM is applied
-                 
-               SLM(I,J) = FLOAT(NINT(XLAND/XNSUM))
-               IF(SLM(I,J).NE.0.) THEN
-                  ORO(I,J)= XL1 / XLAND
-               ELSE
-                  ORO(I,J)= XS1 / XWATR
-               ENDIF
-               VAR(I,J)=SQRT(MAX(XW2/XNSUM-(XW1/XNSUM)**2,0.))
-            DO II1 = 1, IEN(I,j) - IST(I,J) + 1
-               I1 = IST(I,J) + II1 - 1
-                  IF(I1.LE.0.) I1 = I1 + IMN
-                  IF(I1.GT.IMN) I1 = I1 - IMN
-                  DO J1=JST(J),JEN(J)
-                     HEIGHT = FLOAT(ZAVG(I1,J1)) 
-                     IF(HEIGHT.LT.-990.) HEIGHT = 0.0
-                     XW4 = XW4 + (HEIGHT-ORO(I,J)) ** 4
-                  ENDDO
-               ENDDO
-               IF(VAR(I,J).GT.1.) THEN
-!          if ( I .lt. 20  .and. J .ge. JM-19 ) then
-!     print *,'I,J,XW4,XNSUM,VAR(I,J)',I,J,XW4,XNSUM,VAR(I,J)
-!          endif
-                  VAR4(I,J) = MIN(XW4/XNSUM/VAR(I,J) **4,10.)
-               ENDIF
-            ENDIF
-         ENDDO
-      ENDDO
-      WRITE(6,*) "! MAKEMT ORO SLM VAR VAR4 DONE"
-C
-
-      RETURN
-      END
+      END SUBROUTINE TERSUB
 
 !> Determine the location of a cubed-sphere point within
 !! the high-resolution orography data.  The location is
@@ -2250,290 +1764,7 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
       deallocate(hgt_1d)
       deallocate(hgt_1d_all)
       RETURN
-      END
-
-!> Make the principle coordinates - slope of orography, 
-!! anisotropy, angle of mountain range with respect to east.
-!! This routine is used for spectral GFS gaussian grids.
-!!
-!! @param[in] zavg The high-resolution input orography dataset.
-!! @param[in] zslm The high-resolution input land-mask dataset.
-!! @param[out] theta Angle of mountain range with respect to
-!! east for each model point.
-!! @param[out] gamma Anisotropy for each model point.
-!! @param[out] sigma Slope of orography for each model point.
-!! @param[out] glat Latitude of each row of the high-resolution 
-!! orography and land-mask datasets.
-!! @param[out] ist This is the 'i' index of high-resolution data set
-!! at the east edge of the model grid cell.
-!! @param[out] ien This is the 'i' index of high-resolution data set
-!! at the west edge of the model grid cell.
-!! @param[out] jst  This is the 'j' index of high-resolution data set
-!! at the south edge of the model grid cell.
-!! @param[out] jen This is the 'j' index of high-resolution data set
-!! at the north edge of the model grid cell.
-!! @param[in] im "i" dimension of the model grid tile.
-!! @param[in] jm "j" dimension of the model grid tile.
-!! @param[in] imn "i" dimension of the hi-res input orog/mask datasets.
-!! @param[in] jmn "j" dimension of the hi-res input orog/mask datasets.
-!! @param[in] xlat The latitude of each row of the model grid.
-!! @param[in] numi For reduced gaussian grids, the number of 'i' points
-!! for each 'j' row.
-!! @author Jordan Alpert NOAA/EMC
-      SUBROUTINE MAKEPC(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
-     1           GLAT,IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
-C
-C===  PC: principal coordinates of each Z avg orog box for L&M
-C
-      parameter(REARTH=6.3712E+6)
-      DIMENSION GLAT(JMN),XLAT(JM),DELTAX(JMN)
-      INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
-      DIMENSION ORO(IM,JM),SLM(IM,JM),HL(IM,JM),HK(IM,JM)
-      DIMENSION HX2(IM,JM),HY2(IM,JM),HXY(IM,JM),HLPRIM(IM,JM)
-      DIMENSION THETA(IM,JM),GAMMA(IM,JM),SIGMA2(IM,JM),SIGMA(IM,JM)
-      DIMENSION IST(IM,jm),IEN(IM,jm),JST(JM),JEN(JM),numi(jm)
-      LOGICAL FLAG, DEBUG
-C===  DATA DEBUG/.TRUE./
-      DATA DEBUG/.FALSE./
-C
-      PI = 4.0 * ATAN(1.0)
-      CERTH = PI * REARTH
-C---- GLOBAL XLAT AND XLON ( DEGREE )
-C
-      JM1 = JM - 1
-      DELXN = 360./IMN      ! MOUNTAIN DATA RESOLUTION
-      DELTAY =  CERTH / FLOAT(JMN)
-      print *, 'MAKEPC: DELTAY=',DELTAY
-C
-      DO J=1,JMN
-         GLAT(J) = -90. + (J-1) * DELXN + DELXN * 0.5
-         DELTAX(J) = DELTAY * COS(GLAT(J)*PI/180.0)
-      ENDDO
-C
-C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
-C
-      DO J=1,JM
-      DO I=1,numi(j)
-C        IM1 = numi(j) - 1
-         DELX  = 360./numi(j)       ! GAUSSIAN GRID RESOLUTION
-         FACLON  = DELX / DELXN
-         IST(I,j) = FACLON * FLOAT(I-1) - FACLON * 0.5
-         IST(I,j) = IST(I,j) + 1
-         IEN(I,j) = FACLON * FLOAT(I) - FACLON * 0.5
-C              if (debug) then
-C          if ( I .lt. 10 .and. J .lt. 10 )
-C    1     PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
-C              endif
-!        IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
-!        IEN(I,j) = FACLON * FLOAT(I)   + 0.0001
-         IF (IST(I,j) .LE. 0)      IST(I,j) = IST(I,j) + IMN
-         IF (IEN(I,j) .LT. IST(I,j)) IEN(I,j) = IEN(I,j) + IMN
-               if (debug) then
-           if ( I .lt. 10 .and. J .lt. 10 )
-     1     PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
-               endif
-         IF (IEN(I,j) .LT. IST(I,j)) 
-     1 print *,' MAKEPC: IEN < IST: I,J,IST(I,J),IEN(I,J)',
-     2 I,J,IST(I,J),IEN(I,J)
-      ENDDO
-      ENDDO
-      DO J=1,JM-1
-         FLAG=.TRUE.
-         DO J1=1,JMN
-            XXLAT = (XLAT(J)+XLAT(J+1))/2.
-            IF(FLAG.AND.GLAT(J1).GT.XXLAT) THEN
-              JST(J) = J1
-              JEN(J+1) = J1 - 1
-              FLAG = .FALSE.
-            ENDIF
-         ENDDO
-      ENDDO
-      JST(JM) = MAX(JST(JM-1) - (JEN(JM-1)-JST(JM-1)),1)
-      JEN(1)  = MIN(JEN(2) + (JEN(2)-JST(2)),JMN)
-               if (debug) then
-        PRINT*, ' IST,IEN(1,1-numi(1,JM))',IST(1,1),IEN(1,1), 
-     1  IST(numi(JM),JM),IEN(numi(JM),JM), numi(JM)
-        PRINT*, ' JST,JEN(1,JM) ',JST(1),JEN(1),JST(JM),JEN(JM) 
-                endif
-C
-C... DERIVITIVE TENSOR OF HEIGHT
-C
-      DO J=1,JM
-        DO I=1,numi(j)
-            ORO(I,J)  = 0.0
-            HX2(I,J) = 0.0
-            HY2(I,J) = 0.0
-            HXY(I,J) = 0.0
-            XNSUM = 0.0
-            XLAND = 0.0
-            XWATR = 0.0
-            XL1 = 0.0
-            XS1 = 0.0
-            xfp = 0.0
-            yfp = 0.0
-            xfpyfp = 0.0
-            xfp2 = 0.0
-            yfp2 = 0.0
-            HL(I,J) = 0.0
-            HK(I,J) = 0.0
-            HLPRIM(I,J) = 0.0
-            THETA(I,J) = 0.0 
-            GAMMA(I,J) = 0.
-            SIGMA2(I,J) = 0.
-            SIGMA(I,J) = 0.
-C
-            DO II1 = 1, IEN(I,J) - IST(I,J) + 1
-               I1 = IST(I,J) + II1 - 1
-               IF(I1.LE.0.)  I1 = I1 + IMN
-               IF(I1.GT.IMN) I1 = I1 - IMN
-C
-C===  set the rest of the indexs for ave: 2pt staggered derivitive
-C
-                i0 = i1 - 1
-                if (i1 - 1 .le. 0 )   i0 = i0 + imn
-                if (i1 - 1 .gt. imn)  i0 = i0 - imn
-C
-                ip1 = i1 + 1
-                if (i1 + 1 .le. 0 )   ip1 = ip1 + imn
-                if (i1 + 1 .gt. imn)  ip1 = ip1 - imn
-C
-               DO J1=JST(J),JEN(J)
-                if (debug) then
-                   if ( I1 .eq. IST(I,J) .and. J1 .eq. JST(J) ) 
-     1   PRINT*, ' J, J1,IST,JST,DELTAX,GLAT ',
-     2             J,J1,IST(I,J),JST(J),DELTAX(J1),GLAT(J1)  
-                   if ( I1 .eq. IEN(I,J) .and. J1 .eq. JEN(J) ) 
-     1   PRINT*, ' J, J1,IEN,JEN,DELTAX,GLAT ',
-     2             J,J1,IEN(I,J),JEN(J),DELTAX(J1),GLAT(J1)  
-                endif
-                  XLAND = XLAND + FLOAT(ZSLM(I1,J1))
-                  XWATR = XWATR + FLOAT(1-ZSLM(I1,J1))
-                  XNSUM = XNSUM + 1.
-C
-                  HEIGHT = FLOAT(ZAVG(I1,J1))
-                  hi0 =  float(zavg(i0,j1))
-                  hip1 =  float(zavg(ip1,j1))
-C
-                  IF(HEIGHT.LT.-990.) HEIGHT = 0.0
-                  if(hi0 .lt. -990.)  hi0 = 0.0
-                  if(hip1 .lt. -990.)  hip1 = 0.0
-C........           xfp = xfp + 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
-           xfp = 0.5 * ( hip1 - hi0 ) / DELTAX(J1)
-          xfp2 = xfp2 + 0.25 * ( ( hip1 - hi0 )/DELTAX(J1) )** 2 
-C
-! --- not at boundaries
-!RAB                 if ( J1 .ne. JST(1)  .and. J1 .ne. JEN(JM) ) then
-                 if ( J1 .ne. JST(JM)  .and. J1 .ne. JEN(1) ) then
-                  hj0 =  float(zavg(i1,j1-1))
-                  hjp1 =  float(zavg(i1,j1+1))
-                  if(hj0 .lt. -990.)  hj0 = 0.0
-                  if(hjp1 .lt. -990.)  hjp1 = 0.0
-C.......          yfp = yfp + 0.5 * ( hjp1 - hj0 ) / DELTAY
-                  yfp = 0.5 * ( hjp1 - hj0 ) / DELTAY
-                  yfp2 = yfp2 + 0.25 * ( ( hjp1 - hj0 )/DELTAY )**2   
-C
-C..............elseif ( J1 .eq. JST(J) .or. J1 .eq. JEN(JM) ) then
-C ===     the NH pole: NB J1 goes from High at NP to Low toward SP
-C
-!RAB                 elseif ( J1 .eq. JST(1) ) then
-                 elseif ( J1 .eq. JST(JM) ) then
-		 ijax = i1 + imn/2 
-               if (ijax .le. 0 )   ijax = ijax + imn
-               if (ijax .gt. imn)  ijax = ijax - imn
-C..... at N pole we stay at the same latitude j1 but cross to opp side
-                 hijax = float(zavg(ijax,j1))
-                 hi1j1 = float(zavg(i1,j1))
-                  if(hijax .lt. -990.)  hijax = 0.0
-                  if(hi1j1 .lt. -990.)  hi1j1 = 0.0
-C.......        yfp = yfp + 0.5 * ( ( 0.5 * ( hijax + hi1j1) ) - hi1j1 )/DELTAY
-        yfp = 0.5 * ( ( 0.5 * ( hijax - hi1j1 ) ) )/DELTAY
-        yfp2 = yfp2 + 0.25 * ( ( 0.5 *  ( hijax - hi1j1) ) 
-     1                                              / DELTAY )**2
-C
-C ===     the SH pole: NB J1 goes from High at NP to Low toward SP
-C
-!RAB                 elseif ( J1 .eq. JEN(JM) ) then
-                 elseif ( J1 .eq. JEN(1) ) then
-		 ijax = i1 + imn/2 
-               if (ijax .le. 0 )   ijax = ijax + imn
-               if (ijax .gt. imn)  ijax = ijax - imn
-                 hijax = float(zavg(ijax,j1))
-                 hi1j1 = float(zavg(i1,j1))
-                  if(hijax  .lt. -990.)  hijax = 0.0
-                  if(hi1j1  .lt. -990.)  hi1j1 = 0.0
-             if ( i1 .lt. 5 )print *,' S.Pole i1,j1 :',i1,j1,hijax,hi1j1
-C.....        yfp = yfp + 0.5 *  (0.5 * ( hijax - hi1j1) )/DELTAY  
-        yfp = 0.5 *  (0.5 * ( hijax - hi1j1) )/DELTAY  
-        yfp2 = yfp2 + 0.25 * (  (0.5 * (hijax - hi1j1) )
-     1                                                 / DELTAY )**2  
-                 endif
-C
-C ===    The above does an average across the pole for the bndry in j.
-C23456789012345678901234567890123456789012345678901234567890123456789012......
-C
-                  xfpyfp = xfpyfp + xfp * yfp
-                  XL1 = XL1 + HEIGHT * FLOAT(ZSLM(I1,J1))
-                  XS1 = XS1 + HEIGHT * FLOAT(1-ZSLM(I1,J1))
-C
-C === average the HX2, HY2 and HXY
-C === This will be done over all land
-C
-               ENDDO
-            ENDDO
-C
-C ===  HTENSR 
-C
-           IF(XNSUM.GT.1.) THEN
-               SLM(I,J) = FLOAT(NINT(XLAND/XNSUM))
-               IF(SLM(I,J).NE.0.) THEN
-                  ORO(I,J)= XL1 / XLAND
-                  HX2(I,J) =  xfp2  / XLAND
-                  HY2(I,J) =  yfp2  / XLAND
-		  HXY(I,J) =  xfpyfp / XLAND
-               ELSE
-                  ORO(I,J)= XS1 / XWATR
-               ENDIF
-C=== degub testing
-      if (debug) then
-          print *," I,J,i1,j1,HEIGHT:", I,J,i1,j1,HEIGHT,
-     1         XLAND,SLM(i,j)
-          print *," xfpyfp,xfp2,yfp2:",xfpyfp,xfp2,yfp2
-          print *," HX2,HY2,HXY:",HX2(I,J),HY2(I,J),HXY(I,J)
-      ENDIF
-C
-C === make the principal axes, theta, and the degree of anisotropy, 
-C === and sigma2, the slope parameter
-C
-               HK(I,J) = 0.5 * ( HX2(I,J) + HY2(I,J) )
-               HL(I,J) = 0.5 * ( HX2(I,J) - HY2(I,J) )
-               HLPRIM(I,J) = SQRT(HL(I,J)*HL(I,J) + HXY(I,J)*HXY(I,J))
-           IF( HL(I,J).NE. 0. .AND. SLM(I,J) .NE. 0. ) THEN
-C
-             THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J)) * 180.0 / PI
-C ===   for testing print out in degrees
-C            THETA(I,J) = 0.5 * ATAN2(HXY(I,J),HL(I,J))
-            ENDIF
-             SIGMA2(I,J) =  ( HK(I,J) + HLPRIM(I,J) )
-        if ( SIGMA2(I,J) .GE. 0. ) then 
-             SIGMA(I,J) =  SQRT(SIGMA2(I,J) )
-             if (sigma2(i,j) .ne. 0. .and. 
-     &        HK(I,J) .GE. HLPRIM(I,J) ) 
-     1       GAMMA(I,J) = sqrt( (HK(I,J) - HLPRIM(I,J)) / SIGMA2(I,J) )
-        else
-             SIGMA(I,J)=0.
-        endif
-           ENDIF
-                  if (debug) then
-       print *," I,J,THETA,SIGMA,GAMMA,",I,J,THETA(I,J),
-     1                                       SIGMA(I,J),GAMMA(I,J)
-       print *," HK,HL,HLPRIM:",HK(I,J),HL(I,J),HLPRIM(I,J)
-                  endif
-        ENDDO
-      ENDDO
-      WRITE(6,*) "! MAKE Principal Coord  DONE"
-C
-      RETURN
-      END
+      END SUBROUTINE MAKEMT2
 
 !> Make the principle coordinates - slope of orography, 
 !! anisotropy, angle of mountain range with respect to east.
@@ -2578,7 +1809,7 @@ SUBROUTINE MAKEPC2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
       integer i,j,i1,j1,i2,jst,jen,numx,i0,ip1,ijax
       integer ilist(IMN)
       logical inside_a_polygon
-      LOGICAL FLAG, DEBUG
+      LOGICAL DEBUG
 C===  DATA DEBUG/.TRUE./
       DATA DEBUG/.FALSE./
 C
@@ -2782,366 +2013,6 @@ SUBROUTINE MAKEPC2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
       RETURN
       END SUBROUTINE MAKEPC2
       
-!> Create orographic asymmetry and orographic length scale on
-!! the model grid.  This routine is used for the spectral
-!! GFS gaussian grid.
-!!
-!! @param[in] zavg The high-resolution input orography dataset.
-!! @param[in] var Standard deviation of orography on the model grid.
-!! @param[out] glat Latitude of each row of input terrain dataset.
-!! @param[out] oa4 Orographic asymmetry on the model grid. Four
-!! directional components - W/S/SW/NW
-!! @param[out] ol Orographic length scale on the model grid. Four
-!! directional components - W/S/SW/NW
-!! @param[out] ioa4 Count of oa4 values between certain thresholds.
-!! @param[out] elvmax Maximum elevation on the model grid.
-!! @param[in] oro Orography on the model grid.
-!! @param[out] oro1 Save array for model grid orography.
-!! @param[out] xnsum Number of high-resolution orography points
-!! higher than the model grid box average.
-!! @param[out] xnsum1 Number of high-resolution orography points
-!! higher than the critical height.
-!! @param[out] xnsum2 Total number of high-resolution orography points
-!! within a model grid box.
-!! @param[out] xnsum3 Same as xnsum1, except shifted by half a 
-!! model grid box.
-!! @param[out] xnsum4 Same as xnsum2, except shifted by half a
-!! model grid box.
-!! @param[out] ist This is the 'i' index of high-resolution data set
-!! at the east edge of the model grid cell.
-!! @param[out] ien This is the 'i' index of high-resolution data set
-!! at the west edge of the model grid cell.
-!! @param[out] jst  This is the 'j' index of high-resolution data set
-!! at the south edge of the model grid cell.
-!! @param[out] jen This is the 'j' index of high-resolution data set
-!! at the north edge of the model grid cell.
-!! @param[in] im "i" dimension of the model grid.
-!! @param[in] jm "j" dimension of the model grid.
-!! @param[in] imn "i" dimension of the input terrain dataset.
-!! @param[in] jmn "j" dimension of the input terrain dataset.
-!! @param[in] xlat The latitude of each row of the model grid.
-!! @param[in] numi For reduced gaussian grids, the number of 'i' points
-!! for each 'j' row.
-!! @author Jordan Alpert NOAA/EMC
-      SUBROUTINE MAKEOA(ZAVG,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
-     1           ORO,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
-     2           IST,IEN,JST,JEN,IM,JM,IMN,JMN,XLAT,numi)
-      DIMENSION GLAT(JMN),XLAT(JM)
-      INTEGER ZAVG(IMN,JMN)
-      DIMENSION ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
-      DIMENSION OA4(IM,JM,4),IOA4(IM,JM,4)
-      DIMENSION IST(IM,jm),IEN(IM,jm),JST(JM),JEN(JM)
-      DIMENSION XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
-      DIMENSION XNSUM3(IM,JM),XNSUM4(IM,JM)
-      DIMENSION VAR(IM,JM),OL(IM,JM,4),numi(jm)
-      LOGICAL FLAG
-C
-C---- GLOBAL XLAT AND XLON ( DEGREE )
-C
-! --- IM1 = IM - 1 removed (not used in this sub)
-      JM1 = JM - 1
-      DELXN = 360./IMN      ! MOUNTAIN DATA RESOLUTION
-C
-      DO J=1,JMN
-         GLAT(J) = -90. + (J-1) * DELXN + DELXN * 0.5
-      ENDDO
-      print *,' IM=',IM,' JM=',JM,' IMN=',IMN,' JMN=',JMN
-C
-C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
-C
-      DO j=1,jm
-      DO I=1,numi(j)
-         DELX  = 360./numi(j)       ! GAUSSIAN GRID RESOLUTION
-         FACLON  = DELX / DELXN
-C ---  minus sign here in IST and IEN as in MAKEMT!
-         IST(I,j) = FACLON * FLOAT(I-1) - FACLON * 0.5
-         IST(I,j) = IST(I,j) + 1
-         IEN(I,j) = FACLON * FLOAT(I) - FACLON * 0.5
-!        IST(I,j) = FACLON * FLOAT(I-1) + 1.0001
-!        IEN(I,j) = FACLON * FLOAT(I)   + 0.0001
-         IF (IST(I,j) .LE. 0)      IST(I,j) = IST(I,j) + IMN
-         IF (IEN(I,j) .LT. IST(I,j)) IEN(I,j) = IEN(I,j) + IMN
-cx         PRINT*, ' I j IST IEN ',I,j,IST(I,j),IEN(I,j)
-           if ( I .lt. 3  .and. J .lt. 3 )
-     1PRINT*,' MAKEOA: I j IST IEN ',I,j,IST(I,j),IEN(I,j)
-           if ( I .lt. 3  .and. J .ge. JM-1 )
-     1PRINT*,' MAKEOA: I j IST IEN ',I,j,IST(I,j),IEN(I,j)
-      ENDDO
-      ENDDO
-      print *,'MAKEOA: DELXN,DELX,FACLON',DELXN,DELX,FACLON
-      print *, '  ***** ready to start JST JEN section '
-      DO J=1,JM-1
-         FLAG=.TRUE.
-         DO J1=1,JMN
-! --- XXLAT added as in MAKEMT and in next line as well
-            XXLAT = (XLAT(J)+XLAT(J+1))/2.
-            IF(FLAG.AND.GLAT(J1).GT.XXLAT) THEN
-              JST(J) = J1
-! ---         JEN(J+1) = J1 - 1
-              FLAG = .FALSE.
-           if ( J .eq. 1 )
-     1PRINT*,' MAKEOA: XX j JST JEN ',j,JST(j),JEN(j)
-            ENDIF
-         ENDDO
-           if ( J .lt. 3 )
-     1PRINT*,' MAKEOA: j JST JEN ',j,JST(j),JEN(j)
-           if ( J .ge. JM-2 )
-     1PRINT*,' MAKEOA: j JST JEN ',j,JST(j),JEN(j)
-C        FLAG=.TRUE.
-C        DO J1=JST(J),JMN
-C           IF(FLAG.AND.GLAT(J1).GT.XLAT(J)) THEN
-C             JEN(J) = J1 - 1
-C             FLAG = .FALSE.
-C           ENDIF
-C        ENDDO
-      ENDDO
-      JST(JM) = MAX(JST(JM-1) - (JEN(JM-1)-JST(JM-1)),1)
-      JEN(1)  = MIN(JEN(2) + (JEN(2)-JST(2)),JMN)
-      print *,' ***** JST(1) JEN(1) ',JST(1),JEN(1)
-      print *,' ***** JST(JM) JEN(JM) ',JST(JM),JEN(JM)
-C
-      DO J=1,JM
-        DO I=1,numi(j)
-          XNSUM(I,J) = 0.0
-          ELVMAX(I,J) = ORO(I,J)
-          ZMAX(I,J)   = 0.0
-        ENDDO
-      ENDDO
-!
-! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
-! ---  to JM or to JM1
-      DO J=1,JM
-        DO I=1,numi(j)
-            DO II1 = 1, IEN(I,J) - IST(I,J) + 1
-               I1 = IST(I,J) + II1 - 1
-! --- next line as in makemt (I1 not II1) (*j*) 20070701
-             IF(I1.LE.0.)  I1 = I1 + IMN
-             IF (I1 .GT. IMN) I1 = I1 - IMN
-             DO J1=JST(J),JEN(J)
-                HEIGHT = FLOAT(ZAVG(I1,J1))
-                IF(HEIGHT.LT.-990.) HEIGHT = 0.0
-                IF ( HEIGHT .gt. ORO(I,J) ) then
-                   if ( HEIGHT .gt. ZMAX(I,J) )ZMAX(I,J) = HEIGHT
-                   XNSUM(I,J) = XNSUM(I,J) + 1
-                  ENDIF
-             ENDDO
-            ENDDO
-           if ( I .lt. 5  .and. J .ge. JM-5 ) then
-      print *,' I,J,ORO(I,J),XNSUM(I,J),ZMAX(I,J):',
-     1         I,J,ORO(I,J),XNSUM(I,J),ZMAX(I,J)
-           endif
-        ENDDO
-      ENDDO
-!
-C....     make ELVMAX  ORO from MAKEMT sub
-C
-! ---  this will make work1 array take on oro's values on return
-      DO J=1,JM
-        DO I=1,numi(j)
-
-          ORO1(I,J) = ORO(I,J)
-          ELVMAX(I,J) = ZMAX(I,J) 
-        ENDDO
-      ENDDO
-C........
-C      The MAX elev peak (no averaging)
-C........
-!     DO J=1,JM
-!       DO I=1,numi(j)
-!           DO II1 = 1, IEN(I,J) - IST(I,J) + 1
-!              I1 = IST(I,J) + II1 - 1
-!              IF(I1.LE.0.)  I1 = I1 + IMN
-!              IF(I1.GT.IMN) I1 = I1 - IMN
-!              DO J1=JST(J),JEN(J)
-!              if ( ELVMAX(I,J) .lt. ZMAX(I1,J1)) 
-!    1              ELVMAX(I,J)   =  ZMAX(I1,J1)
-!              ENDDO
-!           ENDDO
-!        ENDDO
-!     ENDDO
-C
-C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
-C     IN A GRID BOX
-      DO J=1,JM
-        DO I=1,numi(j)
-          XNSUM1(I,J) = 0.0
-          XNSUM2(I,J) = 0.0
-          XNSUM3(I,J) = 0.0
-          XNSUM4(I,J) = 0.0
-        ENDDO
-      ENDDO
-! ---                 loop 
-      DO J=1,JM1
-        DO I=1,numi(j)
-           HC = 1116.2 - 0.878 * VAR(I,J)
-!          print *,' I,J,HC,VAR:',I,J,HC,VAR(I,J)
-            DO II1 = 1, IEN(I,J) - IST(I,J) + 1
-               I1 = IST(I,J) + II1 - 1
-!     IF (I1.LE.0.) print *,' I1 less than 0',I1,II1,IST(I,J),IEN(I,J) 
-!               if ( J .lt. 3 .or. J .gt. JM-2 ) then 
-!     IF(I1 .GT. IMN)print *,' I1 > IMN',J,I1,II1,IMN,IST(I,J),IEN(I,J) 
-!               endif
-             IF(I1.GT.IMN) I1 = I1 - IMN
-             DO J1=JST(J),JEN(J)
-               IF(FLOAT(ZAVG(I1,J1)) .GT. HC)
-     1            XNSUM1(I,J) = XNSUM1(I,J) + 1
-               XNSUM2(I,J) = XNSUM2(I,J) + 1
-             ENDDO
-           ENDDO
-C
-           INCI = NINT((IEN(I,j)-IST(I,j)) * 0.5)
-           ISTTT = MIN(MAX(IST(I,j)-INCI,1),IMN)
-           IEDDD = MIN(MAX(IEN(I,j)-INCI,1),IMN)
-C
-           INCJ = NINT((JEN(J)-JST(J)) * 0.5)
-           JSTTT = MIN(MAX(JST(J)-INCJ,1),JMN)
-           JEDDD = MIN(MAX(JEN(J)-INCJ,1),JMN)
-!               if ( J .lt. 3 .or. J .gt. JM-3 )  then
-!                 if(I .lt. 3 .or. I .gt. IM-3) then
-!        print *,' INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD:',
-!    1  I,J,INCI,ISTTT,IEDDD,INCJ,JSTTT,JEDDD  
-!                 endif
-!               endif
-C
-           DO I1=ISTTT,IEDDD
-             DO J1=JSTTT,JEDDD
-               IF(FLOAT(ZAVG(I1,J1)) .GT. HC)
-     1            XNSUM3(I,J) = XNSUM3(I,J) + 1
-               XNSUM4(I,J) = XNSUM4(I,J) + 1
-             ENDDO
-           ENDDO
-cx         print*,' i j hc var ',i,j,hc,var(i,j)
-cx         print*,'xnsum12 ',xnsum1(i,j),xnsum2(i,j)
-cx         print*,'xnsum34 ',xnsum3(i,j),xnsum4(i,j)
-        ENDDO
-      ENDDO
-C
-C---- CALCULATE THE 3D OROGRAPHIC ASYMMETRY FOR 4 WIND DIRECTIONS
-C---- AND THE 3D OROGRAPHIC SUBGRID OROGRAPHY FRACTION
-C     (KWD = 1  2  3  4)
-C     ( WD = W  S SW NW)
-C
-C
-      DO KWD = 1, 4
-        DO J=1,JM
-          DO I=1,numi(j)
-            OA4(I,J,KWD) = 0.0
-          ENDDO
-        ENDDO
-      ENDDO
-C
-      DO J=1,JM-2
-        DO I=1,numi(j)
-        II = I + 1
-        IF (II .GT. numi(j)) II = II - numi(j)
-          XNPU = XNSUM(I,J)    + XNSUM(I,J+1)
-          XNPD = XNSUM(II,J)   + XNSUM(II,J+1)
-          IF (XNPD .NE. XNPU) OA4(II,J+1,1) = 1. - XNPD / MAX(XNPU , 1.)
-          OL(II,J+1,1) = (XNSUM3(I,J+1)+XNSUM3(II,J+1))/
-     1                   (XNSUM4(I,J+1)+XNSUM4(II,J+1))
-!          if ( I .lt. 20  .and. J .ge. JM-19 ) then
-!        PRINT*,' MAKEOA: I J IST IEN ',I,j,IST(I,J),IEN(I,J)
-!        PRINT*,' HC VAR ',HC,VAR(i,j)
-!        PRINT*,' MAKEOA: XNSUM(I,J)=',XNSUM(I,J),XNPU, XNPD
-!        PRINT*,' MAKEOA: XNSUM3(I,J+1),XNSUM3(II,J+1)',
-!    1                    XNSUM3(I,J+1),XNSUM3(II,J+1)
-!        PRINT*,' MAKEOA: II, OA4(II,J+1,1), OL(II,J+1,1):',
-!    1                    II, OA4(II,J+1,1), OL(II,J+1,1)
-!          endif
-        ENDDO
-      ENDDO
-      DO J=1,JM-2
-        DO I=1,numi(j)
-        II = I + 1
-        IF (II .GT. numi(j)) II = II - numi(j)
-          XNPU = XNSUM(I,J+1)   + XNSUM(II,J+1)
-          XNPD = XNSUM(I,J)     + XNSUM(II,J)
-          IF (XNPD .NE. XNPU) OA4(II,J+1,2) = 1. - XNPD / MAX(XNPU , 1.)
-          OL(II,J+1,2) = (XNSUM3(II,J)+XNSUM3(II,J+1))/
-     1                   (XNSUM4(II,J)+XNSUM4(II,J+1))
-        ENDDO
-      ENDDO
-      DO J=1,JM-2
-        DO I=1,numi(j)
-        II = I + 1
-        IF (II .GT. numi(j)) II = II - numi(j)
-          XNPU = XNSUM(I,J+1)  + ( XNSUM(I,J) + XNSUM(II,J+1) )*0.5
-          XNPD = XNSUM(II,J)   + ( XNSUM(I,J) + XNSUM(II,J+1) )*0.5
-          IF (XNPD .NE. XNPU) OA4(II,J+1,3) = 1. - XNPD / MAX(XNPU , 1.)
-          OL(II,J+1,3) = (XNSUM1(II,J)+XNSUM1(I,J+1))/
-     1                   (XNSUM2(II,J)+XNSUM2(I,J+1))
-        ENDDO
-      ENDDO
-      DO J=1,JM-2
-        DO I=1,numi(j)
-        II = I + 1
-        IF (II .GT. numi(j)) II = II - numi(j)
-          XNPU = XNSUM(I,J)    + ( XNSUM(II,J) + XNSUM(I,J+1) )*0.5
-          XNPD = XNSUM(II,J+1) + ( XNSUM(II,J) + XNSUM(I,J+1) )*0.5
-          IF (XNPD .NE. XNPU) OA4(II,J+1,4) = 1. - XNPD / MAX(XNPU , 1.)
-          OL(II,J+1,4) = (XNSUM1(I,J)+XNSUM1(II,J+1))/
-     1                   (XNSUM2(I,J)+XNSUM2(II,J+1))
-        ENDDO
-      ENDDO
-C
-      DO KWD = 1, 4
-        DO I=1,numi(j)
-          OL(I,1,KWD)  = OL(I,2,KWD)
-          OL(I,JM,KWD) = OL(I,JM-1,KWD)
-        ENDDO
-      ENDDO
-C
-      DO KWD=1,4
-        DO J=1,JM
-          DO I=1,numi(j)
-            T = OA4(I,J,KWD)
-            OA4(I,J,KWD) = SIGN( MIN( ABS(T), 1. ), T )
-          ENDDO
-        ENDDO
-      ENDDO
-C
-      NS0 = 0
-      NS1 = 0
-      NS2 = 0
-      NS3 = 0
-      NS4 = 0
-      NS5 = 0
-      NS6 = 0
-      DO KWD=1,4
-      DO J=1,JM
-      DO I=1,numi(j)
-         T = ABS( OA4(I,J,KWD) )
-         IF(T .EQ. 0.) THEN
-            IOA4(I,J,KWD) = 0
-            NS0 = NS0 + 1
-         ELSE IF(T .GT. 0. .AND. T .LE. 1.) THEN
-            IOA4(I,J,KWD) = 1
-            NS1 = NS1 + 1
-         ELSE IF(T .GT. 1. .AND. T .LE. 10.) THEN
-            IOA4(I,J,KWD) = 2
-            NS2 = NS2 + 1
-         ELSE IF(T .GT. 10. .AND. T .LE. 100.) THEN
-            IOA4(I,J,KWD) = 3
-            NS3 = NS3 + 1
-         ELSE IF(T .GT. 100. .AND. T .LE. 1000.) THEN
-            IOA4(I,J,KWD) = 4
-            NS4 = NS4 + 1
-         ELSE IF(T .GT. 1000. .AND. T .LE. 10000.) THEN
-            IOA4(I,J,KWD) = 5
-            NS5 = NS5 + 1
-         ELSE IF(T .GT. 10000.) THEN
-            IOA4(I,J,KWD) = 6
-            NS6 = NS6 + 1
-         ENDIF
-      ENDDO
-      ENDDO
-      ENDDO
-C
-      WRITE(6,*) "! MAKEOA EXIT"
-C
-      RETURN
-      END SUBROUTINE MAKEOA
-
 !> Convert the 'x' direction distance of a cubed-sphere grid
 !! point to the corresponding distance in longitude.
 !!
@@ -3238,22 +2109,20 @@ SUBROUTINE MAKEOA2(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
       real    XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
       real    XNSUM3(IM,JM),XNSUM4(IM,JM)
       real    VAR(IM,JM),OL(IM,JM,4)
-      LOGICAL FLAG
       integer i,j,ilist(IMN),numx,i1,j1,ii1
-      integer KWD,II,npts
+      integer KWD
       real    LONO(4),LATO(4),LONI,LATI
       real    DELXN,HC,HEIGHT,XNPU,XNPD,T
       integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
       logical inside_a_polygon
       real    lon,lat,dlon,dlat,dlat_old
       real    lon1,lat1,lon2,lat2
-      real    xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
+      real    xnsum11,xnsum12,xnsum21,xnsum22
       real    HC_11, HC_12, HC_21, HC_22
       real    xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
       real    xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
       real    get_lon_angle, get_lat_angle, get_xnsum
-      integer ist, ien, jst, jen
-      real    xland,xwatr,xl1,xs1,oroavg,slm
+      integer jst, jen
 C   
 C---- GLOBAL XLAT AND XLON ( DEGREE )
 C
@@ -3748,7 +2617,6 @@ end subroutine interpolate_mismatch
 !! is computed from the high-resolution orography data.
 !!
 !! @param[in] zavg High-resolution orography data.
-!! @param[in] zslm High-resolution land-mask data.  Not used.
 !! @param[in] var Standard deviation of orography on the model grid.
 !! @param[out] glat Latitude of each row of input terrain dataset.
 !! @param[out] oa4 Orographic asymmetry on the model grid. Four
@@ -3775,8 +2643,6 @@ end subroutine interpolate_mismatch
 !! @param[in] lat_c Corner point latitudes of the model grid points.
 !! @param[in] lon_t Center point longitudes of the model grid points.
 !! @param[in] lat_t Center point latitudes of the model grid points.
-!! @param[in] is_south_pole Not used.
-!! @param[in] is_north_pole Not used.
 !! @param[in] imi 'i' dimension of input gfs orography data.
 !! @param[in] jmi 'j' dimension of input gfs orography data.
 !! @param[in] oa_in Asymmetry on the input gfs orography data.
@@ -3785,10 +2651,10 @@ end subroutine interpolate_mismatch
 !! @param[in] lon_in Longitude on the input gfs orography data.
 !! @param[in] lat_in Latitude on the input gfs orography data.
 !! @author Jordan Alpert NOAA/EMC
-      SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
+      SUBROUTINE MAKEOA3(ZAVG,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
      1           ORO,SLM,oro1,XNSUM,XNSUM1,XNSUM2,XNSUM3,XNSUM4,
      2           IM,JM,IMN,JMN,lon_c,lat_c,lon_t,lat_t,
-     3           is_south_pole,is_north_pole,IMI,JMI,OA_IN,OL_IN,
+     3           IMI,JMI,OA_IN,OL_IN,
      4           slm_in,lon_in,lat_in)
 
 ! Required when using iplib v4.0 or higher.
@@ -3802,7 +2668,7 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
       real, PARAMETER :: R2D=180./3.14159265358979
       integer IM,JM,IMN,JMN,IMI,JMI
       real    GLAT(JMN)
-      INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
+      INTEGER ZAVG(IMN,JMN)
       real    SLM(IM,JM)
       real    ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
       real    OA4(IM,JM,4)
@@ -3812,26 +2678,16 @@ SUBROUTINE MAKEOA3(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
       real    lon_in(IMI,JMI), lat_in(IMI,JMI)
       real    lon_c(IM+1,JM+1), lat_c(IM+1,JM+1)
       real    lon_t(IM,JM), lat_t(IM,JM)
-      logical is_south_pole(IM,JM), is_north_pole(IM,JM)
       real    XNSUM(IM,JM),XNSUM1(IM,JM),XNSUM2(IM,JM)
       real    XNSUM3(IM,JM),XNSUM4(IM,JM)
       real    VAR(IM,JM),OL(IM,JM,4)
-      LOGICAL FLAG
       integer i,j,ilist(IMN),numx,i1,j1,ii1
-      integer KWD,II,npts
+      integer KWD
       real    LONO(4),LATO(4),LONI,LATI
-      real    DELXN,HC,HEIGHT,XNPU,XNPD,T
+      real    DELXN,HC,HEIGHT,T
       integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
       logical inside_a_polygon
-      real    lon,lat,dlon,dlat,dlat_old
-      real    lon1,lat1,lon2,lat2
-      real    xnsum11,xnsum12,xnsum21,xnsum22,xnsumx
-      real    HC_11, HC_12, HC_21, HC_22
-      real    xnsum1_11,xnsum1_12,xnsum1_21,xnsum1_22
-      real    xnsum2_11,xnsum2_12,xnsum2_21,xnsum2_22
-      real    get_lon_angle, get_lat_angle, get_xnsum
-      integer ist, ien, jst, jen
-      real    xland,xwatr,xl1,xs1,oroavg
+      integer jst, jen
       integer int_opt, ipopt(20), kgds_input(200), kgds_output(200)
       integer count_land_output
       integer ij, ijmdl_output, iret, num_mismatch_land, num
@@ -4182,15 +3038,13 @@ END SUBROUTINE MAKEOA3
 !!
 !! @param [in] im 'i' dimension of the 2-d array.
 !! @param [in] jm 'j' dimension of the 2-d array.
-!! @param [in] numi Not used.
 !! @param [inout] f The two-dimensional array to
 !! be processed.
 !! @param [out] wrk Two-dimensional work array.
 !! @author Jordan Alpert NOAA/EMC
-      SUBROUTINE REVERS(IM, JM, numi, F, WRK)
+      SUBROUTINE REVERS(IM, JM, F, WRK)
 !
       REAL F(IM,JM), WRK(IM,JM)
-      integer numi(jm)
       imb2 = im / 2
       do i=1,im*jm
          WRK(i,1) = F(i,1)
@@ -4339,91 +3193,11 @@ SUBROUTINE mnmxja(IM,JM,A,imax,jmax,title)
       RETURN
       END
 
-!> Perform multiple fast fourier transforms.
-!!
-!! This subprogram performs multiple fast fourier transforms
-!! between complex amplitudes in fourier space and real values
-!! in cyclic physical space.
-!!
-!! Subprograms called (NCEPLIB SP Library):
-!!  - scrft Complex to real fourier transform
-!!  - dcrft Complex to real fourier transform
-!!  - srcft Real to complex fourier transform
-!!  - drcft Real to complex fourier transform
-!!
-!! Program history log:
-!! 1998-12-18  Mark Iredell
-!!
-!! @param[in] imax Integer number of values in the cyclic physical
-!! space. See limitations on imax in remarks below.
-!! @param[in] incw Integer first dimension of the complex amplitude array.
-!! (incw >= imax/2+1).
-!! @param[in] incg Integer first dimension of the real value array.
-!! (incg >= imax).
-!! @param[in] kmax Integer number of transforms to perform.
-!! @param[in] w Complex amplitudes on input if idir>0, and on output
-!! if idir<0.
-!! @param[in] g Real values on input if idir<0, and on output if idir>0.
-!! @param[in] idir Integer direction flag. idir>0 to transform from
-!! fourier to physical space. idir<0 to transform from physical to
-!! fourier space.
-!!
-!! @note The restrictions on imax are that it must be a multiple
-!! of 1 to 25 factors of two, up to 2 factors of three,
-!! and up to 1 factor of five, seven and eleven.
-!!
-!! @author Mark Iredell ORG: W/NMC23 @date 96-02-20
-      SUBROUTINE SPFFT1(IMAX,INCW,INCG,KMAX,W,G,IDIR)
-        IMPLICIT NONE
-        INTEGER,INTENT(IN):: IMAX,INCW,INCG,KMAX,IDIR
-        COMPLEX,INTENT(INOUT):: W(INCW,KMAX)
-        REAL,INTENT(INOUT):: G(INCG,KMAX)
-        REAL:: AUX1(25000+INT(0.82*IMAX))
-        REAL:: AUX2(20000+INT(0.57*IMAX))
-        INTEGER:: NAUX1,NAUX2
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-        NAUX1=25000+INT(0.82*IMAX)
-        NAUX2=20000+INT(0.57*IMAX)
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-C  FOURIER TO PHYSICAL TRANSFORM.
-        SELECT CASE(IDIR)
-        CASE(1:)
-          SELECT CASE(DIGITS(1.))
-          CASE(DIGITS(1._4))
-            CALL SCRFT(1,W,INCW,G,INCG,IMAX,KMAX,-1,1.,
-     &                 AUX1,NAUX1,AUX2,NAUX2,0.,0)
-            CALL SCRFT(0,W,INCW,G,INCG,IMAX,KMAX,-1,1.,
-     &                 AUX1,NAUX1,AUX2,NAUX2,0.,0)
-          CASE(DIGITS(1._8))
-            CALL DCRFT(1,W,INCW,G,INCG,IMAX,KMAX,-1,1.,
-     &                 AUX1,NAUX1,AUX2,NAUX2,0.,0)
-            CALL DCRFT(0,W,INCW,G,INCG,IMAX,KMAX,-1,1.,
-     &                 AUX1,NAUX1,AUX2,NAUX2,0.,0)
-          END SELECT
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-C  PHYSICAL TO FOURIER TRANSFORM.
-        CASE(:-1)
-          SELECT CASE(DIGITS(1.))
-          CASE(DIGITS(1._4))
-            CALL SRCFT(1,G,INCG,W,INCW,IMAX,KMAX,+1,1./IMAX,
-     &               AUX1,NAUX1,AUX2,NAUX2,0.,0)
-            CALL SRCFT(0,G,INCG,W,INCW,IMAX,KMAX,+1,1./IMAX,
-     &               AUX1,NAUX1,AUX2,NAUX2,0.,0)
-          CASE(DIGITS(1._8))
-            CALL DRCFT(1,G,INCG,W,INCW,IMAX,KMAX,+1,1./IMAX,
-     &               AUX1,NAUX1,AUX2,NAUX2,0.,0)
-            CALL DRCFT(0,G,INCG,W,INCW,IMAX,KMAX,+1,1./IMAX,
-     &               AUX1,NAUX1,AUX2,NAUX2,0.,0)
-          END SELECT
-        END SELECT
-      END SUBROUTINE
-
 !> Read input global 30-arc second orography data.
 !!
 !! @param[out] glob The orography data.
-!! @param[in] itopo Not used.
 !! @author Jordan Alpert NOAA/EMC
-      subroutine read_g(glob,ITOPO)
+      subroutine read_g(glob)
       implicit none
 cc
       integer*2    glob(360*120,180*120)
@@ -4434,10 +3208,7 @@ subroutine read_g(glob,ITOPO)
       parameter   (ix=40*120,jx=50*120)
       parameter   (ia=60*120,ja=30*120)
 cc
-      integer*2    idat(ix,jx)
-      integer      itopo
-cc
-      integer      i,j,inttyp
+      integer      inttyp
 cc
       real(kind=8) dloin,dlain,rlon,rlat
 cc
@@ -4480,7 +3251,7 @@ subroutine maxmin(ia,len,tile)
 ccmr
       integer*2 ia(len)
       character*7 tile
-      integer iaamax, iaamin, len, j, m, ja, kount
+      integer iaamax, iaamin, len, m, ja, kount
       integer(8) sum2,std,mean,isum
       integer i_count_notset,kount_9
 ! --- missing is -9999
@@ -4756,7 +3527,6 @@ function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
         implicit none
 
         real get_xnsum
-        logical verbose
         real, intent(in) :: lon1,lat1,lon2,lat2,delxn
         integer, intent(in) ::  IMN,JMN
         real, intent(in) :: glat(JMN)
@@ -4783,7 +3553,6 @@ function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
         ien = lon2/delxn
         if(ist .le.0) ist = ist + IMN
         if(ien < ist) ien = ien + IMN
-!        if(verbose) print*, "ist,ien=",ist,ien,jst,jen
 
         !--- compute average oro
           oro = 0.0
@@ -4826,7 +3595,6 @@ function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
                IF ( HEIGHT .gt. ORO ) get_xnsum = get_xnsum + 1
             enddo       
          enddo
-!         if(verbose) print*, "get_xnsum=", get_xnsum, oro
          
       end function get_xnsum  
       
@@ -4862,14 +3630,13 @@ subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
         implicit none
 
         real, intent(out) :: xnsum1,xnsum2,HC
-        logical verbose
-        real lon1,lat1,lon2,lat2,oro,delxn
+        real lon1,lat1,lon2,lat2,delxn
         integer IMN,JMN
         real    glat(JMN)
         integer zavg(IMN,JMN)
         integer i, j, ist, ien, jst, jen, i1
         real    HEIGHT, var
-        real    XW1,XW2,slm,xnsum
+        real    XW1,XW2,xnsum
         !---figure out ist,ien,jst,jen
         do j = 1, JMN
            if( GLAT(J) .GT. lat1 ) then
@@ -4889,7 +3656,6 @@ subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
         ien = lon2/delxn
         if(ist .le.0) ist = ist + IMN
         if(ien < ist) ien = ien + IMN
-!        if(verbose) print*, "ist,ien=",ist,ien,jst,jen
 
         !--- compute average oro
           xnsum = 0
@@ -4957,13 +3723,12 @@ subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
         implicit none
 
         real, intent(out) :: xnsum1,xnsum2
-        real lon1,lat1,lon2,lat2,oro,delxn
+        real lon1,lat1,lon2,lat2,delxn
         integer IMN,JMN
         real    glat(JMN)
         integer zavg(IMN,JMN)
         integer i, j, ist, ien, jst, jen, i1
         real    HEIGHT, HC
-        real    XW1,XW2,slm,xnsum
         !---figure out ist,ien,jst,jen
         ! if lat1 or lat 2 is 90 degree. set jst = JMN
         jst = JMN
@@ -4986,7 +3751,6 @@ subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
         ien = lon2/delxn
         if(ist .le.0) ist = ist + IMN
         if(ien < ist) ien = ien + IMN
-!        if(verbose) print*, "ist,ien=",ist,ien,jst,jen
 
          xnsum1 = 0
          xnsum2 = 0
@@ -5035,8 +3799,7 @@ subroutine nanc(a,l,c)
       data inaq3/x'FFC00000'/
       data inaq4/x'FFFFFFFF'/
 c
-      real(kind=8)a(l),rtc,t1,t2
-      character*24 cn
+      real(kind=8)a(l)
       character*(*) c
 c     t1=rtc()
 cgwv        print *, ' nanc call ',c
diff --git a/sorc/orog_mask_tools.fd/orog.fd/netcdf_io.F90 b/sorc/orog_mask_tools.fd/orog.fd/netcdf_io.F90
index 09d994b0b..7ff8ce725 100644
--- a/sorc/orog_mask_tools.fd/orog.fd/netcdf_io.F90
+++ b/sorc/orog_mask_tools.fd/orog.fd/netcdf_io.F90
@@ -25,7 +25,7 @@ subroutine write_netcdf(im, jm, slm, land_frac, oro, orf, hprime, ntiles, tile,
     real, intent(in), dimension(im,jm)  :: slm, oro, orf, geolon, geolat, land_frac
     real, intent(in), dimension(im,jm,14):: hprime
     character(len=128) :: outfile
-    integer            :: error, ncid, i
+    integer            :: error, ncid
     integer            :: header_buffer_val = 16384      
     integer            :: fsize=65536, inital = 0  
     integer            :: dim1, dim2
@@ -245,7 +245,7 @@ subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geola
     integer, intent(in):: im, jm, ntiles, tile
     real, intent(in), dimension(im,jm)  :: slm, geolon, geolat, land_frac
     character(len=128) :: outfile
-    integer            :: error, ncid, i
+    integer            :: error, ncid
     integer            :: header_buffer_val = 16384      
     integer            :: fsize=65536, inital = 0  
     integer            :: dim1, dim2
diff --git a/ush/fv3gfs_make_orog.sh b/ush/fv3gfs_make_orog.sh
index a684cf6e4..4693e47d0 100755
--- a/ush/fv3gfs_make_orog.sh
+++ b/ush/fv3gfs_make_orog.sh
@@ -57,15 +57,9 @@ fi
 if [ ! -s $workdir ]; then mkdir -p $workdir ;fi
 if [ ! -s $outdir ]; then mkdir -p $outdir ;fi
 
-#jcap is for Gaussian grid
-#jcap=`expr $latb - 2 `
-jcap=0
-NF1=0
-NF2=0
 mtnres=1
 efac=0
 blat=0
-NR=0
 
 if [ $is_latlon -eq 1 ]; then
   OUTGRID="none"
@@ -96,7 +90,7 @@ if [ $is_latlon -eq 0 ]; then
 fi
 cp $executable .
 
-echo  $mtnres $lonb $latb $jcap $NR $NF1 $NF2 $efac $blat > INPS
+echo  $mtnres $lonb $latb $efac $blat > INPS
 echo $OUTGRID >> INPS
 echo $orogfile >> INPS
 if [ -z ${ocn+x} ]; then