From 57bd8321764f39a9e1926d104a9207f800bcc80d Mon Sep 17 00:00:00 2001
From: GeorgeGayno-NOAA <52789452+GeorgeGayno-NOAA@users.noreply.github.com>
Date: Tue, 6 Feb 2024 10:04:05 -0500
Subject: [PATCH] Revert "orog.fd - Reduce dependency on SPLIB and some cleanup
 (#889)" (#893)

This reverts commit 6dee3922c34a1a9a1e4ebc5fa5f266ffa07afd7c.
---
 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, 1356 insertions(+), 112 deletions(-)

diff --git a/reg_tests/grid_gen/driver.hercules.sh b/reg_tests/grid_gen/driver.hercules.sh
index 5a323249c..3c46a2d45 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:-/work2/noaa/stmp/$LOGNAME}"
+export WORK_DIR="${WORK_DIR:-/work/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 5fba04abe..369796d61 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:-/work2/noaa/stmp/$LOGNAME}"
+export DATA="${WORK_DIR:-/work/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 bc0e5d5b1..04ab86742 100644
--- a/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
+++ b/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt
@@ -25,6 +25,7 @@ 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 add54f00c..96e3b38d8 100644
--- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
+++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F
@@ -1,9 +1,25 @@
 C> @file
-C> Terrain maker for the ufs weather model.
+C> Terrain maker for global spectral model.
 C> @author Mark Iredell @date 92-04-16
   
-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> 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>
 C> PROGRAM HISTORY LOG:
 C> -  92-04-16  IREDELL
@@ -21,76 +37,90 @@
 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      - 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>  -   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>  -   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>  -   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>  -   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>
 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, imn, jmn
+      integer fsize, ncid, error, id_dim, nx, ny
       character(len=256) :: OUTGRID = "none"
       character(len=256) :: INPUTOROG = "none"
       character(len=256) :: merge_file = "none"
       logical :: mask_only = .false.
-      integer :: MTNRES,IM,JM,EFAC,BLAT
+      integer :: MTNRES,IM,JM,NM,NR,NF0,NF1,EFAC,BLAT,NW
       fsize=65536
-      READ(5,*) MTNRES,IM,JM,EFAC,BLAT
+      READ(5,*) MTNRES,IM,JM,NM,NR,NF0,NF1,EFAC,BLAT
       READ(5,*) OUTGRID
       READ(5,*) INPUTOROG
       READ(5,*) mask_only
       READ(5,*) 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=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)
 ! --- 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.
-
-      print*, "MTNRES,EFAC,BLAT: ",MTNRES,EFAC,BLAT
-
+      MTNRES=1
+      print*, MTNRES,IM,JM,NM,NR,NF0,NF1,EFAC,BLAT
+      NW=(NM+1)*((NR+1)*NM+2)
       IMN = 360*120/MTNRES
       JMN = 180*120/MTNRES
+      print *, ' Starting terr12 mtnlm7_slm30.f  IMN,JMN:',IMN,JMN
 
-! --- Read the model grid resolution from 'grid' file OUTGRID.
-
-      READ_GRID_FILE : if( trim(OUTGRID) .NE. "none" ) then
-
+! --- read the grid resolution if the OUTGRID exists.
+      if( trim(OUTGRID) .NE. "none" ) then
          inquire(file=trim(OUTGRID), exist=fexist)
          if(.not. fexist) then
-            print*, "FATAL ERROR: file "
-     &             //trim(OUTGRID)//" does not exist"
+            print*, "file "//trim(OUTGRID)//" does not exist"
             CALL ERREXIT(4)
          endif
          do ncid = 103, 512
@@ -98,7 +128,7 @@
            if( .NOT.opened )exit
          end do
 
-         print*, "READ outgrid=", trim(outgrid)
+         print*, "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)
@@ -127,15 +157,10 @@
          endif
          error=nf_close(ncid)
          call netcdf_err(error, 'close file '//trim(OUTGRID) )
-
-      else
-
-        print*, "FATAL ERROR: SPECIFY OUTGRID"
-        CALL ERREXIT(8)
          
-      endif READ_GRID_FILE
+      endif         
          
-      CALL TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
+      CALL TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
      &            OUTGRID,INPUTOROG,MASK_ONLY,MERGE_FILE)
       STOP
       END
