mtnlm7_oclsm.F

Go to the documentation of this file.

C> @file
C> Terrain maker for global spectral model.
C> @author Mark Iredell @date 92-04-16
  
C> This program creates 7 terrain-related files computed from the
C> GMTED2010 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
C> -  98-02-02  IREDELL  FILTER
C> -  98-05-31  HONG Modified for subgrid orography used in Kim's scheme
C> -  98-12-31  HONG Modified for high-resolution GTOPO orography
C> -  99-05-31  HONG Modified for getting OL4 (mountain fraction)
C>  - 00-02-10  Moorthi's modifications
C> -  00-04-11  HONG Modified for reduced grids
C> -  00-04-12  Iredell Modified for reduced grids
C> -  02-01-07  (*j*) modified for principal axes of orography
C>             There are now 14 files, 4 additional for lm mb
C>  - 04-04-04  (*j*) re-Test on IST/ilen calc for sea-land mask(*j*)
C>  - 04-09-04   minus sign here in MAKEOA IST and IEN as in MAKEMT!
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>  -   NCID       - GMTED2010 USGS orography (NetCDF)
C>  -   NCID       - 30" UMD land cover mask. (NetCDF)
C>  -   NCID       - GICE Grumbine 30" RAMP Antarctica orog IMNx3601. (NetCDF)
C>  -   UNIT25     - Ocean land-sea mask on gaussian grid         
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>
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.
      include 'netcdf.inc'
      logical fexist, opened
      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,nm,nr,nf0,nf1,efac,blat,nw
      fsize=65536
      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