@@ -146,6 +171,11 @@
 !! @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.
@@ -157,12 +187,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,EFAC,BLAT,
+      SUBROUTINE TERSUB(IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW,EFAC,BLAT,
      &     OUTGRID,INPUTOROG,MASK_ONLY,MERGE_FILE)
       implicit none
       include 'netcdf.inc'
 C
-      integer, intent(in)          :: IMN,JMN,IM,JM
+      integer                      :: IMN,JMN,IM,JM,NM,NR,NF0,NF1,NW
       character(len=*), intent(in) :: OUTGRID
       character(len=*), intent(in) :: INPUTOROG
       character(len=*), intent(in) :: MERGE_FILE
@@ -170,19 +200,20 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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 :: IMT,ios,iosg,latg2,istat,itest,jtest
+      integer :: M,N,IMT,IRET,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 :: numi(:)
+      integer, allocatable :: JST(:),JEN(:),KPDS(:),KGDS(:),numi(:)
       integer, allocatable :: lonsperlat(:)
 
       integer, allocatable :: IST(:,:),IEN(:,:),ZSLMX(:,:)
@@ -192,11 +223,12 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
 
       integer, allocatable :: IWORK(:,:,:)
 
-      real :: maxlat, minlat,timef,tbeg,tend,tbeg1
-      real :: DELXN,RS,RN,slma,oroa,vara,var4a,xn,XS
+      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, allocatable :: WGTCLT(:),XLAT(:)
-      real, allocatable :: XLON(:),oaa(:),ola(:),GLAT(:)
+      real, allocatable :: COSCLT(:),WGTCLT(:),RCLT(:),XLAT(:),DIFFX(:)
+      real, allocatable :: XLON(:),ORS(:),oaa(:),ola(:),GLAT(:)
 
       real, allocatable :: GEOLON(:,:),GEOLON_C(:,:),DX(:,:)
       real, allocatable :: GEOLAT(:,:),GEOLAT_C(:,:),DY(:,:)
@@ -213,23 +245,28 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
       real, allocatable :: OA(:,:,:),OL(:,:,:),HPRIME(:,:,:)
       real, allocatable :: oa_in(:,:,:), ol_in(:,:,:)
 
-      logical :: fexist,opened
-      logical :: REVLAT
+
+      complex :: ffj(im/2+1)
+
+      logical :: grid_from_file,output_binary,fexist,opened
+      logical :: SPECTR, REVLAT, FILTER
 
       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 (numi(jm))
+      allocate (JST(JM),JEN(JM),KPDS(200),KGDS(200),numi(jm))
       allocate (lonsperlat(jm/2))
       allocate (IST(IM,jm),IEN(IM,jm),ZSLMX(2700,1350))
       allocate (glob(IMN,JMN))
 
 ! reals
-      allocate (WGTCLT(JM),XLAT(JM))
-      allocate (XLON(IM),oaa(4),ola(4),GLAT(JMN))
+      allocate (COSCLT(JM),WGTCLT(JM),RCLT(JM),XLAT(JM),DIFFX(JM/2))
+      allocate (XLON(IM),ORS(NW),oaa(4),ola(4),GLAT(JMN))
 
       allocate (ZAVG(IMN,JMN))
       allocate (ZSLM(IMN,JMN))
@@ -238,6 +275,9 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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
@@ -320,7 +360,7 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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)
+          if ( ITOPO .ne. 0 ) call read_g(glob,ITOPO)
 ! --- 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
@@ -343,8 +383,8 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
 !
 
 !  --- IMN,JMN
-      print*, ' IM, JM, EFAC, BLAT'
-      print*, IM,JM,EFAC,BLAT
+      print*, ' IM, JM, NM, NR, NF0, NF1, EFAC, BLAT'
+      print*, IM,JM,NM,NR,NF0,NF1,EFAC,BLAT
        print *,'  imn,jmn,glob(imn,jmn)=',imn,jmn,glob(imn,jmn)
        print *,' UBOUND ZAVG=',UBOUND(ZAVG)
        print *,' UBOUND glob=',UBOUND(glob)