!      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.
      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 grid resolution if the OUTGRID exists.
      if( trim(outgrid) .NE. "none" ) then
         inquire(file=trim(outgrid), exist=fexist)
         if(.not. fexist) then
            print*, "file "//trim(outgrid)//" does not exist"
            CALL errexit(4)
         endif
         do ncid = 103, 512
           inquire( ncid,opened=opened )
           if( .NOT.opened )exit
         end do

         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)
         call netcdf_err(error, 'inquire dimension nx from file '//
     &                   trim(outgrid) )
         error=nf_inq_dimlen(ncid,id_dim,nx)
         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) )
         print*, "nx = ", nx
         if(im .ne. nx/2) then
            print*, "IM=",im, " /= grid file nx/2=",nx/2
            print*, "Set IM = ", nx/2
            im = nx/2
         endif
         if(jm .ne. ny/2) then
            print*, "JM=",jm, " /= grid file ny/2=",ny/2
            print*, "Set JM = ", ny/2
            jm = ny/2
         endif
         error=nf_close(ncid)
         call netcdf_err(error, 'close file '//trim(outgrid) )
         
      endif         
         
      CALL tersub(imn,jmn,im,jm,nm,nr,nf0,nf1,nw,efac,blat,
     &            outgrid,inputorog,mask_only,merge_file)
      stop
      END

      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                      :: 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

      logical, intent(in) :: mask_only

      real, parameter :: missing_value=-9999.
      real, PARAMETER :: pi=3.1415926535897931
      integer, PARAMETER :: nmt=14

      integer :: efac,blat,zsave1,zsave2
      integer :: mskocn,notocn
      integer :: i,j,nx,ny,ncid,js,jn,iw,ie,k,it,jt,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 :: 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 :: lonsperlat(:)

      integer, allocatable :: ist(:,:),ien(:,:),zslmx(:,:)
      integer, allocatable :: zavg(:,:),zslm(:,:)
      integer(1), allocatable :: umd(:,:)
      integer(2), allocatable :: glob(:,:)

      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, 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(:,:)
      real, allocatable :: slm(:,:),oro(:,:),var(:,:),orf(:,:)
      real, allocatable :: land_frac(:,:),lake_frac(:,:)
      real, allocatable :: theta(:,:),gamma(:,:),sigma(:,:),elvmax(:,:)
      real, allocatable :: var4(:,:),slmi(:,:)
      real, allocatable :: work1(:,:),work2(:,:),work3(:,:),work4(:,:)
      real, allocatable :: work5(:,:),work6(:,:)
      real, allocatable :: tmpvar(:,:)
      real, allocatable :: slm_in(:,:), lon_in(:,:), lat_in(:,:)
      real(4), allocatable::  gice(:,:),oclsm(:,:)

      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 :: 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 (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 (zavg(imn,jmn))
      allocate (zslm(imn,jmn))
      allocate (umd(imn,jmn))

!
!  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
      revlat = blat .LT. 0   ! Reverse latitude/longitude for output
      mskocn = 1             ! Ocean land sea mask =1, =0 if not present
      notocn = 1             ! =1 Ocean lsm input reverse: Ocean=1, land=0 
! --- The LSM Gaussian file from the ocean model sometimes arrives with 
! --- 0=Ocean and 1=Land or it arrives with 1=Ocean and 0=land without 
! --- metadata to distinguish its disposition.  The AI below mitigates this.

      print *,' In TERSUB'
                      if (mskocn .eq. 1)then
      print *,' Ocean Model LSM Present and '
      print *, ' Overrides OCEAN POINTS in LSM: mskocn=',mskocn
                    if (notocn .eq. 1)then
      print *,' Ocean LSM Reversed:  NOTOCN=',notocn
                    endif
                      endif

      print *,' Attempt to open/read UMD 30sec slmsk.'

      error=nf__open("./landcover.umd.30s.nc",nf_nowrite,fsize,ncid)
      call netcdf_err(error, 'Open file landcover.umd.30s.nc' )
      error=nf_inq_varid(ncid, 'land_mask', id_var)
      call netcdf_err(error, 'Inquire varid of land_mask')
      error=nf_get_var_int1(ncid, id_var, umd)
      call netcdf_err(error, 'Inquire data of land_mask')
      error = nf_close(ncid)

      print *,' UMD lake, UMD(50,50)=',umd(50,50)
C
C- READ_G for global 30" terrain 
C
       print *,' Call read_g to read global topography'
       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
        jt=jmn - j + 1
        i2save = glob(i,j)
        glob(i,j)=glob(i,jt)
        glob(i,jt) = i2save
       enddo
       enddo 
! --- transpose glob as USGS 30" is from dateline and NCEP std is 0
       do j=1,jmn
       do i=1,imn/2
        it=imn/2 + i 
        i2save = glob(i,j)
        glob(i,j)=glob(it,j)
        glob(it,j) = i2save
       enddo
       enddo 
       print *,' After read_g, glob(500,500)=',glob(500,500)
!

!  --- IMN,JMN
      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)
       print *,' UBOUND ZSLM=',ubound(zslm)
       print *,' UBOUND GICE=',imn+1,3601
       print *,' UBOUND OCLSM=',im,jm
!
! ---  0 is ocean and 1 is land for slm
!
C
! --- ZSLM initialize with all land 1, ocean 0
      zslm=1
! --- ZAVG initialize from glob
      zavg=glob

! --- transpose mask even though glob 30" is from N to S and NCEP std is S to N
      do j=1,jmn/2
      do i=1,imn
        jt=jmn - j + 1
        i3save = umd(i,j)
        umd(i,j)=umd(i,jt)
        umd(i,jt) = i3save
      enddo
      enddo 
! --- transpose UMD as USGS 30" is from dateline and NCEP std is 0
      do j=1,jmn
      do i=1,imn/2
        it=imn/2 + i 
        i3save = umd(i,j)
        umd(i,j)=umd(it,j)
        umd(it,j) = i3save
      enddo
      enddo
! ---  Non-land is 0.
      do j=1,jmn
      do i=1,imn
         if ( umd(i,j) .eq. 0 ) zslm(i,j) = 0
      enddo
      enddo

      deallocate (zslmx,umd,glob)
! ---
! ---  Fixing an error in the topo 30" data set at pole (-9999).  
            do i=1,imn
            zslm(i,1)=0
            zslm(i,jmn)=1
            enddo
!
!      print *,' kount1,2,ZAVG(1,1),ZAVG(imn,jmn),ZAVG(500,500)',
!     & kount,kount2,ZAVG(1,1),ZAVG(imn,jmn),ZAVG(500,500)
! --- The center of pixel (1,1) is 89.9958333N/179.9958333W with dx/dy 
! --- spacing of 1/120 degrees. 
!
!  READ REDUCED GRID EXTENTS IF GIVEN
!
      read(20,*,iostat=ios) latg2,lonsperlat
      if(ios.ne.0.or.2*latg2.ne.jm) then
        do j=1,jm
          numi(j)=im
        enddo
        print *,ios,latg2,'COMPUTE TERRAIN ON A FULL GAUSSIAN GRID'
      else
        do j=1,jm/2
          numi(j)=lonsperlat(j)
        enddo
        do j=jm/2+1,jm
          numi(j)=lonsperlat(jm+1-j)
        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

      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
         glat(j) = -90. + (j-1) * delxn + delxn * 0.5
      ENDDO
      print *,
     & ' Before GICE ZAVG(1,2)=',zavg(1,2),zslm(1,2)
      print *,
     & ' Before GICE ZAVG(1,12)=',zavg(1,12),zslm(1,12)
      print *,
     & ' Before GICE ZAVG(1,52)=',zavg(1,52),zslm(1,52)
      print *,
     & ' Before GICE ZAVG(1,112)=',zavg(1,jmn-112),zslm(1,112)

! Read 30-sec Antarctica RAMP data. Points scan from South
! to North, and from Greenwich to Greenwich.

! The error handling here needs to be cleaned up.
      iosg = 0
      error=nf__open("./topography.antarctica.ramp.30s.nc",
     &               nf_nowrite,fsize,ncid)
      call netcdf_err(error, 'Opening RAMP topo file' )
      error=nf_inq_varid(ncid, 'topo', id_var)
      call netcdf_err(error, 'Inquire varid of RAMP topo')
      error=nf_get_var_real(ncid, id_var, gice)
      iosg=error
      call netcdf_err(error, 'Inquire data of RAMP topo')
      error = nf_close(ncid)
      
      if(iosg .ne. 0 ) then
        print *,' *** Err on reading GICE record, iosg=',iosg
        print *,' exec continues but NO GICE correction done '
      else
        print *,' GICE 30" Antarctica RAMP orog 43201x3601 read OK'
        print *,' Processing! '
        print *,' Processing! '
        print *,' Processing! '
         do j = 1, 3601 
         do i = 1, imn
           zsave1 = zavg(i,j)
           zsave2 = zslm(i,j)
         if( gice(i,j) .ne. -99. .and.  gice(i,j) .ne. -1.0 ) then
           if ( gice(i,j) .gt. 0.) then 
                zavg(i,j) = int( gice(i,j) + 0.5 )
!! --- for GICE values less than or equal to 0 (0, -1, or -99) then
!! --- radar-sat (RAMP) values are not valid and revert back to old orog 
                zslm(i,j) = 1
           endif
         endif
  152 format(1x,' ZAVG(i=',i4,' j=',i4,')=',i5,i3,
     &' orig:',i5,i4,' Lat=',f7.3,f8.2,'E',' GICE=',f8.1)
         enddo
         enddo
      endif

      deallocate (gice)

      allocate (oclsm(im,jm),slmi(im,jm))
!C
C     COMPUTE MOUNTAIN DATA : ORO SLM VAR (Std Dev) OC
C
! --- The coupled ocean model is already on a Guasian grid if (IM,JM)
! --- Attempt to Open the file if mskocn=1
      istat=0
          if (mskocn .eq. 1) then
!     open(25,form='unformatted',iostat=istat)
!     open(25,form='binary',iostat=istat)
! --- open to fort.25 with link to file in script
      open(25,form='formatted',iostat=istat)
                  if (istat.ne.0) then
                  mskocn = 0
      print *,' Ocean lsm file Open failure: mskocn,istat=',mskocn,istat
                  else
                  mskocn = 1
      print *,' Ocean lsm file Opened OK: mskocn,istat=',mskocn,istat
                  endif
! --- Read it in
      ios=0
      oclsm=0.
!      read(25,iostat=ios)OCLSM
       read(25,*,iostat=ios)oclsm
         if (ios.ne.0) then 
         mskocn = 0
! --- did not properly read Gaussian grid ocean land-sea mask, but
!     continue using ZSLMX 
      print *,' Rd fail: Ocean lsm - continue, mskocn,ios=',mskocn,ios
         else
         mskocn = 1
      print *,' Rd OK: ocean lsm:  mskocn,ios=',mskocn,ios
! --- LSM initialized to ocean mask especially for case where Ocean
! --- changed by ocean model to land to cope with its problems
! --- remember, that lake mask is in zslm to be assigned in MAKEMT.
          if ( mskocn .eq. 1 ) then
      DO j = 1,jm
      DO i = 1,numi(j)
            if ( notocn .eq. 0 ) then
            slmi(i,j) = float(nint(oclsm(i,j)))
            else
                     if ( nint(oclsm(i,j)) .eq. 0) then
                     slmi(i,j) = 1
                     else 
                     slmi(i,j) = 0
                     endif
            endif
      enddo
      enddo
      print *,' OCLSM',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
      print *,' SLMI:',slmi(1,1),slmi(50,50),slmi(75,75),slmi(im,jm)
! --- Diag
!        WRITE(27,iostat=ios) REAL(SLMI,4)
!        print *,' write SLMI/OCLSM diag input:',ios
          endif
         endif

          else
          print *,' Not using Ocean model land sea mask'
          endif

      if (mskocn .eq. 1)then
      print *,' LSM:',oclsm(1,1),oclsm(50,50),oclsm(75,75),oclsm(im,jm)
      endif

      allocate (geolon(im,jm),geolon_c(im+1,jm+1),dx(im,jm))
      allocate (geolat(im,jm),geolat_c(im+1,jm+1),dy(im,jm))
      allocate (slm(im,jm),oro(im,jm),var(im,jm),var4(im,jm))
      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.
         inquire(file=trim(outgrid), exist=fexist)
         if(.not. fexist) then
            print*, "file "//trim(outgrid)//" does not exist"
            CALL errexit(4)
         endif
         do ncid = 103, 512
           inquire( ncid,opened=opened )
           if( .NOT.opened )exit
         end do

         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)
         call netcdf_err(error, 'inquire dimension nx from file '//
     &                   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))

         error=nf_inq_varid(ncid, 'x', id_var)
         call netcdf_err(error, 'inquire varid of x from file '
     &                   //trim(outgrid) )
         error=nf_get_var_double(ncid, id_var, tmpvar)
         call netcdf_err(error, 'inquire data of x from file '
     &                   //trim(outgrid) )
         !--- adjust lontitude to be between 0 and 360.
         do j = 1,ny+1; do i = 1,nx+1
            if(tmpvar(i,j) .NE. missing_value) then
              if(tmpvar(i,j) .GT. 360) tmpvar(i,j) = tmpvar(i,j) - 360
              if(tmpvar(i,j) .LT. 0) tmpvar(i,j) = tmpvar(i,j) + 360
            endif
         enddo; enddo

         geolon(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
         geolon_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)
         
         error=nf_inq_varid(ncid, 'y', id_var)
         call netcdf_err(error, 'inquire varid of y from file '
     &                   //trim(outgrid) )
         error=nf_get_var_double(ncid, id_var, tmpvar)
         call netcdf_err(error, 'inquire data of y from file '
     &                   //trim(outgrid) )
         geolat(1:im,1:jm) = tmpvar(2:nx:2,2:ny:2)
         geolat_c(1:im+1,1:jm+1) = tmpvar(1:nx+1:2,1:ny+1:2)

         !--- figure out pole location.
         maxlat = -90
         minlat = 90
         i_north_pole = 0
         j_north_pole = 0
         i_south_pole = 0
         j_south_pole = 0
         do j = 1, ny+1; do i = 1, nx+1
            if( tmpvar(i,j) > maxlat ) then
               i_north_pole=i
               j_north_pole=j
               maxlat = tmpvar(i,j)
            endif
            if( tmpvar(i,j) < minlat ) then
               i_south_pole=i
               j_south_pole=j
               minlat = tmpvar(i,j)
            endif
         enddo ; enddo
         !--- only when maxlat is close to 90. the point is north pole
         if(maxlat < 89.9 ) then
            i_north_pole = 0
            j_north_pole = 0
         endif
         if(minlat > -89.9 ) then
            i_south_pole = 0
            j_south_pole = 0
         endif
         print*, "minlat=", minlat, "maxlat=", maxlat
         print*, "north pole supergrid index is ",
     &           i_north_pole, j_north_pole
         print*, "south pole supergrid index is ",
     &           i_south_pole, j_south_pole
         deallocate(tmpvar)

         if(i_south_pole >0 .and. j_south_pole > 0) then
           if(mod(i_south_pole,2)==0) then ! stretched grid
             do j = 1, jm; do i = 1, im
               if(i==i_south_pole/2 .and. (j==j_south_pole/2 
     &              .or. j==j_south_pole/2+1) ) then
                 is_south_pole(i,j) = .true.
                 print*, "south pole at i,j=", i, j
               endif
             enddo; enddo      
            else
               do j = 1, jm; do i = 1, im
                  if((i==i_south_pole/2 .or. i==i_south_pole/2+1)
     &             .and. (j==j_south_pole/2 .or.
     &                j==j_south_pole/2+1) ) then
                    is_south_pole(i,j) = .true.
                    print*, "south pole at i,j=", i, j
                  endif
               enddo; enddo
            endif            
         endif
         if(i_north_pole >0 .and. j_north_pole > 0) then
            if(mod(i_north_pole,2)==0) then ! stretched grid
               do j = 1, jm; do i = 1, im
                  if(i==i_north_pole/2 .and. (j==j_north_pole/2 .or.
     &                j==j_north_pole/2+1) ) then
                    is_north_pole(i,j) = .true.
                    print*, "north pole at i,j=", i, j
                  endif
               enddo; enddo      
            else
               do j = 1, jm; do i = 1, im
                  if((i==i_north_pole/2 .or. i==i_north_pole/2+1)
     &             .and. (j==j_north_pole/2 .or.
     &                j==j_north_pole/2+1) ) then
                    is_north_pole(i,j) = .true.
                    print*, "north pole at i,j=", i, j
                  endif
               enddo; enddo
            endif            
         endif
         
         
         allocate(tmpvar(nx,ny))
         error=nf_inq_varid(ncid, 'area', id_var)
         call netcdf_err(error, 'inquire varid of area from file '
     &                   //trim(outgrid) )
         error=nf_get_var_double(ncid, id_var, tmpvar)
         call netcdf_err(error, 'inquire data of area from file '
     &                   //trim(outgrid) )

         do j = 1, jm
            do i = 1, im
               dx(i,j) = sqrt(tmpvar(2*i-1,2*j-1)+tmpvar(2*i,2*j-1)
     &                       +tmpvar(2*i-1,2*j  )+tmpvar(2*i,2*j  ))
               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
           CALL make_mask(zslm,slm,land_frac,glat,
     &      im,jm,imn,jmn,geolon_c,geolat_c)
            lake_frac=9999.9
         ELSE
           print*,'got here - read in external mask ',merge_file
           CALL read_mask(merge_file,slm,land_frac,lake_frac,im,jm)
         ENDIF

         IF (mask_only) THEN
           print*,'got here computing mask only.'
           CALL write_mask_netcdf(im,jm,slm,land_frac,
     1                  1,1,geolon,geolat)

           print*,' DONE.'
           stop
         END IF

         CALL makemt2(zavg,zslm,oro,slm,var,var4,glat,
     & im,jm,imn,jmn,geolon_c,geolat_c,lake_frac,land_frac)

      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')
       call minmxj(im,jm,var,'     VAR')
       call minmxj(im,jm,var4,'    VAR4')
C --- check for nands in above
!      call nanc(ORO,IM*JM,"MAKEMT_ORO")
!      call nanc(SLM,IM*JM,"MAKEMT_SLM")
!      call nanc(VAR,IM*JM,"MAKEMT_VAR")
!      call nanc(VAR4,IM*JM,"MAKEMT_VAR4")
!
C   check antarctic pole 
!     DO J = 1,JM
!     DO I = 1,numi(j)
!        if ( i .le. 100 .and. i .ge. 1 )then
!           if (j .ge. JM-1 )then
!      if (height .eq. 0.) print *,'I,J,SLM:',I,J,SLM(I,J)
!      write(6,153)i,j,ORO(i,j),HEIGHT,SLM(i,j)
!             endif
!        endif    
!     ENDDO
!     ENDDO
C
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')
       call minmxj(im,jm,sigma,'   SIGMA')

C --- check for nands in above
!      call nanc(THETA,IM*JM,"MAKEPC_THETA")
!      call nanc(GAMMA,IM*JM,"MAKEPC_GAMMA")
!      call nanc(SIGMA,IM*JM,"MAKEPC_SIGMA")
C
C     COMPUTE MOUNTAIN DATA : OA OL
C
       allocate (iwork(im,jm,4))
       allocate (oa(im,jm,4),ol(im,jm,4),hprime(im,jm,14))
       allocate (work1(im,jm),work2(im,jm),work3(im,jm),work4(im,jm))
       allocate (work5(im,jm),work6(im,jm))

       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()
           CALL makeoa2(zavg,zslm,var,glat,oa,ol,iwork,elvmax,oro,
     1            work1,work2,work3,work4,work5,work6,
     2            im,jm,imn,jmn,geolon_c,geolat_c,
     3            geolon,geolat,dx,dy,is_south_pole,is_north_pole)
      tend=timef()
      write(6,*)' MAKEOA2 time= ',tend-tbeg
         else
           !-- read the data from INPUTOROG file.
           error=nf__open(trim(inputorog),nf_nowrite,fsize,ncid)
           call netcdf_err(error, 'Open file '//trim(inputorog) )
           error=nf_inq_dimid(ncid, 'lon', id_dim)
           call netcdf_err(error, 'inquire dimension lon from file '//
     &                   trim(inputorog) )
           error=nf_inq_dimlen(ncid,id_dim,nx_in)
           call netcdf_err(error, 'inquire dimension lon length '//
     &       'from file '//trim(inputorog) )
           error=nf_inq_dimid(ncid, 'lat', id_dim)
           call netcdf_err(error, 'inquire dimension lat from file '//
     &                   trim(inputorog) )
           error=nf_inq_dimlen(ncid,id_dim,ny_in)
           call netcdf_err(error, 'inquire dimension lat length '//
     &       'from file '//trim(inputorog) )
           
           print*, "extrapolate OA, OL from Gaussian grid with nx=",
     &              nx_in, ", ny=", ny_in
           allocate(oa_in(nx_in,ny_in,4), ol_in(nx_in,ny_in,4))
           allocate(slm_in(nx_in,ny_in) )
           allocate(lon_in(nx_in,ny_in), lat_in(nx_in,ny_in) )
           
           error=nf_inq_varid(ncid, 'oa1', id_var)
           call netcdf_err(error, 'inquire varid of oa1 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, oa_in(:,:,1))
           call netcdf_err(error, 'inquire data of oa1 from file '
     &                   //trim(inputorog) )
           error=nf_inq_varid(ncid, 'oa2', id_var)
           call netcdf_err(error, 'inquire varid of oa2 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, oa_in(:,:,2))
           call netcdf_err(error, 'inquire data of oa2 from file '
     &                   //trim(inputorog) )
           error=nf_inq_varid(ncid, 'oa3', id_var)
           call netcdf_err(error, 'inquire varid of oa3 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, oa_in(:,:,3))
           call netcdf_err(error, 'inquire data of oa3 from file '
     &                   //trim(inputorog) )           
           error=nf_inq_varid(ncid, 'oa4', id_var)
           call netcdf_err(error, 'inquire varid of oa4 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, oa_in(:,:,4))
           call netcdf_err(error, 'inquire data of oa4 from file '
     &                   //trim(inputorog) )

           error=nf_inq_varid(ncid, 'ol1', id_var)
           call netcdf_err(error, 'inquire varid of ol1 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, ol_in(:,:,1))
           call netcdf_err(error, 'inquire data of ol1 from file '
     &                   //trim(inputorog) )
           error=nf_inq_varid(ncid, 'ol2', id_var)
           call netcdf_err(error, 'inquire varid of ol2 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, ol_in(:,:,2))
           call netcdf_err(error, 'inquire data of ol2 from file '
     &                   //trim(inputorog) )
           error=nf_inq_varid(ncid, 'ol3', id_var)
           call netcdf_err(error, 'inquire varid of ol3 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, ol_in(:,:,3))
           call netcdf_err(error, 'inquire data of ol3 from file '
     &                   //trim(inputorog) )           
           error=nf_inq_varid(ncid, 'ol4', id_var)
           call netcdf_err(error, 'inquire varid of ol4 from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, ol_in(:,:,4))
           call netcdf_err(error, 'inquire data of ol4 from file '
     &                   //trim(inputorog) )

           error=nf_inq_varid(ncid, 'slmsk', id_var)
           call netcdf_err(error, 'inquire varid of slmsk from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, slm_in)
           call netcdf_err(error, 'inquire data of slmsk from file '
     &                   //trim(inputorog) )

           error=nf_inq_varid(ncid, 'geolon', id_var)
           call netcdf_err(error, 'inquire varid of geolon from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, lon_in)
           call netcdf_err(error, 'inquire data of geolon from file '
     &                   //trim(inputorog) )

           error=nf_inq_varid(ncid, 'geolat', id_var)
           call netcdf_err(error, 'inquire varid of geolat from file '
     &                   //trim(inputorog) )
           error=nf_get_var_double(ncid, id_var, lat_in)
           call netcdf_err(error, 'inquire data of geolat from file '
     &                   //trim(inputorog) )
           
           ! set slmsk=2 to be ocean (0)
           do j=1,ny_in; do i=1,nx_in
              if(slm_in(i,j) == 2) slm_in(i,j) = 0
           enddo; enddo
           
           error=nf_close(ncid)
           call netcdf_err(error, 'close file '
     &                   //trim(inputorog) )
           
           print*, "calling MAKEOA3 to compute OA, OL"
           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,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)
       deallocate (zslm,zavg)
       deallocate (dx,dy)
       deallocate (work2,work3,work4,work5,work6)

!      Deallocate 3d vars
       deallocate(iwork)

       tbeg=timef()
       call minmxj(im,jm,oa,'      OA')
       call minmxj(im,jm,ol,'      OL')
       call minmxj(im,jm,elvmax,'  ELVMAX')
       call minmxj(im,jm,oro,'     ORO')
C --- check for nands in above
!      call nanc(OA(1,1,1), IM*JM,"MAKEOA_OA(1,1,1)")
!      call nanc(OA(1,1,2), IM*JM,"MAKEOA_OA(1,1,2)")
!      call nanc(OA(1,1,3), IM*JM,"MAKEOA_OA(1,1,3)")
!      call nanc(OA(1,1,4), IM*JM,"MAKEOA_OA(1,1,4)")
!      call nanc(OL(1,1,1), IM*JM,"MAKEOA_OL(1,1,1)")
!      call nanc(OL(1,1,2), IM*JM,"MAKEOA_OL(1,1,2)")
!      call nanc(OL(1,1,3), IM*JM,"MAKEOA_OL(1,1,3)")
!      call nanc(OL(1,1,4), IM*JM,"MAKEOA_OL(1,1,4)")
!      call nanc(ELVMAX, IM*JM,"MAKEPC_ELVMAX")

        maxc3 = 0
        maxc4 = 0
        maxc5 = 0
        maxc6 = 0
        maxc7 = 0
        maxc8 = 0
      DO j = 1,jm
      DO i = 1,numi(j)
         if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
         if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
         if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
         if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
         if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
         if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
      ENDDO
      ENDDO
      print *,' MAXC3:',maxc3,maxc4,maxc5,maxc6,maxc7,maxc8
!
c      itest=151
c      jtest=56
C      
       print *,' ===> Replacing ELVMAX with ELVMAX-ORO <=== ' 
       print *,' ===> if ELVMAX<=ORO replace with proxy <=== ' 
       print *,' ===> the sum of mean orog (ORO) and std dev <=== ' 
      DO j = 1,jm
      DO i = 1,numi(j)
        if (elvmax(i,j) .lt. oro(i,j) ) then
C---  subtracting off ORO leaves std dev (this should never happen)
       elvmax(i,j) = max(  3. * var(i,j),0.)
        else
       elvmax(i,j) = max( elvmax(i,j) - oro(i,j),0.)
        endif
      ENDDO
      ENDDO
        maxc3 = 0
        maxc4 = 0
        maxc5 = 0
        maxc6 = 0
        maxc7 = 0
        maxc8 = 0
      DO j = 1,jm
      DO i = 1,numi(j)
         if (elvmax(i,j) .gt. 3000.) maxc3 = maxc3 +1
         if (elvmax(i,j) .gt. 4000.) maxc4 = maxc4 +1
         if (elvmax(i,j) .gt. 5000.) maxc5 = maxc5 +1
         if (elvmax(i,j) .gt. 6000.) maxc6 = maxc6 +1
         if (elvmax(i,j) .gt. 7000.) maxc7 = maxc7 +1
         if (elvmax(i,j) .gt. 8000.) maxc8 = maxc8 +1
      ENDDO
      ENDDO
      print *,' after MAXC 3-6 km:',maxc3,maxc4,maxc5,maxc6
c
       call mnmxja(im,jm,elvmax,itest,jtest,'  ELVMAX')
!     if (JM .gt. 0) stop
C
C     ZERO OVER OCEAN
C
      print *,' Testing at point (itest,jtest)=',itest,jtest
      print *,' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
      print *,' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest
      DO j = 1,jm
        DO i = 1,numi(j)
          IF(slm(i,j).EQ.0.) THEN
C           VAR(I,J) = 0.
            var4(i,j) = 0.
            oa(i,j,1) = 0.
            oa(i,j,2) = 0.
            oa(i,j,3) = 0.
            oa(i,j,4) = 0.
            ol(i,j,1) = 0.
            ol(i,j,2) = 0.
            ol(i,j,3) = 0.
            ol(i,j,4) = 0.
C           THETA(I,J) =0.
C           GAMMA(I,J) =0.
C           SIGMA(I,J) =0.
C           ELVMAX(I,J)=0.
! --- the sub-grid scale parameters for mtn blocking and gwd retain 
! --- properties even if over ocean but there is elevation within the 
! --- gaussian grid box.
          ENDIF
       ENDDO
      ENDDO
C
! --- if mskocn=1 ocean land sea mask given, =0 if not present
! ---  OCLSM is real(*4)  array with fractional values possible
! ---  0 is ocean and 1 is land for slm
! ---  Step 1: Only change SLM after GFS SLM is applied
! ---  SLM is only field that will be altered by OCLSM
! ---  Ocean land sea mask ocean points made ocean in atm model
! ---  Land and Lakes and all other atm elv moments remain unchanged.  

      IF (merge_file == 'none') then

      msk_ocn :      if ( mskocn .eq. 1 ) then

      DO j = 1,jm
        DO i = 1,numi(j)
              if (abs(oro(i,j)) .lt. 1. ) then 
              slm(i,j) = slmi(i,j)
              else
      if ( slmi(i,j) .eq. 1. .and. slm(i,j) .eq. 1) slm(i,j) = 1
      if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
      if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 1) slm(i,j) = 0
      if ( slmi(i,j) .eq. 0. .and. slm(i,j) .eq. 0) slm(i,j) = 0
              endif
        enddo
      enddo
            endif msk_ocn
      endif
      print *,' SLM(itest,jtest)=',slm(itest,jtest),itest,jtest
      print *,' ORO(itest,jtest)=',oro(itest,jtest),itest,jtest

      deallocate(slmi)

      IF (merge_file == 'none') then

C  REMOVE ISOLATED POINTS
      iso_loop : DO j=2,jm-1
        jn=j-1
        js=j+1
        rn=REAL(numi(jn))/REAL(numi(j))
        rs=REAL(numi(js))/REAL(numi(j))
        DO i=1,numi(j)
          iw=mod(i+im-2,im)+1
          ie=mod(i,im)+1
          slma=slm(iw,j)+slm(ie,j)
          oroa=oro(iw,j)+oro(ie,j)
          vara=var(iw,j)+var(ie,j)
          var4a=var4(iw,j)+var4(ie,j)
          DO k=1,4
            oaa(k)=oa(iw,j,k)+oa(ie,j,k)
! --- (*j*) fix typo:
            ola(k)=ol(iw,j,k)+ol(ie,j,k)
          ENDDO
          wgta=2
          xn=rn*(i-1)+1
          IF(abs(xn-nint(xn)).LT.1.e-2) THEN
            in=mod(nint(xn)-1,numi(jn))+1
            inw=mod(in+numi(jn)-2,numi(jn))+1
            ine=mod(in,numi(jn))+1
            slma=slma+slm(inw,jn)+slm(in,jn)+slm(ine,jn)
            oroa=oroa+oro(inw,jn)+oro(in,jn)+oro(ine,jn)
            vara=vara+var(inw,jn)+var(in,jn)+var(ine,jn)
            var4a=var4a+var4(inw,jn)+var4(in,jn)+var4(ine,jn)
            DO k=1,4
              oaa(k)=oaa(k)+oa(inw,jn,k)+oa(in,jn,k)+oa(ine,jn,k)
              ola(k)=ola(k)+ol(inw,jn,k)+ol(in,jn,k)+ol(ine,jn,k)
            ENDDO
            wgta=wgta+3
          ELSE
            inw=int(xn)
            ine=mod(inw,numi(jn))+1
            slma=slma+slm(inw,jn)+slm(ine,jn)
            oroa=oroa+oro(inw,jn)+oro(ine,jn)
            vara=vara+var(inw,jn)+var(ine,jn)
            var4a=var4a+var4(inw,jn)+var4(ine,jn)
            DO k=1,4
              oaa(k)=oaa(k)+oa(inw,jn,k)+oa(ine,jn,k)
              ola(k)=ola(k)+ol(inw,jn,k)+ol(ine,jn,k)
            ENDDO
            wgta=wgta+2
          ENDIF
          xs=rs*(i-1)+1
          IF(abs(xs-nint(xs)).LT.1.e-2) THEN
            is=mod(nint(xs)-1,numi(js))+1
            isw=mod(is+numi(js)-2,numi(js))+1
            ise=mod(is,numi(js))+1
            slma=slma+slm(isw,js)+slm(is,js)+slm(ise,js)
            oroa=oroa+oro(isw,js)+oro(is,js)+oro(ise,js)
            vara=vara+var(isw,js)+var(is,js)+var(ise,js)
            var4a=var4a+var4(isw,js)+var4(is,js)+var4(ise,js)
            DO k=1,4
              oaa(k)=oaa(k)+oa(isw,js,k)+oa(is,js,k)+oa(ise,js,k)
              ola(k)=ola(k)+ol(isw,js,k)+ol(is,js,k)+ol(ise,js,k)
            ENDDO
            wgta=wgta+3
          ELSE
            isw=int(xs)
            ise=mod(isw,numi(js))+1
            slma=slma+slm(isw,js)+slm(ise,js)
            oroa=oroa+oro(isw,js)+oro(ise,js)
            vara=vara+var(isw,js)+var(ise,js)
            var4a=var4a+var4(isw,js)+var4(ise,js)
            DO k=1,4
              oaa(k)=oaa(k)+oa(isw,js,k)+oa(ise,js,k)
              ola(k)=ola(k)+ol(isw,js,k)+ol(ise,js,k)
            ENDDO
            wgta=wgta+2
          ENDIF
          oroa=oroa/wgta
          vara=vara/wgta
          var4a=var4a/wgta
          DO k=1,4
            oaa(k)=oaa(k)/wgta
            ola(k)=ola(k)/wgta
          ENDDO
          IF(slm(i,j).EQ.0..AND.slma.EQ.wgta) THEN
            print '("SEA ",2F8.0," MODIFIED TO LAND",2F8.0," AT ",2I8)',
     &       oro(i,j),var(i,j),oroa,vara,i,j
            slm(i,j)=1.
            oro(i,j)=oroa
            var(i,j)=vara
            var4(i,j)=var4a
            DO k=1,4
              oa(i,j,k)=oaa(k)
              ol(i,j,k)=ola(k)
            ENDDO
          ELSEIF(slm(i,j).EQ.1..AND.slma.EQ.0.) THEN
            print '("LAND",2F8.0," MODIFIED TO SEA ",2F8.0," AT ",2I8)',
     &       oro(i,j),var(i,j),oroa,vara,i,j
            slm(i,j)=0.
            oro(i,j)=oroa
            var(i,j)=vara
            var4(i,j)=var4a
            DO k=1,4
              oa(i,j,k)=oaa(k)
              ol(i,j,k)=ola(k)
            ENDDO
          ENDIF
        ENDDO
      ENDDO iso_loop
C--- print for testing after isolated points removed
      print *,' after isolated points removed'
       call minmxj(im,jm,oro,'     ORO')
C     print *,' JM=',JM,' numi=',numi
      print *,' ORO(itest,jtest)=',oro(itest,jtest)
      print *,' VAR(itest,jtest)=',var(itest,jtest)
      print *,' VAR4(itest,jtest)=',var4(itest,jtest)
      print *,' OA(itest,jtest,1)=',oa(itest,jtest,1)
      print *,' OA(itest,jtest,2)=',oa(itest,jtest,2)
      print *,' OA(itest,jtest,3)=',oa(itest,jtest,3)
      print *,' OA(itest,jtest,4)=',oa(itest,jtest,4)
      print *,' OL(itest,jtest,1)=',ol(itest,jtest,1)
      print *,' OL(itest,jtest,2)=',ol(itest,jtest,2)
      print *,' OL(itest,jtest,3)=',ol(itest,jtest,3)
      print *,' OL(itest,jtest,4)=',ol(itest,jtest,4)
      print *,' Testing at point (itest,jtest)=',itest,jtest
      print *,' THETA(itest,jtest)=',theta(itest,jtest)
      print *,' GAMMA(itest,jtest)=',gamma(itest,jtest)
      print *,' SIGMA(itest,jtest)=',sigma(itest,jtest)
      print *,' ELVMAX(itest,jtest)=',elvmax(itest,jtest)
      print *,' EFAC=',efac

      endif

C
      DO j=1,jm
        DO i=1,numi(j)
          oro(i,j) = oro(i,j) + efac*var(i,j)
          hprime(i,j,1) = var(i,j)
          hprime(i,j,2) = var4(i,j)
          hprime(i,j,3) = oa(i,j,1)
          hprime(i,j,4) = oa(i,j,2)
          hprime(i,j,5) = oa(i,j,3)
          hprime(i,j,6) = oa(i,j,4)
          hprime(i,j,7) = ol(i,j,1)
          hprime(i,j,8) = ol(i,j,2)
          hprime(i,j,9) = ol(i,j,3)
          hprime(i,j,10)= ol(i,j,4)
          hprime(i,j,11)= theta(i,j)
          hprime(i,j,12)= gamma(i,j)
          hprime(i,j,13)= sigma(i,j)
          hprime(i,j,14)= elvmax(i,j)
        ENDDO
      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
      ENDIF

      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
C  USE NEAREST NEIGHBOR INTERPOLATION TO FILL FULL GRIDS
      call rg2gg(im,jm,numi,slm)
      call rg2gg(im,jm,numi,oro)
      call rg2gg(im,jm,numi,orf)
C ---   not apply to new prin coord and ELVMAX (*j*)
      do imt=1,10 
        call rg2gg(im,jm,numi,hprime(1,1,imt))
      enddo
C 
      print *,' after nearest neighbor interpolation applied '
       call minmxj(im,jm,oro,'     ORO')
       call minmxj(im,jm,orf,'     ORF')
       call mnmxja(im,jm,elvmax,itest,jtest,'  ELVMAX')
       print *,' ORO,ORF(itest,jtest),itest,jtest:',
     &          oro(itest,jtest),orf(itest,jtest),itest,jtest
      print *,' ELVMAX(',itest,jtest,')=',elvmax(itest,jtest)


C   check antarctic pole 
      DO j = 1,jm
      DO i = 1,numi(j)
         if ( i .le. 21 .and. i .ge. 1 )then
         if (j .eq. jm )write(6,153)i,j,oro(i,j),elvmax(i,j),slm(i,j)
  153 format(1x,' ORO,ELVMAX(i=',i4,' j=',i4,')=',2e14.5,f5.1)
         endif
      ENDDO
      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)
      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)
      tend=timef()
      write(6,*)' WRITE_NETCDF time= ',tend-tbeg
      print *,' wrote netcdf file out.oro.tile?.nc'

      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 2d vars
      deallocate (oclsm)
      deallocate (geolon,geolon_c,geolat,geolat_c)
      deallocate (slm,oro,var,orf,land_frac)
      deallocate (theta,gamma,sigma,elvmax)


      tend=timef()
      write(6,*)' Total runtime time= ',tend-tbeg1
      RETURN
      END

      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)
      LOGICAL flag, debug
C==== DATA DEBUG/.TRUE./
      DATA debug/.false./
C
! ---- OCLSM holds the ocean (im,jm) grid
      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

      SUBROUTINE get_index(IMN,JMN,npts,lonO,latO,DELXN,
     &           jst,jen,ilist,numx)
        implicit none
        integer, intent(in)  :: imn,jmn
        integer              :: npts
        real,    intent(in)  :: lono(npts), lato(npts)
        real,    intent(in)  :: delxn
        integer, intent(out) :: jst,jen
        integer, intent(out) :: ilist(imn)
        integer, intent(out) :: numx
        real    minlat,maxlat,minlon,maxlon
        integer i2, ii, ist, ien
        
         minlat = minval(lato)
         maxlat = maxval(lato)
         minlon = minval(lono)
         maxlon = maxval(lono)
         ist = minlon/delxn+1
         ien = maxlon/delxn+1
         jst = (minlat+90)/delxn+1 
         jen = (maxlat+90)/delxn 
         !--- add a few points to both ends of j-direction
         jst = jst - 5
         if(jst<1) jst = 1
         jen = jen + 5
         if(jen>jmn) jen = jmn

         !--- when around the pole, just search through all the points.
         if((jst == 1 .OR. jen == jmn) .and. 
     &            (ien-ist+1 > imn/2) )then
            numx = imn
            do i2 = 1, imn
               ilist(i2) = i2
            enddo
         else if( ien-ist+1 > imn/2 ) then  ! cross longitude = 0
            !--- find the minimum that greater than IMN/2 
            !--- and maximum that less than IMN/2
            ist = 0
            ien = imn+1
            do i2 = 1, npts
               ii = lono(i2)/delxn+1
               if(ii <0 .or. ii>imn) print*,"ii=",ii,imn,lono(i2),delxn
               if( ii < imn/2 ) then
                  ist = max(ist,ii)
               else if( ii > imn/2 ) then
                  ien = min(ien,ii)
               endif
            enddo 
            if(ist<1 .OR. ist>imn) then
               print*, .or."ist<1  ist>IMN"
               call abort()
            endif
           if(ien<1 .OR. ien>imn) then
               print*, .or."iend<1  iend>IMN"
               call abort()
            endif

            numx = imn - ien + 1
            do i2 = 1, numx
               ilist(i2) = ien + (i2-1)           
            enddo
            do i2 = 1, ist
               ilist(numx+i2) = i2
            enddo
            numx = numx+ist
         else
            numx = ien-ist+1
            do i2 = 1, numx
               ilist(i2) = ist + (i2-1)
            enddo
         endif

      END   
      
      SUBROUTINE make_mask(zslm,SLM,land_frac,
     1 glat,im,jm,imn,jmn,lon_c,lat_c)
      implicit none
      real, parameter :: d2r = 3.14159265358979/180.
      integer, parameter :: maxsum=20000000
      integer im, jm, imn, jmn, jst, jen
      real glat(jmn), glon(imn)
      INTEGER zslm(imn,jmn)
      real land_frac(im,jm)
      real slm(im,jm)
      real lon_c(im+1,jm+1), lat_c(im+1,jm+1)
      real    lono(4),lato(4),loni,lati
      integer jm1,i,j,nsum,nsum_all,ii,jj,numx,i2
      integer ilist(imn)
      real    delxn,xnsum,xland,xwatr,xl1,xs1,xw1
      real    xnsum_all,xland_all,xwatr_all
      logical inside_a_polygon
C
! ---- OCLSM holds the ocean (im,jm) grid
      print *,' _____ SUBROUTINE MAKE_MASK '
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
      DO i=1,imn
         glon(i) = 0. + (i-1) * delxn + delxn * 0.5
      ENDDO
 
      land_frac(:,:) = 0.0     
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 
!$omp parallel do
!$omp* private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,lono,
!$omp*          lato,jst,jen,ilist,numx,jj,i2,ii,loni,lati,
!$omp*          xnsum_all,xland_all,xwatr_all,nsum_all)
!$omp*          
      DO j=1,jm
!       print*, "J=", J
       DO i=1,im
         xnsum = 0.0
         xland = 0.0
         xwatr = 0.0
         nsum = 0
         xnsum_all = 0.0
         xland_all = 0.0
         xwatr_all = 0.0
         nsum_all = 0
         
         lono(1) = lon_c(i,j) 
         lono(2) = lon_c(i+1,j) 
         lono(3) = lon_c(i+1,j+1) 
         lono(4) = lon_c(i,j+1) 
         lato(1) = lat_c(i,j) 
         lato(2) = lat_c(i+1,j) 
         lato(3) = lat_c(i+1,j+1) 
         lato(4) = lat_c(i,j+1) 
         call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
         do jj = jst, jen; do i2 = 1, numx
            ii = ilist(i2)
            loni = ii*delxn
            lati = -90 + jj*delxn

            xland_all = xland_all + float(zslm(ii,jj))
            xwatr_all = xwatr_all + float(1-zslm(ii,jj))
            xnsum_all = xnsum_all + 1.
            nsum_all = nsum_all+1
            if(nsum_all > maxsum) then
              print*, "nsum_all is greater than MAXSUM, increase MAXSUM"
              call abort()
            endif

            if(inside_a_polygon(loni*d2r,lati*d2r,4,
     &          lono*d2r,lato*d2r))then

               xland = xland + float(zslm(ii,jj))
               xwatr = xwatr + float(1-zslm(ii,jj))
               xnsum = xnsum + 1.
               nsum = nsum+1
               if(nsum > maxsum) then
                 print*, "nsum is greater than MAXSUM, increase MAXSUM"
                 call abort()
               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
               land_frac(i,j) = xland/xnsum  
               slm(i,j) = float(nint(xland/xnsum))
         ELSEIF(xnsum_all.GT.1.) THEN
               land_frac(i,j) = xland_all/xnsum _all 
               slm(i,j) = float(nint(xland_all/xnsum_all))
         ELSE
               print*, "no source points in MAKE_MASK"
               call abort()
         ENDIF
       ENDDO
      ENDDO
!$omp end parallel do
      WRITE(6,*) "! MAKE_MASK DONE"
C
      RETURN
      END SUBROUTINE make_mask
      SUBROUTINE makemt2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4,
     1 glat,im,jm,imn,jmn,lon_c,lat_c,lake_frac,land_frac)
      implicit none
      real, parameter :: d2r = 3.14159265358979/180.
      integer, parameter :: maxsum=20000000
      real, dimension(:), allocatable ::  hgt_1d, hgt_1d_all
      integer im, jm, imn, jmn
      real glat(jmn), glon(imn)
      INTEGER zavg(imn,jmn),zslm(imn,jmn)
      real oro(im,jm),var(im,jm),var4(im,jm)
      real, intent(in) :: slm(im,jm), lake_frac(im,jm),land_frac(im,jm)
      integer jst, jen
      real lon_c(im+1,jm+1), lat_c(im+1,jm+1)
      real    lono(4),lato(4),loni,lati
      real    height
      integer jm1,i,j,nsum,nsum_all,ii,jj,i1,numx,i2
      integer ilist(imn)
      real    delxn,xnsum,xland,xwatr,xl1,xs1,xw1,xw2,xw4
      real    xnsum_all,xland_all,xwatr_all,height_all
      real    xl1_all,xs1_all,xw1_all,xw2_all,xw4_all
      logical inside_a_polygon
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 MAKEMT2 '
      allocate(hgt_1d(maxsum))
      allocate(hgt_1d_all(maxsum))
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
      DO i=1,imn
         glon(i) = 0. + (i-1) * delxn + delxn * 0.5
      ENDDO
 
!     land_frac(:,:) = 0.0     
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 
!$omp parallel do
!$omp* private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,xw2,xw4,lono,
!$omp*          lato,jst,jen,ilist,numx,jj,i2,ii,loni,lati,height,
!$omp*          hgt_1d,
!$omp*          xnsum_all,xland_all,xwatr_all,nsum_all,
!$omp*          xl1_all,xs1_all,xw1_all,xw2_all,xw4_all,
!$omp*          height_all,hgt_1d_all)
      DO j=1,jm
!       print*, "J=", J
       DO i=1,im
         oro(i,j)  = 0.0
         var(i,j)  = 0.0
         var4(i,j) = 0.0
         xnsum = 0.0
         xland = 0.0
         xwatr = 0.0
         nsum = 0
         xl1 = 0.0
         xs1 = 0.0
         xw1 = 0.0
         xw2 = 0.0
         xw4 = 0.0
         xnsum_all = 0.0
         xland_all = 0.0
         xwatr_all = 0.0
         nsum_all = 0
         xl1_all = 0.0
         xs1_all = 0.0
         xw1_all = 0.0
         xw2_all = 0.0
         xw4_all = 0.0
         
         lono(1) = lon_c(i,j) 
         lono(2) = lon_c(i+1,j) 
         lono(3) = lon_c(i+1,j+1) 
         lono(4) = lon_c(i,j+1) 
         lato(1) = lat_c(i,j) 
         lato(2) = lat_c(i+1,j) 
         lato(3) = lat_c(i+1,j+1) 
         lato(4) = lat_c(i,j+1) 
         call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
         do jj = jst, jen; do i2 = 1, numx
            ii = ilist(i2)
            loni = ii*delxn
            lati = -90 + jj*delxn

            xland_all = xland_all + float(zslm(ii,jj))
            xwatr_all = xwatr_all + float(1-zslm(ii,jj))
            xnsum_all = xnsum_all + 1.
            height_all = float(zavg(ii,jj)) 
            nsum_all = nsum_all+1
            if(nsum_all > maxsum) then
              print*, "nsum_all is greater than MAXSUM, increase MAXSUM"
              call abort()
            endif
            hgt_1d_all(nsum_all) = height_all
            IF(height_all.LT.-990.) height_all = 0.0
            xl1_all = xl1_all + height_all * float(zslm(ii,jj))
            xs1_all = xs1_all + height_all * float(1-zslm(ii,jj))
            xw1_all = xw1_all + height_all
            xw2_all = xw2_all + height_all ** 2

            if(inside_a_polygon(loni*d2r,lati*d2r,4,
     &          lono*d2r,lato*d2r))then

               xland = xland + float(zslm(ii,jj))
               xwatr = xwatr + float(1-zslm(ii,jj))
               xnsum = xnsum + 1.
               height = float(zavg(ii,jj)) 
               nsum = nsum+1
               if(nsum > maxsum) then
                 print*, "nsum is greater than MAXSUM, increase MAXSUM"
                 call abort()
               endif
               hgt_1d(nsum) = height
               IF(height.LT.-990.) height = 0.0
               xl1 = xl1 + height * float(zslm(ii,jj))
               xs1 = xs1 + height * float(1-zslm(ii,jj))
               xw1 = xw1 + height
               xw2 = xw2 + height ** 2
            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

               !IF(SLM(I,J).NE.0.) THEN
               IF(slm(i,j) .NE. 0. .OR. land_frac(i,j) > 0.) THEN
                  IF (xland > 0) THEN
                    oro(i,j)= xl1 / xland
                  ELSE
                    oro(i,j)= xs1 / xwatr
                  ENDIF
               ELSE
                  IF (xwatr > 0) THEN
                    oro(i,j)= xs1 / xwatr
                  ELSE
                    oro(i,j)= xl1 / xland
                  ENDIF
               ENDIF

               var(i,j)=sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
               do i1 = 1, nsum
                  xw4 = xw4 + (hgt_1d(i1) - oro(i,j)) ** 4
               enddo   

               IF(var(i,j).GT.1.) THEN
                  var4(i,j) = min(xw4/xnsum/var(i,j) **4,10.)
               ENDIF

         ELSEIF(xnsum_all.GT.1.) THEN

               !IF(SLM(I,J).NE.0.) THEN
               IF(slm(i,j) .NE. 0. .OR. land_frac(i,j) > 0.) THEN
                  IF (xland_all > 0) THEN
                    oro(i,j)= xl1_all / xland_all
                  ELSE
                    oro(i,j)= xs1_all / xwatr_all
                  ENDIF
               ELSE
                  IF (xwatr_all > 0) THEN
                    oro(i,j)= xs1_all / xwatr_all
                  ELSE
                    oro(i,j)= xl1_all / xland_all
                  ENDIF
               ENDIF

               var(i,j)=sqrt(max(xw2_all/xnsum_all-
     &                         (xw1_all/xnsum_all)**2,0.))
               do i1 = 1, nsum_all
                  xw4_all = xw4_all + 
     &                     (hgt_1d_all(i1) - oro(i,j)) ** 4
               enddo   

               IF(var(i,j).GT.1.) THEN
                  var4(i,j) = min(xw4_all/xnsum_all/var(i,j) **4,10.)
               ENDIF
         ELSE
               print*, "no source points in MAKEMT2"
               call abort()
         ENDIF
   
! set orog to 0 meters at ocean.
!        IF (LAKE_FRAC(I,J) .EQ. 0. .AND. SLM(I,J) .EQ. 0.)THEN
         IF (lake_frac(i,j) .EQ. 0. .AND. land_frac(i,j) .EQ. 0.)THEN 
         oro(i,j) = 0.0
         ENDIF

       ENDDO
      ENDDO
!$omp end parallel do
      WRITE(6,*) "! MAKEMT2 ORO SLM VAR VAR4 DONE"
C
      deallocate(hgt_1d)
      deallocate(hgt_1d_all)
      RETURN
      END

      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

      SUBROUTINE makepc2(ZAVG,ZSLM,THETA,GAMMA,SIGMA,
     1           glat,im,jm,imn,jmn,lon_c,lat_c,slm)
C
C===  PC: principal coordinates of each Z avg orog box for L&M
C
      implicit none
      real, parameter :: rearth=6.3712e+6
      real, parameter :: d2r = 3.14159265358979/180. 
      integer :: im,jm,imn,jmn
      real  :: glat(jmn),deltax(jmn)
      INTEGER zavg(imn,jmn),zslm(imn,jmn)
      real lon_c(im+1,jm+1), lat_c(im+1,jm+1)
      real, intent(in) :: slm(im,jm)
      real hl(im,jm),hk(im,jm)
      real hx2(im,jm),hy2(im,jm),hxy(im,jm),hlprim(im,jm)
      real theta(im,jm),gamma(im,jm),sigma2(im,jm),sigma(im,jm)
      real pi,certh,delxn,deltay,xnsum,xland
      real xfp,yfp,xfpyfp,xfp2,yfp2
      real hi0,hip1,hj0,hjp1,hijax,hi1j1
      real lono(4),lato(4),loni,lati
      integer i,j,i1,j1,i2,jst,jen,numx,i0,ip1,ijax
      integer ilist(imn)
      logical inside_a_polygon
      LOGICAL 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
      delxn = 360./imn      ! MOUNTAIN DATA RESOLUTION
      deltay =  certh / float(jmn)
      print *, 'MAKEPC2: DELTAY=',deltay
C
      DO j=1,jmn
         glat(j) = -90. + (j-1) * delxn + delxn * 0.5
         deltax(j) = deltay * cos(glat(j)*d2r)
      ENDDO
C
C---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
C

C... DERIVITIVE TENSOR OF HEIGHT
C
!$omp parallel do
!$omp* private (j,i,xnsum,xland,xfp,yfp,xfpyfp,
!$omp*          xfp2,yfp2,lono,lato,jst,jen,ilist,numx,j1,i2,i1,
!$omp*          loni,lati,i0,ip1,hi0,hip1,hj0,hjp1,ijax,
!$omp*          hijax,hi1j1)
      jloop : DO j=1,jm
!        print*, "J=", J
        iloop : DO i=1,im
          hx2(i,j) = 0.0
          hy2(i,j) = 0.0
          hxy(i,j) = 0.0
          xnsum = 0.0
          xland = 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.

            lono(1) = lon_c(i,j) 
            lono(2) = lon_c(i+1,j) 
            lono(3) = lon_c(i+1,j+1) 
            lono(4) = lon_c(i,j+1) 
            lato(1) = lat_c(i,j) 
            lato(2) = lat_c(i+1,j) 
            lato(3) = lat_c(i+1,j+1) 
            lato(4) = lat_c(i,j+1) 
            call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)

            do j1 = jst, jen; do i2 = 1, numx
              i1 = ilist(i2)         
              loni = i1*delxn
              lati = -90 + j1*delxn
              inside : if(inside_a_polygon(loni*d2r,lati*d2r,4,
     &           lono*d2r,lato*d2r))then

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

                  xland = xland + float(zslm(i1,j1))
                  xnsum = xnsum + 1.
C
                  hi0 =  float(zavg(i0,j1))
                  hip1 =  float(zavg(ip1,j1))
C
                  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. 1  .and. j1 .ne. jmn ) 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. 1 ) 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. jmn ) 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
               ENDIF inside
C
C === average the HX2, HY2 and HXY
C === This will be done over all land
C
               ENDDO
            ENDDO
C
C ===  HTENSR 
C
         xnsum_gt_1 : IF(xnsum.GT.1.) THEN
               IF(slm(i,j).NE.0.) THEN
                  IF (xland > 0) THEN
                    hx2(i,j) =  xfp2  / xland
                    hy2(i,j) =  yfp2  / xland
                    hxy(i,j) =  xfpyfp / xland
                  ELSE
                    hx2(i,j) =  xfp2  / xnsum
                    hy2(i,j) =  yfp2  / xnsum
                    hxy(i,j) =  xfpyfp / xnsum
                  ENDIF
               ENDIF
C=== degub testing
      if (debug) then
          print *," I,J,i1,j1:", i,j,i1,j1,
     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)) / d2r
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 xnsum_gt_1
                  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 iloop
      ENDDO jloop
!$omp end parallel do
      WRITE(6,*) "! MAKE Principal Coord  DONE"
C
      RETURN
      END SUBROUTINE makepc2
      
      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

      function get_lon_angle(dx,lat, DEGRAD)
      implicit none
      real dx, lat, degrad
      
      real get_lon_angle
         real, parameter :: radius = 6371200

         get_lon_angle = 2*asin( sin(dx/radius*0.5)/cos(lat) )*degrad
         
      end function get_lon_angle

      function get_lat_angle(dy, DEGRAD)
      implicit none
      real dy, degrad
      
      real get_lat_angle
         real, parameter :: radius = 6371200

         get_lat_angle = dy/radius*degrad
         
      end function get_lat_angle
      
      SUBROUTINE makeoa2(ZAVG,zslm,VAR,GLAT,OA4,OL,IOA4,ELVMAX,
     1           oro,oro1,xnsum,xnsum1,xnsum2,xnsum3,xnsum4,
     2           im,jm,imn,jmn,lon_c,lat_c,lon_t,lat_t,dx,dy,
     3           is_south_pole,is_north_pole )
      implicit none
      real, parameter :: missing_value = -9999.
      real, parameter :: d2r = 3.14159265358979/180.
      real, PARAMETER :: r2d=180./3.14159265358979
      integer im,jm,imn,jmn
      real    glat(jmn)
      INTEGER zavg(imn,jmn),zslm(imn,jmn)
      real    oro(im,jm),oro1(im,jm),elvmax(im,jm),zmax(im,jm)
      real    oa4(im,jm,4)
      integer ioa4(im,jm,4)
      real    lon_c(im+1,jm+1), lat_c(im+1,jm+1)
      real    lon_t(im,jm), lat_t(im,jm)
      real    dx(im,jm), dy(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)
      integer i,j,ilist(imn),numx,i1,j1,ii1
      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
      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
C   
C---- GLOBAL XLAT AND XLON ( DEGREE )
C
      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
C
      DO j=1,jm
        DO i=1,im
          xnsum(i,j) = 0.0
          elvmax(i,j) = oro(i,j)
          zmax(i,j)   = 0.0
C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
C     IN A GRID BOX
          xnsum1(i,j) = 0.0
          xnsum2(i,j) = 0.0
          xnsum3(i,j) = 0.0
          xnsum4(i,j) = 0.0
          oro1(i,j) = oro(i,j)
          elvmax(i,j) = zmax(i,j) 
        ENDDO
      ENDDO

! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
! ---  to JM or to JM1
!$omp parallel do
!$omp* private (j,i,hc,lono,lato,jst,jen,ilist,numx,j1,ii1,i1,loni,
!$omp*          lati,height)
      DO j=1,jm
!        print*, "J=", J 
        DO i=1,im
          hc = 1116.2 - 0.878 * var(i,j) 
          lono(1) = lon_c(i,j) 
          lono(2) = lon_c(i+1,j) 
          lono(3) = lon_c(i+1,j+1) 
          lono(4) = lon_c(i,j+1) 
          lato(1) = lat_c(i,j) 
          lato(2) = lat_c(i+1,j) 
          lato(3) = lat_c(i+1,j+1) 
          lato(4) = lat_c(i,j+1) 
          call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
          do j1 = jst, jen; do ii1 = 1, numx          
            i1 = ilist(ii1)
            loni = i1*delxn
            lati = -90 + j1*delxn
            if(inside_a_polygon(loni*d2r,lati*d2r,4,
     &          lono*d2r,lato*d2r))then

              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
              ENDIF   
            endif
          ENDDO ; ENDDO
        ENDDO
      ENDDO
!$omp end parallel do      
C
! ---  this will make work1 array take on oro's values on return
! ---  this will make work1 array take on oro's values on return
      DO j=1,jm
        DO i=1,im

          oro1(i,j) = oro(i,j)
          elvmax(i,j) = zmax(i,j) 
        ENDDO
      ENDDO
      
      DO kwd = 1, 4
        DO j=1,jm
          DO i=1,im
            oa4(i,j,kwd) = 0.0
            ol(i,j,kwd) = 0.0
          ENDDO
        ENDDO
      ENDDO
                                !
! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
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
!$omp parallel do
!$omp* private (j,i,lon,lat,kwd,dlon,dlat,lon1,lon2,lat1,lat2,
!$omp*          xnsum11,xnsum12,xnsum21,xnsum22,xnpu,xnpd,
!$omp*          xnsum1_11,xnsum2_11,hc_11, xnsum1_12,xnsum2_12,
!$omp*          hc_12,xnsum1_21,xnsum2_21,hc_21, xnsum1_22,
!$omp*          xnsum2_22,hc_22)
      DO j=1,jm
!       print*, "j = ", j
        DO i=1,im
          lon = lon_t(i,j)
          lat = lat_t(i,j)
          !--- for around north pole, oa and ol are all 0
          
          if(is_north_pole(i,j)) then
             print*, "set oa1 = 0 and ol=0 at i,j=", i,j
             do kwd = 1, 4
                  oa4(i,j,kwd) = 0.
                  ol(i,j,kwd) = 0.
             enddo
          else if(is_south_pole(i,j)) then
             print*, "set oa1 = 0 and ol=1 at i,j=", i,j
             do kwd = 1, 4
                oa4(i,j,kwd) = 0.
                ol(i,j,kwd) = 1.
             enddo    
          else
             
          !--- for each point, find a lat-lon grid box with same dx and dy as the cubic grid box
          dlon = get_lon_angle(dx(i,j), lat*d2r, r2d  )
          dlat = get_lat_angle(dy(i,j), r2d)
          !--- adjust dlat if the points are close to pole.
          if( lat-dlat*0.5<-90.) then
             print*, "at i,j =", i,j, lat, dlat, lat-dlat*0.5
             print*, "ERROR: lat-dlat*0.5<-90."
             call errexit(4)
          endif
          if( lat+dlat*2 > 90.) then
             dlat_old = dlat
             dlat = (90-lat)*0.5
             print*, "at i,j=",i,j," adjust dlat from ",
     &              dlat_old, " to ", dlat
          endif   
          !--- lower left 
          lon1 = lon-dlon*1.5
          lon2 = lon-dlon*0.5
          lat1 = lat-dlat*0.5
          lat2 = lat+dlat*0.5

          if(lat1<-90 .or. lat2>90) then
             print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          xnsum11 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,zslm,delxn)          

          !--- upper left 
          lon1 = lon-dlon*1.5
          lon2 = lon-dlon*0.5
          lat1 = lat+dlat*0.5
          lat2 = lat+dlat*1.5
          if(lat1<-90 .or. lat2>90) then
             print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          xnsum12 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,zslm,delxn)          

          !--- lower right
          lon1 = lon-dlon*0.5
          lon2 = lon+dlon*0.5
          lat1 = lat-dlat*0.5
          lat2 = lat+dlat*0.5
          if(lat1<-90 .or. lat2>90) then
             print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          xnsum21 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,zslm,delxn)          

          !--- upper right 
          lon1 = lon-dlon*0.5
          lon2 = lon+dlon*0.5
          lat1 = lat+dlat*0.5
          lat2 = lat+dlat*1.5
          if(lat1<-90 .or. lat2>90) then
             print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          
          xnsum22 = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,zslm,delxn)          
          
           xnpu = xnsum11 + xnsum12
           xnpd = xnsum21 + xnsum22
           IF (xnpd .NE. xnpu) oa4(i,j,1) = 1. - xnpd / max(xnpu , 1.)

           xnpu = xnsum11 + xnsum21
           xnpd = xnsum12 + xnsum22
           IF (xnpd .NE. xnpu) oa4(i,j,2) = 1. - xnpd / max(xnpu , 1.)

           xnpu = xnsum11 + (xnsum12+xnsum21)*0.5
           xnpd = xnsum22 + (xnsum12+xnsum21)*0.5
           IF (xnpd .NE. xnpu) oa4(i,j,3) = 1. - xnpd / max(xnpu , 1.)

           xnpu = xnsum12 + (xnsum11+xnsum22)*0.5
           xnpd = xnsum21 + (xnsum11+xnsum22)*0.5
           IF (xnpd .NE. xnpu) oa4(i,j,4) = 1. - xnpd / max(xnpu , 1.)

           
          !--- calculate OL3 and OL4
          !--- lower left 
          lon1 = lon-dlon*1.5
          lon2 = lon-dlon*0.5
          lat1 = lat-dlat*0.5
          lat2 = lat+dlat*0.5
          if(lat1<-90 .or. lat2>90) then
             print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
          endif          
          call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_11, xnsum2_11, hc_11)          

          !--- upper left 
          lon1 = lon-dlon*1.5
          lon2 = lon-dlon*0.5
          lat1 = lat+dlat*0.5
          lat2 = lat+dlat*1.5
          if(lat1<-90 .or. lat2>90) then
             print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_12, xnsum2_12, hc_12)          

          !--- lower right
          lon1 = lon-dlon*0.5
          lon2 = lon+dlon*0.5
          lat1 = lat-dlat*0.5
          lat2 = lat+dlat*0.5
          if(lat1<-90 .or. lat2>90) then
             print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_21, xnsum2_21, hc_21)          

          !--- upper right 
          lon1 = lon-dlon*0.5
          lon2 = lon+dlon*0.5
          lat1 = lat+dlat*0.5
          lat2 = lat+dlat*1.5
          if(lat1<-90 .or. lat2>90) then
             print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
          endif          
          call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_22, xnsum2_22, hc_22)           
                  
          ol(i,j,3) = (xnsum1_22+xnsum1_11)/(xnsum2_22+xnsum2_11)
          ol(i,j,4) = (xnsum1_12+xnsum1_21)/(xnsum2_12+xnsum2_21)

          !--- calculate OL1 and OL2
          !--- lower left 
          lon1 = lon-dlon*2.0
          lon2 = lon-dlon
          lat1 = lat
          lat2 = lat+dlat
          if(lat1<-90 .or. lat2>90) then
             print*, "at upper left i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_11, xnsum2_11, hc_11)          

          !--- upper left 
          lon1 = lon-dlon*2.0
          lon2 = lon-dlon
          lat1 = lat+dlat
          lat2 = lat+dlat*2.0
          if(lat1<-90 .or. lat2>90) then
             print*, "at lower left i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          
          call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_12, xnsum2_12, hc_12)          

          !--- lower right
          lon1 = lon-dlon
          lon2 = lon
          lat1 = lat
          lat2 = lat+dlat
          if(lat1<-90 .or. lat2>90) then
             print*, "at upper right i=,j=", i, j, lat, dlat,lat1,lat2
          endif          
          call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_21, xnsum2_21, hc_21)          

          !--- upper right 
          lon1 = lon-dlon
          lon2 = lon
          lat1 = lat+dlat
          lat2 = lat+dlat*2.0
          if(lat1<-90 .or. lat2>90) then
             print*, "at lower right i=,j=", i, j, lat, dlat,lat1,lat2
          endif
          
          call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn,glat,
     &     zavg,delxn, xnsum1_22, xnsum2_22, hc_22)           
                  
          ol(i,j,1) = (xnsum1_11+xnsum1_21)/(xnsum2_11+xnsum2_21)
          ol(i,j,2) = (xnsum1_21+xnsum1_22)/(xnsum2_21+xnsum2_22)         
          ENDIF          
        ENDDO
      ENDDO
!$omp end parallel do
      DO kwd=1,4
        DO j=1,jm
          DO i=1,im
            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,im
         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,*) "! MAKEOA2 EXIT"
C
      RETURN

      END SUBROUTINE makeoa2
 
      function spherical_distance(theta1,phi1,theta2,phi2)

      real, intent(in) :: theta1, phi1, theta2, phi2
      real :: spherical_distance, dot

      if(theta1 == theta2 .and. phi1 == phi2) then
        spherical_distance = 0.0
        return
      endif
  
      dot = cos(phi1)*cos(phi2)*cos(theta1-theta2) + sin(phi1)*sin(phi2)
      if(dot > 1. ) dot = 1.
      if(dot < -1.) dot = -1.
      spherical_distance = acos(dot)

      return

      end function spherical_distance
      
      subroutine get_mismatch_index(im_in, jm_in, geolon_in,geolat_in,
     &           bitmap_in,num_out, lon_out,lat_out, iindx, jindx )
      integer, intent(in) :: im_in, jm_in, num_out
      real,    intent(in) :: geolon_in(im_in,jm_in)
      real,    intent(in) :: geolat_in(im_in,jm_in)
      logical*1, intent(in) :: bitmap_in(im_in,jm_in)   
      real,    intent(in) :: lon_out(num_out), lat_out(num_out)
      integer, intent(out):: iindx(num_out), jindx(num_out)
      real, parameter :: max_dist = 1.e+20
      integer, parameter :: numnbr = 20
      integer :: i_c,j_c,jstart,jend
      real    :: lon,lat
      
      print*, "im_in,jm_in = ", im_in, jm_in
      print*, "num_out = ", num_out
      print*, "max and min of lon_in is", minval(geolon_in),
     &                                    maxval(geolon_in)
      print*, "max and min of lat_in is", minval(geolat_in),
     &                                    maxval(geolat_in)   
      print*, "max and min of lon_out is", minval(lon_out),
     &                                     maxval(lon_out)
      print*, "max and min of lat_out is", minval(lat_out),
     &                                     maxval(lat_out)   
      print*, "count(bitmap_in)= ", count(bitmap_in), max_dist
      
      do n = 1, num_out
        !      print*, "n = ", n
        lon = lon_out(n)
        lat = lat_out(n)
        !--- find the j-index for the nearest point
        i_c = 0; j_c = 0
        do j = 1, jm_in-1
          if(lat .LE. geolat_in(1,j) .and.
     &       lat .GE. geolat_in(1,j+1)) then
            j_c = j
          endif
        enddo
        if(lat > geolat_in(1,1)) j_c = 1
        if(lat < geolat_in(1,jm_in)) j_c = jm_in
        !      print*, "lat =", lat, geolat_in(1,jm_in), geolat_in(1,jm_in-1)
        ! The input is Gaussian grid.
        jstart = max(j_c-numnbr,1)
        jend = min(j_c+numnbr,jm_in)
        dist = max_dist
        iindx(n) = 0
        jindx(n) = 0
        !      print*, "jstart, jend =", jstart, jend
        do j = jstart, jend; do i = 1,im_in   
          if(bitmap_in(i,j) ) then
            !            print*, "bitmap_in is true"
            d = spherical_distance(lon_out(n),lat_out(n),
     &                             geolon_in(i,j), geolat_in(i,j))
            if( dist > d ) then
              iindx(n) = i; jindx(n) = j
              dist = d
            endif
          endif
        enddo; enddo
        if(iindx(n) ==0) then
          print*, "lon,lat=", lon,lat
          print*, "jstart, jend=", jstart, jend, dist
          print*, "ERROR in get mismatch_index: not find nearest points"
          call errexit(4)
        endif
      enddo   

      end subroutine get_mismatch_index
      
      subroutine interpolate_mismatch(im_in, jm_in, data_in,
     &                                num_out, data_out, iindx, jindx)
      integer, intent(in) :: im_in, jm_in, num_out
      real,    intent(in) :: data_in(im_in,jm_in)
      real,    intent(out):: data_out(num_out)
      integer, intent(in) :: iindx(num_out), jindx(num_out)
      
      do n = 1, num_out
        data_out(n) = data_in(iindx(n),jindx(n))
      enddo   

      end subroutine interpolate_mismatch
      
      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           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.
#ifdef IP_V4
      use ipolates_mod
#endif

      implicit none
      real, parameter :: missing_value = -9999.
      real, parameter :: d2r = 3.14159265358979/180.
      real, PARAMETER :: r2d=180./3.14159265358979
      integer im,jm,imn,jmn,imi,jmi
      real    glat(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)
      integer ioa4(im,jm,4)
      real    oa_in(imi,jmi,4), ol_in(imi,jmi,4)
      real    slm_in(imi,jmi)
      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
      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    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
      integer ibo(1), ibi(1)
      logical*1, allocatable :: bitmap_input(:,:)
      logical*1, allocatable :: bitmap_output(:,:)
      integer, allocatable :: ijsav_land_output(:)
      real,    allocatable :: lats_land_output(:)
      real,    allocatable :: lons_land_output(:)
      real,    allocatable :: output_data_land(:,:)
      real,    allocatable :: lons_mismatch_output(:)
      real,    allocatable :: lats_mismatch_output(:)
      real,    allocatable :: data_mismatch_output(:)
      integer, allocatable :: iindx(:), jindx(:)
C   
C---- GLOBAL XLAT AND XLON ( DEGREE )
C
      delxn = 360./imn      ! MOUNTAIN DATA RESOLUTION
C
      ijmdl_output = im*jm
      
      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
C
      DO j=1,jm
        DO i=1,im
          xnsum(i,j) = 0.0
          elvmax(i,j) = oro(i,j)
          zmax(i,j)   = 0.0
C---- COUNT NUMBER OF MODE. HIGHER THAN THE HC, CRITICAL HEIGHT
C     IN A GRID BOX
          xnsum1(i,j) = 0.0
          xnsum2(i,j) = 0.0
          xnsum3(i,j) = 0.0
          xnsum4(i,j) = 0.0
          oro1(i,j) = oro(i,j)
          elvmax(i,j) = zmax(i,j) 
        ENDDO
      ENDDO

! --- # of peaks > ZAVG value and ZMAX(IM,JM) -- ORO is already avg.
! ---  to JM or to JM1
      DO j=1,jm
!        print*, "J=", J 
        DO i=1,im
          hc = 1116.2 - 0.878 * var(i,j) 
          lono(1) = lon_c(i,j) 
          lono(2) = lon_c(i+1,j) 
          lono(3) = lon_c(i+1,j+1) 
          lono(4) = lon_c(i,j+1) 
          lato(1) = lat_c(i,j) 
          lato(2) = lat_c(i+1,j) 
          lato(3) = lat_c(i+1,j+1) 
          lato(4) = lat_c(i,j+1) 
          call get_index(imn,jmn,4,lono,lato,delxn,jst,jen,ilist,numx)
          do j1 = jst, jen; do ii1 = 1, numx          
            i1 = ilist(ii1)
            loni = i1*delxn
            lati = -90 + j1*delxn
            if(inside_a_polygon(loni*d2r,lati*d2r,4,
     &          lono*d2r,lato*d2r))then

              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
              ENDIF   
            endif
          ENDDO ; ENDDO
        ENDDO
      ENDDO
      
C
! ---  this will make work1 array take on oro's values on return
! ---  this will make work1 array take on oro's values on return
      DO j=1,jm
        DO i=1,im

          oro1(i,j) = oro(i,j)
          elvmax(i,j) = zmax(i,j) 
        ENDDO
      ENDDO
      
      DO kwd = 1, 4
        DO j=1,jm
          DO i=1,im
            oa4(i,j,kwd) = 0.0
            ol(i,j,kwd) = 0.0
          ENDDO
        ENDDO
      ENDDO

      !--- use the nearest point to do remapping.
      int_opt = 2
      ipopt=0
      kgds_input = 0
      kgds_input(1) = 4          ! OCT 6 - TYPE OF GRID (GAUSSIAN)
      kgds_input(2) = imi        ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
      kgds_input(3) = jmi        ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
      kgds_input(4) = 90000      ! OCT 11-13 - LAT OF ORIGIN
      kgds_input(5) = 0          ! OCT 14-16 - LON OF ORIGIN
      kgds_input(6) = 128        ! OCT 17 - RESOLUTION FLAG
      kgds_input(7) = -90000     ! OCT 18-20 - LAT OF EXTREME POINT
      kgds_input(8) = nint(-360000./imi)  ! OCT 21-23 - LON OF EXTREME POINT
      kgds_input(9)  = nint((360.0 / float(imi))*1000.0)
                                 ! OCT 24-25 - LONGITUDE DIRECTION INCR.
      kgds_input(10) = jmi /2    ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
      kgds_input(12) = 255       ! OCT 29 - RESERVED
      kgds_input(20) = 255       ! OCT 5  - NOT USED, SET TO 255


      kgds_output = -1
!      KGDS_OUTPUT(1) = IDRT       ! OCT 6 - TYPE OF GRID (GAUSSIAN)
      kgds_output(2) = im        ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
      kgds_output(3) = jm        ! OCT 9-10 - # PTS ON LONGITUDE CIRCLE
      kgds_output(4) = 90000      ! OCT 11-13 - LAT OF ORIGIN
      kgds_output(5) = 0          ! OCT 14-16 - LON OF ORIGIN
      kgds_output(6) = 128        ! OCT 17 - RESOLUTION FLAG
      kgds_output(7) = -90000     ! OCT 18-20 - LAT OF EXTREME POINT
      kgds_output(8) = nint(-360000./im)  ! OCT 21-23 - LON OF EXTREME POINT
      kgds_output(9)  = nint((360.0 / float(im))*1000.0)
                                  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
      kgds_output(10) = jm /2    ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
      kgds_output(12) = 255       ! OCT 29 - RESERVED
      kgds_output(20) = 255       ! OCT 5  - NOT USED, SET TO 255

      count_land_output=0    
      print*, "sum(slm) = ", sum(slm)
      do ij = 1, ijmdl_output
        j = (ij-1)/im + 1
        i = mod(ij-1,im) + 1         
        if (slm(i,j) > 0.0) then
          count_land_output=count_land_output+1
        endif
      enddo
      allocate(bitmap_input(imi,jmi))
      bitmap_input=.false.
      print*, "number of land input=", sum(slm_in)
      where(slm_in > 0.0) bitmap_input=.true.   
      print*, "count(bitmap_input)", count(bitmap_input)
      
      allocate(bitmap_output(count_land_output,1))
      allocate(output_data_land(count_land_output,1))
      allocate(ijsav_land_output(count_land_output))
      allocate(lats_land_output(count_land_output))
      allocate(lons_land_output(count_land_output))

      
      
      count_land_output=0
      do ij = 1, ijmdl_output
        j = (ij-1)/im + 1
        i = mod(ij-1,im) + 1 
        if (slm(i,j) > 0.0) then
          count_land_output=count_land_output+1
          ijsav_land_output(count_land_output)=ij
          lats_land_output(count_land_output)=lat_t(i,j)
          lons_land_output(count_land_output)=lon_t(i,j)
        endif
      enddo

      oa4 = 0.0
      ol = 0.0
      ibi = 1
      
      do kwd=1,4
        bitmap_output = .false.
         output_data_land = 0.0
        call ipolates(int_opt, ipopt, kgds_input, kgds_output,   
     &         (imi*jmi), count_land_output,               
     &          1, ibi, bitmap_input, oa_in(:,:,kwd),  
     &          count_land_output, lats_land_output,
     &          lons_land_output, ibo,  
     &          bitmap_output, output_data_land, iret)
        if (iret /= 0) then
          print*,'- ERROR IN IPOLATES ',iret
          call errexit(4)
        endif

        num_mismatch_land = 0     
        do ij = 1, count_land_output
          if (bitmap_output(ij,1)) then
            j = (ijsav_land_output(ij)-1)/im + 1
            i = mod(ijsav_land_output(ij)-1,im) + 1
            oa4(i,j,kwd)=output_data_land(ij,1)
          else  ! default value
            num_mismatch_land =  num_mismatch_land + 1
          endif
        enddo  
        print*, "num_mismatch_land = ", num_mismatch_land
        
        if(.not. allocated(data_mismatch_output)) then
          allocate(lons_mismatch_output(num_mismatch_land))
          allocate(lats_mismatch_output(num_mismatch_land))      
          allocate(data_mismatch_output(num_mismatch_land))
          allocate(iindx(num_mismatch_land))
          allocate(jindx(num_mismatch_land))
          
          num = 0     
          do ij = 1, count_land_output
            if (.not. bitmap_output(ij,1)) then   
              num = num+1
              lons_mismatch_output(num) = lons_land_output(ij)
              lats_mismatch_output(num) = lats_land_output(ij)
            endif
          enddo

          ! For those land points that with bitmap_output=.false. use
          ! the nearest land points to interpolate.
          print*,"before get_mismatch_index", count(bitmap_input)
          call get_mismatch_index(imi,jmi,lon_in*d2r,
     &        lat_in*d2r,bitmap_input,num_mismatch_land, 
     &        lons_mismatch_output*d2r,lats_mismatch_output*d2r,
     &        iindx, jindx )
        endif                
 
        data_mismatch_output = 0
        call interpolate_mismatch(imi,jmi,oa_in(:,:,kwd),
     &        num_mismatch_land,data_mismatch_output,iindx,jindx)  
       
        num = 0
        do ij = 1, count_land_output
          if (.not. bitmap_output(ij,1)) then   
            num = num+1
            j = (ijsav_land_output(ij)-1)/im + 1
            i = mod(ijsav_land_output(ij)-1,im) + 1
            oa4(i,j,kwd) = data_mismatch_output(num)
            if(i==372 .and. j== 611) then
            print*, "ij=",ij, num, oa4(i,j,kwd),iindx(num),jindx(num)
            endif
          endif
        enddo     
       
       
      enddo

      !OL
      do kwd=1,4
        bitmap_output = .false.
        output_data_land = 0.0
        call ipolates(int_opt, ipopt, kgds_input, kgds_output,   
     &         (imi*jmi), count_land_output,               
     &          1, ibi, bitmap_input, ol_in(:,:,kwd),  
     &          count_land_output, lats_land_output,
     &          lons_land_output, ibo,  
     &          bitmap_output, output_data_land, iret)
        if (iret /= 0) then
          print*,'- ERROR IN IPOLATES ',iret
          call errexit(4)
        endif

        num_mismatch_land = 0     
        do ij = 1, count_land_output
          if (bitmap_output(ij,1)) then
            j = (ijsav_land_output(ij)-1)/im + 1
            i = mod(ijsav_land_output(ij)-1,im) + 1
            ol(i,j,kwd)=output_data_land(ij,1)
          else  ! default value
            num_mismatch_land =  num_mismatch_land + 1
          endif
        enddo  
        print*, "num_mismatch_land = ", num_mismatch_land
        
        data_mismatch_output = 0
        call interpolate_mismatch(imi,jmi,ol_in(:,:,kwd),
     &        num_mismatch_land,data_mismatch_output,iindx,jindx)  
       
        num = 0
        do ij = 1, count_land_output
          if (.not. bitmap_output(ij,1)) then   
            num = num+1
            j = (ijsav_land_output(ij)-1)/im + 1
            i = mod(ijsav_land_output(ij)-1,im) + 1
            ol(i,j,kwd) = data_mismatch_output(num)
            if(i==372 .and. j== 611) then
            print*, "ij=",ij, num, ol(i,j,kwd),iindx(num),jindx(num)
            endif
          endif
        enddo     
       
       
      enddo
     
      deallocate(lons_mismatch_output,lats_mismatch_output)
      deallocate(data_mismatch_output,iindx,jindx)
      deallocate(bitmap_output,output_data_land)
      deallocate(ijsav_land_output,lats_land_output)
      deallocate(lons_land_output)
       
      DO kwd=1,4
        DO j=1,jm
          DO i=1,im
            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,im
         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,*) "! MAKEOA3 EXIT"
C
      RETURN
      END SUBROUTINE makeoa3
      
      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)
      enddo
      do j=1,jm
         jr = jm - j + 1
         do i=1,im
            ir = i + imb2
            if (ir .gt. im) ir = ir - im
            f(ir,jr) = wrk(i,j)
         enddo
      enddo
!
      tem = 0.0
      do i=1,im
        tem= tem + f(i,1)
      enddo
      tem = tem / im
      do i=1,im
         f(i,1) = tem
      enddo
!
      RETURN
      END

      subroutine rg2gg(im,jm,numi,a)
        implicit none
        integer,intent(in):: im,jm,numi(jm)
        real,intent(inout):: a(im,jm)
        integer j,ir,ig
        real r,t(im)
        do j=1,jm
          r=real(numi(j))/real(im)
          do ig=1,im
            ir=mod(nint((ig-1)*r),numi(j))+1
            t(ig)=a(ir,j)
          enddo
          do ig=1,im
            a(ig,j)=t(ig)
          enddo
        enddo
      end subroutine

      subroutine gg2rg(im,jm,numi,a)
        implicit none
        integer,intent(in):: im,jm,numi(jm)
        real,intent(inout):: a(im,jm)
        integer j,ir,ig
        real r,t(im)
        do j=1,jm
          r=real(numi(j))/real(im)
          do ir=1,numi(j)
            ig=nint((ir-1)/r)+1
            t(ir)=a(ig,j)
          enddo
          do ir=1,numi(j)
            a(ir,j)=t(ir)
          enddo
        enddo
      end subroutine

      SUBROUTINE minmxj(IM,JM,A,title)
      implicit none

      real a(im,jm),rmin,rmax
      integer i,j,im,jm
      character*8 title

      rmin=1.e+10
      rmax=-rmin
csela....................................................
csela if(rmin.eq.1.e+10)return
csela....................................................
      DO j=1,jm
      DO i=1,im
        if(a(i,j).ge.rmax)rmax=a(i,j)
        if(a(i,j).le.rmin)rmin=a(i,j)
      ENDDO
      ENDDO
      write(6,150)rmin,rmax,title
150   format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
C
      RETURN
      END

      SUBROUTINE mnmxja(IM,JM,A,imax,jmax,title)
      implicit none

      real a(im,jm),rmin,rmax
      integer i,j,im,jm,imax,jmax
      character*8 title

      rmin=1.e+10
      rmax=-rmin
csela....................................................
csela if(rmin.eq.1.e+10)return
csela....................................................
      DO j=1,jm
      DO i=1,im
        if(a(i,j).ge.rmax)then
                         rmax=a(i,j)
                        imax=i
                        jmax=j
        endif
        if(a(i,j).le.rmin)rmin=a(i,j)
      ENDDO
      ENDDO
      write(6,150)rmin,rmax,title
150   format('rmin=',e13.4,2x,'rmax=',e13.4,2x,a8,' ')
C
      RETURN
      END

      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

      subroutine read_g(glob)
      implicit none

      include 'netcdf.inc'

      integer*2, intent(out) :: glob(360*120,180*120)

      integer :: ncid, error, id_var, fsize

      fsize=65536

      error=nf__open("./topography.gmted2010.30s.nc",
     &               nf_nowrite,fsize,ncid)
      call netcdf_err(error, 'Open file topography.gmted2010.30s.nc' )
      error=nf_inq_varid(ncid, 'topo', id_var)
      call netcdf_err(error, 'Inquire varid of topo')
      error=nf_get_var_int2(ncid, id_var, glob)
      call netcdf_err(error, 'Read topo')
      error = nf_close(ncid)

      print*,' '
      call maxmin(glob,360*120*180*120,'global0')

      return
      end subroutine read_g

      subroutine maxmin(ia,len,tile)
ccmr
      implicit none
ccmr
      integer*2 ia(len)
      character*7 tile
      integer iaamax, iaamin, len, j, m, ja, kount
      integer(8) sum2,std,mean,isum
      integer i_count_notset,kount_9
! --- missing is -9999
c
      isum = 0
      sum2 = 0
      kount = 0
      kount_9 = 0
      iaamax = -9999999
ccmr  iaamin = 1
      iaamin =  9999999
      i_count_notset=0
           do 10 m=1,len
      ja=ia(m)
ccmr  if ( ja .lt. 0 ) print *,' ja < 0:',ja
ccmr  if ( ja .eq. -9999 ) goto 10
      if ( ja .eq. -9999 ) then
           kount_9=kount_9+1
           goto 10
      endif
      if ( ja .eq. -12345 ) i_count_notset=i_count_notset+1
ccmr  if ( ja .eq. 0 ) goto 11
      iaamax = max0( iaamax, ja )
      iaamin = min0( iaamin, ja )
!     iaamax = max0( iaamax, ia(m,j) )
!     iaamin = min0( iaamin, ia(m,j) )
  11  continue
      kount = kount + 1
      isum = isum + ja
ccmr  sum2 = sum2 + ifix( float(ja) * float(ja) )
      sum2 = sum2 + ja*ja
  10  continue
!
      mean = isum/kount
      std = ifix(sqrt(float((sum2/(kount))-mean**2)))
      print*,tile,' max=',iaamax,' min=',iaamin,' sum=',isum,
     &       ' i_count_notset=',i_count_notset
      print*,tile,' mean=',mean,' std.dev=',std,
     &       ' ko9s=',kount,kount_9,kount+kount_9 
      return
      end

      SUBROUTINE minmaxj(IM,JM,A,title)
      implicit none

      real(kind=4) a(im,jm),rmin,rmax,undef
      integer i,j,im,jm,imax,jmax,imin,jmin,iundef
      character*8 title,chara
      data chara/'        '/
      chara=title
      rmin=1.e+10
      rmax=-rmin
      imax=0
      imin=0
      jmax=0
      jmin=0
      iundef=0
      undef=-9999.
csela....................................................
csela if(rmin.eq.1.e+10)return
csela....................................................
      DO j=1,jm
      DO i=1,im
        if(a(i,j).ge.rmax)then
                        rmax=a(i,j)
                       imax=i
                       jmax=j
        endif
         if(a(i,j).le.rmin)then
            if ( a(i,j) .eq. undef ) then
                iundef = iundef + 1
            else
                        rmin=a(i,j)
                       imin=i
                       jmin=j
            endif
        endif
      ENDDO
      ENDDO
      write(6,150)chara,rmin,imin,jmin,rmax,imax,jmax,iundef
150   format(1x,a8,2x,'rmin=',e13.4,2i6,2x,'rmax=',e13.4,3i6)
C
      RETURN
      END
      
      subroutine latlon2xyz(siz,lon, lat, x, y, z)
      implicit none
      integer, intent(in) :: siz
      real, intent(in) :: lon(siz), lat(siz)
      real, intent(out) :: x(siz), y(siz), z(siz)
      
      integer n

      do n = 1, siz
        x(n) = cos(lat(n))*cos(lon(n))
        y(n) = cos(lat(n))*sin(lon(n))
        z(n) = sin(lat(n))
      enddo
      end

      FUNCTION spherical_angle(v1, v2, v3)
        implicit none
        real, parameter :: epsln30 = 1.e-30
        real, parameter :: pi=3.1415926535897931
        real v1(3), v2(3), v3(3)
        real  spherical_angle
 
        real px, py, pz, qx, qy, qz, ddd;
  
        ! vector product between v1 and v2 
        px = v1(2)*v2(3) - v1(3)*v2(2)
        py = v1(3)*v2(1) - v1(1)*v2(3)
        pz = v1(1)*v2(2) - v1(2)*v2(1)
        ! vector product between v1 and v3 
        qx = v1(2)*v3(3) - v1(3)*v3(2);
        qy = v1(3)*v3(1) - v1(1)*v3(3);
        qz = v1(1)*v3(2) - v1(2)*v3(1);

        ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
        if ( ddd <= 0.0 ) then
          spherical_angle = 0. 
        else 
          ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd);
          if( abs(ddd-1) < epsln30 ) ddd = 1;
          if( abs(ddd+1) < epsln30 ) ddd = -1;
          if ( ddd>1. .or. ddd<-1. ) then
            !FIX to correctly handle co-linear points (angle near pi or 0) */
            if (ddd < 0.) then
              spherical_angle = pi
            else
              spherical_angle = 0.
            endif
          else
            spherical_angle = acos( ddd )
          endif
        endif  

        return
      END  
      
      FUNCTION inside_a_polygon(lon1, lat1, npts, lon2, lat2)
        implicit none
        real, parameter :: epsln10 = 1.e-10
        real, parameter :: epsln8 = 1.e-8
        real, parameter :: pi=3.1415926535897931
        real, parameter :: range_check_criteria=0.05
        real :: anglesum, angle, spherical_angle
        integer i, ip1
        real lon1, lat1
        integer npts
        real lon2(npts), lat2(npts)
        real x2(npts), y2(npts), z2(npts)
        real lon1_1d(1), lat1_1d(1)
        real x1(1), y1(1), z1(1)
        real pnt0(3),pnt1(3),pnt2(3)
        logical inside_a_polygon
        real max_x2,min_x2,max_y2,min_y2,max_z2,min_z2
        !first convert to cartesian grid */
        call latlon2xyz(npts,lon2, lat2, x2, y2, z2);
        lon1_1d(1) = lon1
        lat1_1d(1) = lat1
        call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
        inside_a_polygon = .false.
        max_x2 = maxval(x2)
        if( x1(1) > max_x2+range_check_criteria ) return
        min_x2 = minval(x2)
        if( x1(1)+range_check_criteria < min_x2 ) return
        max_y2 = maxval(y2)
        if( y1(1) > max_y2+range_check_criteria ) return
        min_y2 = minval(y2)
        if( y1(1)+range_check_criteria < min_y2 ) return
        max_z2 = maxval(z2)
        if( z1(1) > max_z2+range_check_criteria ) return
        min_z2 = minval(z2)
        if( z1(1)+range_check_criteria < min_z2 ) return

        pnt0(1) = x1(1)
        pnt0(2) = y1(1)
        pnt0(3) = z1(1)
        
        anglesum = 0;
        do i = 1, npts
           if(abs(x1(1)-x2(i)) < epsln10 .and.
     &          abs(y1(1)-y2(i)) < epsln10 .and.
     &         abs(z1(1)-z2(i)) < epsln10 ) then ! same as the corner point
              inside_a_polygon = .true.
              return
           endif
           ip1 = i+1
           if(ip1>npts) ip1 = 1
           pnt1(1) = x2(i)
           pnt1(2) = y2(i)
           pnt1(3) = z2(i)
           pnt2(1) = x2(ip1)
           pnt2(2) = y2(ip1)
           pnt2(3) = z2(ip1)

           angle = spherical_angle(pnt0, pnt2, pnt1);
!           anglesum = anglesum + spherical_angle(pnt0, pnt2, pnt1);
           anglesum = anglesum + angle
        enddo

        if(abs(anglesum-2*pi) < epsln8) then
           inside_a_polygon = .true.
        else
           inside_a_polygon = .false.
        endif

        return
        
      end

      function get_xnsum(lon1,lat1,lon2,lat2,IMN,JMN,
     &                   glat,zavg,zslm,delxn)
        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)
        integer, intent(in) ::  zavg(imn,jmn),zslm(imn,jmn)
        integer i, j, ist, ien, jst, jen, i1
        real    oro, height
        real    xland,xwatr,xl1,xs1,slm,xnsum
        !---figure out ist,ien,jst,jen
        do j = 1, jmn
           if( glat(j) .GT. lat1 ) then
              jst = j
              exit
           endif
        enddo
        do j = 1, jmn
           if( glat(j) .GT. lat2 ) then
              jen = j
              exit
           endif
        enddo

        
        ist = lon1/delxn + 1
        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
          xnsum = 0
          xland = 0
          xwatr = 0
          xl1 = 0
          xs1 = 0
          do j = jst,jen
             do i1 = 1, ien - ist + 1
                i = ist + i1 -1
                if( i .LE. 0) i = i + imn
                if( i .GT. imn) i = i - imn
                xland = xland + float(zslm(i,j))
                xwatr = xwatr + float(1-zslm(i,j))
                xnsum = xnsum + 1.
                height = float(zavg(i,j)) 
                IF(height.LT.-990.) height = 0.0
                xl1 = xl1 + height * float(zslm(i,j))
                xs1 = xs1 + height * float(1-zslm(i,j))
             enddo
          enddo
          if( xnsum > 1.) THEN
             slm = float(nint(xland/xnsum))
               IF(slm.NE.0.) THEN
                  oro= xl1 / xland
               ELSE
                  oro = xs1 / xwatr
               ENDIF
          ENDIF
          
         get_xnsum = 0
         do j = jst, jen
            do i1= 1, ien-ist+1
               i = ist + i1 -1
               if( i .LE. 0) i = i + imn
               if( i .GT. imn) i = i - imn
               height = float(zavg(i,j))
               IF(height.LT.-990.) height = 0.0
               IF ( height .gt. oro ) get_xnsum = get_xnsum + 1
            enddo       
         enddo
!         if(verbose) print*, "get_xnsum=", get_xnsum, oro
         
      end function get_xnsum  
      
      subroutine get_xnsum2(lon1,lat1,lon2,lat2,IMN,JMN,
     &                   glat,zavg,delxn,xnsum1,xnsum2,hc)
        implicit none

        real, intent(out) :: xnsum1,xnsum2,hc
        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,slm,xnsum
        !---figure out ist,ien,jst,jen
        do j = 1, jmn
           if( glat(j) .GT. lat1 ) then
              jst = j
              exit
           endif
        enddo
        do j = 1, jmn
           if( glat(j) .GT. lat2 ) then
              jen = j
              exit
           endif
        enddo

        
        ist = lon1/delxn + 1
        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
          xw1 = 0
          xw2 = 0
          do j = jst,jen
             do i1 = 1, ien - ist + 1
                i = ist + i1 -1
                if( i .LE. 0) i = i + imn
                if( i .GT. imn) i = i - imn
                xnsum = xnsum + 1.
                height = float(zavg(i,j)) 
                IF(height.LT.-990.) height = 0.0
                xw1 = xw1 + height
                xw2 = xw2 + height ** 2
             enddo
          enddo
          var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
          hc = 1116.2 - 0.878 * var
         xnsum1 = 0
         xnsum2 = 0
         do j = jst, jen
            do i1= 1, ien-ist+1
               i = ist + i1 -1
               if( i .LE. 0) i = i + imn
               if( i .GT. imn) i = i - imn
               height = float(zavg(i,j))
               IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
                xnsum2 = xnsum2 + 1
            enddo       
         enddo
         
      end subroutine get_xnsum2 

      subroutine get_xnsum3(lon1,lat1,lon2,lat2,IMN,JMN,
     &                   glat,zavg,delxn,xnsum1,xnsum2,hc)
        implicit none

        real, intent(out) :: xnsum1,xnsum2
        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
        jen = jmn
        do j = 1, jmn
           if( glat(j) .GT. lat1 ) then
              jst = j
              exit
           endif
        enddo
        do j = 1, jmn
           if( glat(j) .GT. lat2 ) then
              jen = j
              exit
           endif
        enddo

        
        ist = lon1/delxn + 1
        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
         do j = jst, jen
            do i1= 1, ien-ist+1
               i = ist + i1 -1
               if( i .LE. 0) i = i + imn
               if( i .GT. imn) i = i - imn
               height = float(zavg(i,j))
               IF ( height .gt. hc ) xnsum1 = xnsum1 + 1
                xnsum2 = xnsum2 + 1
            enddo       
         enddo
         
      end subroutine get_xnsum3
      
      subroutine nanc(a,l,c)
      integer inan1,inan2,inan3,inan4,inaq1,inaq2,inaq3,inaq4
       real word
       integer itest
       equivalence(itest,word)
c
c signaling NaN
      data inan1/x'7F800001'/
      data inan2/x'7FBFFFFF'/
      data inan3/x'FF800001'/
      data inan4/x'FFBFFFFF'/
c
c  quiet NaN
c
      data inaq1/x'7FC00000'/
      data inaq2/x'7FFFFFFF'/
      data inaq3/x'FFC00000'/
      data inaq4/x'FFFFFFFF'/
c
      real(kind=8)a(l),rtc,t1,t2
      character*(*) c
c     t1=rtc()
cgwv        print *, ' nanc call ',c
      do k=1,l
      word=a(k)
       if( (itest  .GE. inan1 .AND. itest .LE. inan2) .OR.
     *      (itest  .GE. inan3 .AND. itest .LE. inan4) ) then
       print *,' NaNs detected at  word',k,' ',c
       return
         endif
        if( (itest  .GE. inaq1 .AND. itest .LE. inaq2) .OR.
     *      (itest  .GE. inaq3 .AND. itest .LE. inaq4) ) then
       print *,' NaNq detected at  word',k,' ',c
         return
           endif

 101  format(e20.10)
      end do
c     t2=rtc()
cgwv      print 102,l,t2-t1,c
 102  format(' time to check ',i9,' words is ',f10.4,' ',a24)
      return
       end

      real function timef()
      character(8) :: date
      character(10) :: time
      character(5) :: zone
      integer,dimension(8) :: values
      integer :: total
      real :: elapsed
      call date_and_time(date,time,zone,values)
      total=(3600*values(5))+(60*values(6))
     *      +values(7)
      elapsed=float(total) + (1.0e-3*float(values(8)))
      timef=elapsed
      return
      end