@@ -467,14 +507,48 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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 *,' REVLAT=',REVLAT,' ** with GICE-07 **'
+
+      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
 
       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
@@ -592,15 +666,15 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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
-
-! Read 'grid' file
-
+      if( trim(OUTGRID) .NE. "none" ) then
+         grid_from_file = .true.
          inquire(file=trim(OUTGRID), exist=fexist)
          if(.not. fexist) then
             print*, "file "//trim(OUTGRID)//" does not exist"
@@ -619,7 +693,25 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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))
@@ -726,6 +818,7 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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 '
@@ -741,10 +834,31 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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
@@ -770,6 +884,10 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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')
@@ -796,12 +914,16 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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')
@@ -821,7 +943,7 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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()
@@ -935,15 +1057,21 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
      &                   //trim(INPUTOROG) )
            
            print*, "calling MAKEOA3 to compute OA, OL"
-           CALL MAKEOA3(ZAVG,VAR,GLAT,OA,OL,IWORK,ELVMAX,ORO,SLM,
+           CALL MAKEOA3(ZAVG,zslm,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,nx_in,ny_in,
+     3            geolon,geolat,is_south_pole,is_north_pole,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)
@@ -1238,23 +1366,67 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,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 (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)
+      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
       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
@@ -1287,14 +1459,152 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
       ENDDO
       tend=timef()
       write(6,*)' Timer 5 time= ',tend-tbeg
-
-      do j = 1, jm
-         xlat(j) = geolat(1,j)
-      enddo
-      do i = 1, im
-         xlon(i) = geolon(i,1)
+      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)
       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
       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)
@@ -1305,8 +1615,8 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
       print *,' ===== Deallocate Arrays and ENDING MTN VAR OROG program'
 
 !     Deallocate 1d vars
-      deallocate(numi,lonsperlat)
-      deallocate(WGTCLT,XLAT,XLON,oaa,ola,GLAT)
+      deallocate(JST,JEN,KPDS,KGDS,numi,lonsperlat)
+      deallocate(COSCLT,WGTCLT,RCLT,XLAT,DIFFX,XLON,ORS,oaa,ola,GLAT)
 
 !     Deallocate 2d vars
       deallocate (OCLSM)
@@ -1318,7 +1628,183 @@ SUBROUTINE TERSUB(IMN,JMN,IM,JM,EFAC,BLAT,
       tend=timef()
       write(6,*)' Total runtime time= ',tend-tbeg1
       RETURN
-      END SUBROUTINE TERSUB
+      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
 
 !> Determine the location of a cubed-sphere point within
 !! the high-resolution orography data.  The location is
@@ -1764,7 +2250,290 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
       deallocate(hgt_1d)
       deallocate(hgt_1d_all)
       RETURN
-      END SUBROUTINE MAKEMT2
+      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
 
 !> Make the principle coordinates - slope of orography, 
 !! anisotropy, angle of mountain range with respect to east.
@@ -1809,7 +2578,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 DEBUG
+      LOGICAL FLAG, DEBUG
 C===  DATA DEBUG/.TRUE./
       DATA DEBUG/.FALSE./
 C
@@ -2013,6 +2782,366 @@ 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.
 !!
@@ -2109,20 +3238,22 @@ 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
+      integer KWD,II,npts
       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
+      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 jst, jen
+      integer ist, ien, jst, jen
+      real    xland,xwatr,xl1,xs1,oroavg,slm
 C   
 C---- GLOBAL XLAT AND XLON ( DEGREE )
 C
@@ -2617,6 +3748,7 @@ 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
@@ -2643,6 +3775,8 @@ 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.
@@ -2651,10 +3785,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,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
+      SUBROUTINE MAKEOA3(ZAVG,zslm,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           IMI,JMI,OA_IN,OL_IN,
+     3           is_south_pole,is_north_pole,IMI,JMI,OA_IN,OL_IN,
      4           slm_in,lon_in,lat_in)
 
 ! Required when using iplib v4.0 or higher.
@@ -2668,7 +3802,7 @@ SUBROUTINE MAKEOA3(ZAVG,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)
+      INTEGER ZAVG(IMN,JMN),ZSLM(IMN,JMN)
       real    SLM(IM,JM)
       real    ORO(IM,JM),ORO1(IM,JM),ELVMAX(IM,JM),ZMAX(IM,JM)
       real    OA4(IM,JM,4)
@@ -2678,16 +3812,26 @@ SUBROUTINE MAKEOA3(ZAVG,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
+      integer KWD,II,npts
       real    LONO(4),LATO(4),LONI,LATI
-      real    DELXN,HC,HEIGHT,T
+      real    DELXN,HC,HEIGHT,XNPU,XNPD,T
       integer NS0,NS1,NS2,NS3,NS4,NS5,NS6
       logical inside_a_polygon
-      integer jst, jen
+      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 int_opt, ipopt(20), kgds_input(200), kgds_output(200)
       integer count_land_output
       integer ij, ijmdl_output, iret, num_mismatch_land, num
@@ -3038,13 +4182,15 @@ 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, F, WRK)
+      SUBROUTINE REVERS(IM, JM, numi, 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)
@@ -3193,11 +4339,91 @@ 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)
+      subroutine read_g(glob,ITOPO)
       implicit none
 cc
       integer*2    glob(360*120,180*120)
@@ -3208,7 +4434,10 @@ subroutine read_g(glob)
       parameter   (ix=40*120,jx=50*120)
       parameter   (ia=60*120,ja=30*120)
 cc
-      integer      inttyp
+      integer*2    idat(ix,jx)
+      integer      itopo
+cc
+      integer      i,j,inttyp
 cc
       real(kind=8) dloin,dlain,rlon,rlat
 cc
@@ -3251,7 +4480,7 @@ subroutine maxmin(ia,len,tile)
 ccmr
       integer*2 ia(len)
       character*7 tile
-      integer iaamax, iaamin, len, m, ja, kount
+      integer iaamax, iaamin, len, j, m, ja, kount
       integer(8) sum2,std,mean,isum
       integer i_count_notset,kount_9
 ! --- missing is -9999
@@ -3527,6 +4756,7 @@ 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)
@@ -3553,6 +4783,7 @@ 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
@@ -3595,6 +4826,7 @@ 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  
       
@@ -3630,13 +4862,14 @@ subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
         implicit none
 
         real, intent(out) :: xnsum1,xnsum2,HC
-        real lon1,lat1,lon2,lat2,delxn
+        logical verbose
+        real lon1,lat1,lon2,lat2,oro,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,xnsum
+        real    XW1,XW2,slm,xnsum
         !---figure out ist,ien,jst,jen
         do j = 1, JMN
            if( GLAT(J) .GT. lat1 ) then
@@ -3656,6 +4889,7 @@ 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
@@ -3723,12 +4957,13 @@ subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
         implicit none
 
         real, intent(out) :: xnsum1,xnsum2
-        real lon1,lat1,lon2,lat2,delxn
+        real lon1,lat1,lon2,lat2,oro,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
@@ -3751,6 +4986,7 @@ 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
@@ -3799,7 +5035,8 @@ subroutine nanc(a,l,c)
       data inaq3/x'FFC00000'/
       data inaq4/x'FFFFFFFF'/
 c
-      real(kind=8)a(l)
+      real(kind=8)a(l),rtc,t1,t2
+      character*24 cn
       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 7ff8ce725..09d994b0b 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
+    integer            :: error, ncid, i
     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
+    integer            :: error, ncid, i
     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 4693e47d0..a684cf6e4 100755
--- a/ush/fv3gfs_make_orog.sh
+++ b/ush/fv3gfs_make_orog.sh
@@ -57,9 +57,15 @@ 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"
@@ -90,7 +96,7 @@ if [ $is_latlon -eq 0 ]; then
 fi
 cp $executable .
 
-echo  $mtnres $lonb $latb $efac $blat > INPS
+echo  $mtnres $lonb $latb $jcap $NR $NF1 $NF2 $efac $blat > INPS
 echo $OUTGRID >> INPS
 echo $orogfile >> INPS
 if [ -z ${ocn+x} ]; then