text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis

	Program gridgen_model
C**********************************************************************
c Portions of the following code were taken from the RAMS software and
c   were used by permission of the *ASTER Division/Mission Research
c   Corporation.
C**********************************************************************
C*	This program will be used to map model grids based on the RAMS*
C*  version 2b which uses polar stereographic projections. Other      *
C*  projections have since been added.                                *
c*                                                                    *
C*********************************************************************
c Beginning July 2000, this software was updated to accomodate the    *
c WRFSI (Weather Research and Forecasting Standard Initialization)    *
c grid requirements.  This included adding stagger calculations, and  *
c map projection factors (map factor and components of coriolis parm, *
c landuse, soil type, vegetation greenness fraction, mean annual soil *
c temperature, albedo, and sea surface temperature.  Significant work *
c has occurred to reduce memory requirements.  Many sections of code  *
c previously in this file now reside in file gridgen_utils.f          *
c For further information contact                                     *
c Mr. John Smart NOAA/FSL, smart@fsl.noaa.gov (303-497-6590)          *
c*********************************************************************

        use horiz_interp

        integer istag
        integer n_staggers
        integer lf
        character c6_maproj*6
        character c8_maproj*8
        character c10_grid_f*10 ! Type of domain (nest7grid, wrfsi)/
        character c_dataroot*200
        character cstaticdir*200
        logical   localize
        real,     allocatable  ::  mother_lat(:,:)
        real,     allocatable  ::  mother_lon(:,:)
        real,     allocatable  ::  dum2d(:,:)
        real      La1,Lo1,La2,Lo2
        real      LoV,Latin1,Latin2
        real      Dx,Dy

        integer lli_orig
        integer llj_orig
        integer uri_orig
        integer urj_orig

        include 'grid_fname.cmn'

        write(6,*)
        write(6,*)' gridgen_model: start'

        call get_config(istatus)
        if(istatus.ne.1)then
           print*,'Error returned: get_config'
           print*,'Abort'
           stop
        endif

        call find_domain_name(c_dataroot,c10_grid_f,istatus)
        if(istatus .ne. 1)then
            write(6,*) 'Error: find domain name - stop'
            stop
        endif
        call s_len(c10_grid_f,lf)

        if(c10_grid_f(1:lf).eq.'wrfsi')then
           call get_num_domains(num_domains,istatus)
           if(istatus.ne.1)then
              print *,'Error returned: get_num_domains'
              print *,'Abort'
              stop
           endif
        else
           num_domains=1
        endif
c -----------------------------------------------------------------
c variable "nest" is common in src/include/grid_fname.cmn
c and is used to communicate the domain # throughout the code.
c -----------------------------------------------------------------
        do nest=1,num_domains
c
c determine if this domain requires localizing or re-localizing
c
           call get_directory('static',cstaticdir,len_dir)

           call rd_static_attr(cstaticdir,nest
     .,c10_grid_f,nx_dom,ny_dom,Dx,Dy,La1,Lo1,Latin1
     .,Latin2,LoV,c8_maproj,istatus)

c -------------------------------------------------------------------------
c routine below accesses environment variable $FORCE_LOCALIZATION which, if
c set, will return variable "localize" = .true.
c -------------------------------------------------------------------------

           call eval_localization(cstaticdir,nest,localize
     .,c10_grid_f,La1,Lo1,istatus)

           call get_parent_id(iparent_id,istatus)

           if(localize)then

c get_mother_dims,  uses parent_id as index to determine mother domain specs
           if(nest.gt.1)then
              call get_mother_dims(nx_dom,ny_dom,istatus)
              allocate (mother_lat(nx_dom,ny_dom),
     +                  mother_lon(nx_dom,ny_dom),
     +                  dum2d(nx_dom,ny_dom),stat=istat)
              if(istat.ne.0)then
                 print*,'cannot deallocate mother arrays: ',istat
                 stop
              endif

              call get_mother_domain(nx_dom,ny_dom,mother_lat
     +,mother_lon,istatus)
              if(istatus.lt.1)then
                 print*
                 print *,'Error reading MOTHER lat, lon data.'
                 print *,'Abort'
                 print*
                 stop
              endif
           endif

           call get_grid_dim_xy(nx_dom,ny_dom,istatus)
           call get_ratio_to_parent(iratio_to_parent,istatus)
           call get_domain_origin(lli_orig,llj_orig,uri_orig,urj_orig
     +,istatus)
c
c use nx_dom and ny_dom generically regardless of a nest or not.
c
           call get_grid_spacing(Dx,istatus)
           call get_c6_maproj(c6_maproj,istatus)
           call downcase(c6_maproj,c6_maproj)

           if(c6_maproj.ne.'rotlat')Dy=Dx

           print*
           print*,'---------------------------------------------'
           print*,'Domain #         = ',nest
           print*,'Parent of nest   = ',iparent_id
           print*,'Ratio to parent  = ',iratio_to_parent
           print*,'Grid_spacing(m)  = ',Dx
           print*,'Grid dims (nx/ny)= ',nx_dom,ny_dom
           print*,'Lower left  i/j  = ',lli_orig,llj_orig
           print*,'Upper right i/j  = ',uri_orig,urj_orig

           if(nest.gt.1)then
              print*,'Calculate center lat/lon '
              La1=mother_lat(lli_orig,llj_orig)
              Lo1=mother_lon(lli_orig,llj_orig)
              La2=mother_lat(uri_orig,urj_orig)
              Lo2=mother_lon(uri_orig,urj_orig)
              print*,'*** SW lat/lon = ',La1,Lo1,'***'
              print*,'*** NE lat/lon = ',La2,Lo2,'***'
              call get_earth_radius(erad,istatus)
c Get X/Y for nested grid center
              call latlon_to_xy(La1,Lo1,ERAD,Xsw,Ysw)
              call latlon_to_xy(La2,Lo2,ERAD,Xne,Yne)
              Xcen=(Xne+Xsw)/2.
              Ycen=(Yne+Ysw)/2.
              call xy_to_latlon(Xcen,Ycen,erad,rlatcen,rloncen)
              print*,'*** Center lat/lon= ',rlatcen,rloncen,' ***'
           endif

           call get_n_staggers(n_staggers,istatus)
           call get_stagger_index(istag,istatus)

           if(istatus.ne.1)then
              print*,'Error getting namelist info for',
     +': ',c10_grid_f(1:lf)
              stop
           endif

           if(allocated(mother_lat))then
              deallocate (mother_lat, mother_lon, dum2d)
           endif

           print*
c
c -------------------------------------------------------------------
c
           call Gridmap_sub(nx_dom,ny_dom,n_staggers,istag,nest
     +,rlatcen,rloncen,lli_orig,llj_orig,uri_orig,urj_orig
     +,iparent_id,iratio_to_parent,istatus)
c
c -------------------------------------------------------------------
c
           print*,' Gridmap_sub finished '
           print*
           if(istatus.eq.1)then

              print*,' Gridmap_sub finish successfully'

              if( allocated(mother_lat) )then
               deallocate (mother_lat,mother_lon,dum2d,stat=istat)
               if(istat.ne.0)then
                  print*
                  print*,'Cannot deallocate mother arrays: ',istat
                  print*,'Abort'
                  print*
                  stop
               endif
              endif

           else
              print*,'Problem in Gridmap_sub. Terminating'
              stop
           endif

           else  !localize or not!

           print*
           print*,' No need to re-localize this domain '
           print*,' Localization namelist parms havent changed'
           print*

           if( allocated(mother_lat) )then
              deallocate (mother_lat,mother_lon,dum2d,stat=istat)
              if(istat.ne.0)then
                 print*
                 print*,'Cannot deallocate mother arrays: ',istat
                 print*,'Abort'
                 print*
                 stop
              endif
           endif

           endif !localize case

        enddo

 999	print*,' gridgen_model finish: istatus = ',istatus
        print*

        stop
        end
c
c --------------------------------------------------------------------
c
        subroutine Gridmap_sub(nnxp,nnyp,n_staggers,istag,nest
     +,mdlat,mdlon,lli,llj,iuri,iurj,iparent_id,iratio_to_parent
     +,istatus)

        include 'trigd.inc'

        use mem_namelist, ONLY: laps_cycle_time, l_fsf_gridgen

        logical exist,new_DEM
        logical lforce_ter_2_zero
!mp
	logical categorical, useland
!mp
        logical l_topo_wps, l_parse, l_topo_1s

        integer nnxp,nnyp,mode
        integer ngrids
        integer n_staggers
        integer nest
        integer iter
        integer istag

        Real  mdlat,mdlon
        Real  xmn(nnxp),ymn(nnyp)
        Real  xtn(nnxp,n_staggers)
        Real  ytn(nnyp,n_staggers)

        real, allocatable ::  topt_10   (:,:)
        real, allocatable ::  topt_10_s (:,:)
        real, allocatable ::  topt_10_ln(:,:)
        real, allocatable ::  topt_10_lt(:,:)
        real, allocatable ::  topt_30   (:,:)
        real, allocatable ::  topt_30_s (:,:)
        real, allocatable ::  topt_30_ln(:,:)
        real, allocatable ::  topt_30_lt(:,:)

        real  topt_out   (nnxp,nnyp)
        real  topt_out_s (nnxp,nnyp)
        real  topt_out_ln(nnxp,nnyp)
        real  topt_out_lt(nnxp,nnyp)

        integer maxdatacat
        parameter (maxdatacat=24)

        real, allocatable ::  adum2d  (:,:)
        real, allocatable ::  adum3d  (:,:,:)

        real, allocatable ::  GEODAT2D(:,:)
        real, allocatable ::  GEODAT3D(:,:,:)

        real lats(nnxp,nnyp,n_staggers)
        real lons(nnxp,nnyp,n_staggers)

	real*8 r8lat(nnxp,nnyp), r8lon(nnxp,nnyp)

	real, allocatable, dimension(:,:):: hlat,hlon,vlat,vlon

        character (len=3),   allocatable :: var(:)
        character (len=125), allocatable :: comment(:)

        character*2   cnest
        character*30  ctopo,clatlon,corners,ctopog
        character*30  ctopos,clatlons,cornerss,ctopogs
        character*30  clatlon2d,clatlons2d
        character*131 model

        character*200 path_to_topt30s
        character*200 path_to_topt10m
        character*200 path_to_pctl10m
        character*200 path_to_soiltype_top_30s
        character*200 path_to_soiltype_bot_30s
        character*200 path_to_green_frac    !no reference to res since there is both
                                            !10m and 8.64m
        character*200 path_to_soiltemp_1deg
        character*200 path_to_luse_30s
        character*200 path_to_albedo        !albedo = 0.144 deg res or 8.64m
        character*200 path_to_maxsnoalb
        character*200 path_to_islope

        character*255 filename
        character*255 filename_wps
        character*200 c_dataroot
        character*200 cdl_dir
        character*180 static_dir
        character*10  c10_grid_f           ! Type of domain (nest7grid, wrfsi)
        character*10  c10_grid_fname       ! Actual filename, cdl (nest7grid for now)
        character*6   c6_maproj
        character*1   cdatatype
        integer len,lf,lfn,ns,avgelem,zinelem
        integer ishow_timer,init_timer
        integer itstatus

        real,  allocatable ::  rmapdata(:,:,:)
        real,  allocatable ::  data(:,:,:)     !primary output array storage unit.

        interface

          subroutine adjust_geog(nnxp,nnyp,ncat,ctype
     &,istat_dat,lat,topt_out,landmask,geog_data,istatus)

          integer nnxp,nnyp
          integer ncat
          integer istat_dat
          integer istatus

          character*(*) ctype
          real    lat(nnxp,nnyp)
          real    landmask(nnxp,nnyp)
          real    topt_out(nnxp,nnyp)
          real    geog_data(nnxp,nnyp,ncat)
          end subroutine

          subroutine proc_geodat(nx_dom,ny_dom,ncat
     1,path_to_tile_data,dom_lats_in,dom_lons_in,lmask_out
     1,geodat,istatus)

          integer nx_dom
          integer ny_dom
          integer ncat
          character*(*) path_to_tile_data
          real    dom_lats_in(nx_dom,ny_dom)
          real    dom_lons_in(nx_dom,ny_dom)
          real    lmask_out(nx_dom,ny_dom)
          real    geodat(nx_dom,ny_dom,ncat)
          integer istatus
          end subroutine

        end interface

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

        itstatus = init_timer()

        call find_domain_name(c_dataroot,c10_grid_f,istatus)
        if(istatus .ne. 1)then
            write(6,*) 'Error: returned fro find_domain_name '
            return
        endif

        call s_len(c10_grid_f,lf)

!mp  - need to know map projection information from the start
        call get_c6_maproj(c6_maproj,istatus)
        if(istatus .ne. 1)then
            write(6,*)' Error calling laps routine'
            return
        endif
        write(6,*)' c6_maproj = ',c6_maproj
!mp

c add 12 for albedo month 1 - 12.
c add 1 more for islope JS: 2-25-03
        if(c10_grid_f(1:lf).eq.'wrfsi')then
!mp
            IF (c6_maproj .eq. 'rotlat') THEN
               write(6,*) 'ngrids set to 103'
               ngrids=103 !!!
            ELSE !normal
               ngrids=112     !97+12   ! 2d grids (including %dist for landuse
                                       ! and two soiltype categories, green frac, albedo and others).
	    ENDIF

        else
           ngrids=38      !26+12   doesn't include terrain slope index and other stuff.
        endif

        allocate (data(nnxp,nnyp,ngrids),stat=istat)

        if(istat.ne.0)then
           print*,'Error allocating data object: array data ',istat
           print*,'Terminating: Maybe not enough memory'
           stop
        endif

        model = 'MODEL 4 delta x smoothed filter\0'

        icount_10 = 0
        icount_30 = 0
        icount_ramp = 0

        call get_directory('static',static_dir,lens)

c iplttopo is 1 if you want to plot the topography
c   the 30s topo covers the continental US
cc        itoptfn_30=static_dir(1:len)//'model/topo_30s/U'
c   the 10m topo covers the world
cc        itoptfn_10=static_dir(1:len)//'model/topo_10m/H'

        call get_path_to_topo_10m(path_to_topt10m,istatus)
        if(istatus .ne. 1)then
           write(6,*) 'Error getting path_to_topt10m'
           return
        endif

        call get_path_to_topo_30s(path_to_topt30s,istatus)
        if(istatus .ne. 1)then
           write(6,*) 'Error getting path_to_topt30s'
           return
        endif

!       l_topo_wps = l_parse(path_to_topt30s,'wps')
        l_topo_1s = l_parse(path_to_topt30s,'topo_1s')

        call get_path_to_soiltype_top(path_to_soiltype_top_30s
     +,istatus)
        if(istatus .ne. 1)then
           write(6,*)'Error getting path_to_soiltype_top_30s'
           return
        endif

        call get_path_to_soiltype_bot(path_to_soiltype_bot_30s
     +,istatus)
        if(istatus .ne. 1)then
           write(6,*)'Error getting path_to_soiltype_bot_30s'
           return
        endif

        call get_path_to_landuse_30s(path_to_luse_30s,istatus)
        if(istatus .ne. 1)then
           write(6,*)'Error getting path_to_luse_30s'
           return
        endif

        call get_path_to_green_frac(path_to_green_frac,istatus)
        if(istatus .ne. 1)then
           print*,'Error getting path_to_green_frac'
           return
        endif

        call get_path_to_soiltemp_1deg(path_to_soiltemp_1deg
     &,istatus)
        if(istatus .ne. 1)then
           print*, 'Error getting path_to_soiltemp_1deg'
           return
        endif

        call get_path_to_albedo(path_to_albedo,istatus)
        if(istatus .ne. 1)then
           print*, 'Error getting path_to_albedo'
           return
        endif

        call get_path_to_maxsnoalb(path_to_maxsnoalb,istatus)
        if(istatus .ne. 1)then
           print*, 'Error getting path_to_maxsnoalb'
           return
        endif

        call get_path_to_islope(path_to_islope,istatus)
        if(istatus .ne. 1)then
           print*, 'Error getting path_to_islope'
           return
        endif

        call s_len(path_to_topt10m,len)
        if(len.gt.0)then
           print*,'path to toptl0m:        ',path_to_topt10m(1:len)
        endif
        path_to_topt10m(len+1:len+2)='/H'

        call s_len(path_to_soiltype_top_30s,len)
        print*,'path to soiltype_top:   '
     .,path_to_soiltype_top_30s(1:len)
        path_to_soiltype_top_30s(len+1:len+2)='/O'
        print*,'path to soiltype_bot:   '
     .,path_to_soiltype_bot_30s(1:len)
        path_to_soiltype_bot_30s(len+1:len+2)='/O'

        call s_len(path_to_luse_30s,len)
        print*,'path to landuse 30s: ',path_to_luse_30s(1:len)
        path_to_luse_30s(len+1:len+2)='/V'

        call s_len(path_to_green_frac,len)
        print*,'path to green frac:     ',path_to_green_frac(1:len)
        path_to_green_frac(len+1:len+2)='/G'

        call s_len(path_to_soiltemp_1deg,len)
        print*,'path to soiltemp 1deg:  ',path_to_soiltemp_1deg(1:len)
        path_to_soiltemp_1deg(len+1:len+2)='/T'

        call s_len(path_to_albedo,len)
        print*,'path to albedo:         ',path_to_albedo(1:len)
        path_to_albedo(len+1:len+2)='/A'

        call s_len(path_to_maxsnoalb,len)
        print*,'path to max snow albedo: ',path_to_maxsnoalb(1:len)
        path_to_maxsnoalb(len+1:len+2)='/M'

        call s_len(path_to_islope,len)
        print*,'path to islope categorical: ',path_to_islope(1:len)
        path_to_islope(len+1:len+2)='/I'

        call get_topo_parms(silavwt_parm,toptwvl_parm,istatus)
	if (istatus .ne. 1) then
           write (6,*) 'Error getting terrain smoothing parms'
	   return
	endif

        call get_r_missing_data(r_missing_data,istatus)
	if (istatus .ne. 1) then
           write (6,*) 'Error getting r_missing_data'
	   return
	endif

!       Silhouette weighting parameter
        silavwt=silavwt_parm

!       Terrain wavelength for filtering
        toptwvl=toptwvl_parm

        iplttopo=1

C*********************************************************************
        call get_gridnl(mode)

c calculate delta x and delta y using grid and map projection parameters
        call get_standard_latitudes(std_lat,std_lat2,istatus)
        if(istatus .ne. 1)then
            write(6,*)' Error calling laps routine'
            return
        endif
        write(6,*)' standard_lats = ',std_lat,std_lat2

	IF (c6_maproj .ne. 'rotlat') THEN

          call get_grid_spacing(grid_spacing_m,istatus)
          if(istatus .ne. 1)then
            write(6,*)' Error calling laps routine'
            return
          endif

	ELSE

          write(6,*)
     1  'ERROR: get_nmm_grid_spacing not currently supported for rotlat'
          stop
!         call get_nmm_grid_spacing(dlmd,dphd,grid_spacing_m)

	ENDIF

        write(6,*)' grid_spacing = ',grid_spacing_m

        if(nest.eq.1)then
           print*,'get grid center from namelist'
           call get_grid_center(mdlat,mdlon,istatus)
           if(istatus .ne. 1)then
              write(6,*)' Error calling laps routine'
              return
           endif
        endif

        print*,' Gridmap_sub: grid_center = ',mdlat,mdlon

        call get_c6_maproj(c6_maproj,istatus)
        if(istatus .ne. 1)then
            write(6,*)' Error calling laps routine'
            return
        endif
        write(6,*)' c6_maproj = ',c6_maproj

        call get_standard_longitude(std_lon,istatus)
        if(istatus .ne. 1)then
            write(6,*)' Error calling laps routine'
            return
        endif
        write(6,*)' std_lon = ',std_lon

        call get_earth_radius(erad,istatus)
        if(istatus .ne. 1)then
            write(6,*)' Error calling laps routine'
            return
        endif
        write(6,*)' Earth Radius = ',erad

        if(c6_maproj .eq. 'plrstr')then
            call get_ps_parms(std_lat,std_lat2,grid_spacing_m,phi0
     1                       ,grid_spacing_proj_m)

            deltax = grid_spacing_proj_m

            if(phi0 .eq. 90.)then
                write(6,*)' Projection is tangent to earths surface'
            else
                write(6,*)' Projection is secant to earths surface'
            endif

            if(std_lat2 .eq. +90.)then
                write(6,*)' Note, grid spacing will equal '
     1                    ,grid_spacing_m,' at a latitude of ',std_lat
!               write(6,*)' deltax, deltay ',deltax,deltay
!    1                   ,' at the north pole'

            elseif(std_lat2 .eq. -90.)then
                write(6,*)' Note, grid spacing will equal '
     1                    ,grid_spacing_m,' at a latitude of ',std_lat
!               write(6,*)' deltax, deltay ',deltax,deltay
!    1                   ,' at the south pole'

            else ! abs(std_lat2) .ne. 90. (local stereographic)
                write(6,*)' The standard latitude ',std_lat,' is'
     1                   ,' relative to where the pole'
                write(6,*)' of the map projection is: lat/lon '
     1                   ,std_lat2,std_lon
!               write(6,*)' deltax, deltay ',deltax,deltay
!    1                   ,' at the projection pole'
                if(std_lat .ne. +90.)then
                    write(6,*)' Note: std_lat should usually be set'
     1                       ,' to +90. for local stereographic'
                endif

            endif

        else ! c6_maproj .ne. 'plrstr'
            deltax = grid_spacing_m

        endif ! c6_maproj .eq. 'plrstr'

        deltay = deltax
        write(6,*)' deltax, deltay (in projection plane) '
     1             ,deltax,deltay

c*********************************************************************
c in arrays lats/lons, the first stagger is actually not staggered
c but is the usual lat/lon values at non-staggered grid points. These
c are used for the analysis (LAPS-nest7grid). For WRFSI, we use the
c mass c-stagger (#4); staggered in both x and y.

!mp

	IF (c6_maproj .eq. 'rotlat') THEN
          write(6,*) 'dlmd, dphd, mdlat, mdlon: ',dlmd,dphd,mdlat,mdlon

          allocate(hlat(nnxp,nnyp),hlon(nnxp,nnyp))
          allocate(vlat(nnxp,nnyp),vlon(nnxp,nnyp))

          call etall(nnxp,nnyp,mdlat,mdlon,dlmd,dphd,hlat,hlon
     &                  ,vlat,vlon)

	  write(6,*) 'EGRID CORNERS'
          write(6,*) 'hlat(1,1),hlon(1,1): ', hlat(1,1),hlon(1,1)
          write(6,*) 'hlat(1,jm),hlon(1,jm): ',hlat(1,nnyp),hlon(1,nnyp)
          write(6,*) 'hlat(im,1),hlon(im,1): ', hlat(nnxp,1),
     &			hlon(nnxp,1)
          write(6,*) 'hlat(im,jm),hlon(im,jm): ', hlat(nnxp,nnyp),
     &			hlon(nnxp,nnyp)

          DO J=1,NNYP
            DO I=1,NNXP
                data(I,J,1)=hlat(I,J)
                data(I,J,2)=hlon(I,J)
                data(I,J,3)=vlat(I,J)
                data(I,J,4)=vlon(I,J)
            ENDDO
          ENDDO

	ELSEIF(c6_maproj .eq. 'icshdr')then ! Read Icosahedral grid from a file
          write(6,*)' Read Icosahedral grid from a file'

          filename=static_dir(1:lens)//'ltln5C.txt'
          open(7,file=filename)
          read(7,*) lats, lons
          close(7)

          rpd = 3.1415926535897932 / 180.

          DO J=1,NNYP
          DO I=1,NNXP
          DO IS = 1,n_staggers
              lats(i,j,is) = lats(i,j,is) / rpd
              lons(i,j,is) = lons(i,j,is) / rpd
              if(lons(i,j,is) .lt. -180.)then
                  lons(i,j,is) = lons(i,j,is) + 360.
              endif
              if(lons(i,j,is) .gt. +180.)then
                  lons(i,j,is) = lons(i,j,is) - 360.
              endif
          ENDDO
          ENDDO
          ENDDO

        ELSE ! Conformal grid ('plrstr','merctr',lambrt')
          call compute_latlon(nnxp,nnyp,n_staggers,mdlat,mdlon
     +                         ,deltax,xtn,ytn,lats,lons,istatus)
          if(istatus.ne.1)then
             print*,'Error returned: compute_latlon'
             return
          endif

        ENDIF

        IF (c6_maproj .ne. 'rotlat' .and. c6_maproj .ne. 'icshdr') THEN

        ns = istag

        print*,'ns [staggers array index] = ',ns
        if(ns.ne.1)
     &print*,' Using C-stagger grid for LSM fields and other output'

C*****************************************************************
        write(6,*)
        write(6,*)'Corner points.'

        call latlon_to_xy(lats(1,1,ns),lons(1,1,ns)      ,erad,xco,yco)
        write(6,701,err=702)1,1,   lats(1,1,ns)   ,lons(1,1,ns),xco,yco

        call latlon_to_xy(lats(nnxp,1,ns),lons(nnxp,1,ns),erad,xco,yco)
        write(6,701,err=702)nnxp,1,lats(nnxp,1,ns),lons(nnxp,1,ns)
     +                                                         ,xco,yco

        call latlon_to_xy(lats(1,nnyp,ns),lons(1,nnyp,ns),erad,xco,yco)
        write(6,701,err=702)1,nnyp,lats(1,nnyp,ns),lons(1,nnyp,ns)
     +                                                         ,xco,yco

        call latlon_to_xy(lats(nnxp,nnyp,ns),lons(nnxp,nnyp,1)
     +                                                   ,erad,xco,yco)
        write(6,701,err=702)nnxp,nnyp,lats(nnxp,nnyp,ns)
     +                                       ,lons(nnxp,nnyp,1),xco,yco

 701    format(' lat/lon/x/y at ',i5,',',i5,' =',2f12.5,2f12.0)
 702    continue

        call check_domain(lats(1,1,ns),lons(1,1,ns),nnxp,nnyp
     +,grid_spacing_m,1,istat_chk)
        if(istat_chk.ne.1)then
           print*,'Error returned from check_domain'
           istatus = istat_chk
           return
        endif

	ENDIF ! NOTE: totally avoided check_domain for rotlat & icshdr

       if(.true.)then
           print*,'get perimeter of grid'
           call get_domain_perimeter_grid(nnxp,nnyp,c10_grid_f
     1                  ,lats(1,1,ns),lons(1,1,ns)
     1                  ,0.0,rnorth,south,east,west,istatus)
           print*,'static dir = ',static_dir(1:lens)
           open(10,file=static_dir(1:lens)//'/llbounds.dat'
     +         ,status='unknown')

           print*,'write llbounds.dat'
           write(10,*)rnorth
           write(10,*)south
           write(10,*)east
           write(10,*)west
           close(10)
       endif

       write(6,*)'deltax = ',deltax
       write(6,*)'check_domain:status = ',istat_chk
c
C*****************************************************************
c calculate surface static fields
c
c ------------------------------------------------------
c ----------------- Terrain USGS 30 sec ----------------
c type = U
c
       itstatus=ishow_timer()
       print*
       print*,' Processing 30s topo data, l_fsf_gridgen = '
     1                                   ,l_fsf_gridgen

       allocate (topt_30(nnxp,nnyp),
     +           topt_30_s(nnxp,nnyp),
     +           topt_30_ln(nnxp,nnyp),
     +           topt_30_lt(nnxp,nnyp))

       l_topo_wps = l_fsf_gridgen

 600   if(.not. l_topo_wps .and. .not. l_topo_1s)then

        call s_len(path_to_topt30s,len)
        print*,'path to topt30s:        ',path_to_topt30s(1:len)
        path_to_topt30s(len+1:len+2)='/U'

        IF (c6_maproj .eq. 'rotlat') THEN
	 categorical=.false.

	 NCAT=1
	 allocate(adum3d(NNXP,NNYP,NCAT))
         allocate(adum2d(nnxp,nnyp))

         CALL alt_10by10_all(nnxp,nnyp,1200,grid_spacing_m/1000.,
     &   hlat,hlon,PATH_TO_TOPT30S,topt_30,adum3d,NCAT,topt_30_lt,
     &   topt_30_ln,topt_30_s,categorical,2,0)

	 write(6,*) 'topo before smooth'
         do J=nnyp,1,-nnyp/30
           write(6,633) (topt_30(I,J),I=1,nnxp,nnxp/15)
         enddo

         deallocate(adum3d)

  633    format(20(f5.0,1x))

!mp	 bogus status value to keep things smooth below
         istatus_30s=1
!mp

        ELSEIF (c6_maproj .eq. 'icshdr') THEN
         allocate  (GEODAT3D(nnxp,nnyp,1))
         call proc_geodat(nnxp,nnyp,1,path_to_topt30s
     +       ,lats(1,1,1),lons(1,1,1),data(1,1,1)
     +       ,GEODAT3D,istatus_30s)

        ELSE

         CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,1)
     +   ,ytn(1,1),deltax,deltay,TOPT_30,TOPT_30_S,TOPT_30_LN
     +   ,TOPT_30_LT,PATH_TO_TOPT30S,TOPTWVL,SILAVWT,new_DEM,1
     +   ,istatus_30s)

        ENDIF

        if(istatus_30s .ne. 1)then

          print*,'WARNING: File(s) missing for 30s terrain data'
          print*,' >>>> Try Processing 10m topo data <<<<'
          print*

          allocate (topt_10(nnxp,nnyp),
     +              topt_10_s(nnxp,nnyp),
     +              topt_10_ln(nnxp,nnyp),
     +              topt_10_lt(nnxp,nnyp))

          print*
          print*,' Calling GEODAT:  10m terrain data ....'

          CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,1)
     +,ytn(1,1),deltax,deltay,TOPT_10,TOPT_10_S,TOPT_10_LN
     +,TOPT_10_LT,PATH_TO_TOPT10M,TOPTWVL,SILAVWT,new_DEM,1
     +,istatus_10m)

          if(istatus_10m .ne. 1)then

             print*
             print*,'Error: File(s) missing for both 10m and
     +30s terrain data'
             print*,'Error: Aborting --- No static file created'
             return

          else

             print *,'topt_10    =',topt_10(1,1),topt_10(nnxp,nnyp)
             print*
             print*,'try blending 10min and 30s terrain data'
             call blend_topo(nnxp,nnyp,lats(1,1,1),lons(1,1,1)
     1,topt_10,topt_10_s,topt_10_ln,topt_10_lt
     1,topt_30,topt_30_s,topt_30_ln,topt_30_lt
     1,topt_out,topt_out_s,topt_out_ln,topt_out_lt)

             deallocate (topt_10,
     1                   topt_10_s,
     1                   topt_10_ln,
     1                   topt_10_lt)

          endif

        else ! go with 30s topo data since all tiles available

          print*,'Topo SW and NE corner values'
          print *,'topt_30    =',topt_30(1,1),topt_30(nnxp,nnyp)

          topt_out=topt_30
          topt_out_s=topt_30_s
          topt_out_ln=topt_30_ln
          topt_out_lt=topt_30_lt
          icount_30=nnyp*nnxp

        endif !istatus_30s ... data processed ok

       elseif(l_topo_1s .eqv. .true.)then
          r8lat(:,:) = lats(:,:,1)
	  r8lon(:,:) = lons(:,:,1)
          call get_topo_1s(nnxp,nnyp,grid_spacing_m,r8lat,r8lon,topt_out
     1                    ,path_to_topt30s,istatus)
          if(istatus .ne. 1)then
              write(6,*)' Error in get_topo_1s routine: returning'
              return
          endif

       elseif(.false.)then ! read topo data from wps output netCDF file
          write(6,*)' calling read_wrfstatic for wps topo'
          call s_len(path_to_topt30s,lenp)
          filename_wps = path_to_topt30s(1:lenp)//'/geo_em.d01.nc'
          call read_wrfstatic(nnxp,nnyp,lat,lon,filename_wps,topt_out
     1                       ,istatus)
          if(istatus .ne. 1)then
              write(6,*)' Error - no wrf static data: returning'
              return
          endif

       else                ! read topo data from LAPS/WRF FSF file
          call get_systime_i4(i4time_sys,istatus)
          if(istatus .ne. 1)then
              write(6,*)' Error - no systime: returning'
              return
          endif

          istatus = 0
          ntrys = 40
          itry = 0
          do while (itry .le. ntrys .AND. istatus .eq. 0)
              i4time_try = i4time_sys - (itry * laps_cycle_time)
              write(6,*)' Trying for FSF data at',itry,i4time_try
              call get_modelfg_2d(i4time_try,'TER',nnxp,nnyp,topt_out
     1                                                      ,istatus)
              itry = itry + 1
          enddo

          if(itry .eq. ntrys)then
              write(6,*)' Reached maximum number of tries'
              write(6,*)' Use geog topo data instead'
              l_topo_wps = .false.
              goto 600
          endif

       endif

       deallocate (topt_30,
     1             topt_30_s,
     1             topt_30_ln,
     1             topt_30_lt)

       print*
c
c ---------------------------------------------------------
c ------------------- 30s Landuse -------------------------
c type = V
c
       allocate  (GEODAT2D(nnxp,nnyp))
       allocate  (GEODAT3D(nnxp,nnyp,24))

        IF (c6_maproj .ne. 'rotlat') THEN

       if(.not.allocated(adum2d))allocate (adum2d(nnxp,nnyp))

       itstatus=ishow_timer()
       print*
       print*,' Calling GEODAT: Processing 30s landuse data.'
       print*,' Re-allocate GEODAT3D ',nnxp,nnyp,' 24'
       print*

       CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns)
     +,ytn(1,ns),deltax,deltay,GEODAT2D,GEODAT3D
     +,adum2d,adum2d,PATH_TO_LUSE_30S,2.0,0.0,new_DEM,24
     +,istatus)

       if(istatus.ne.1)then
         print*,'********* ERROR ***********'
         print*,'ERROR:  File(s) missing for landuse data'
         print*,'ERROR:  Static file not created'
         print*
         return
       endif

	ELSE

	categorical=.true.
        call alt_10by10_all(nnxp,nnyp,1200,grid_spacing_m/1000.,
     &		hlat,hlon,PATH_TO_LUSE_30S,GEODAT2D,
     &          GEODAT3D,24,adum2d,adum2d,adum2d,categorical,1,9999)

        do J=nnyp,1,-nnyp/25
        write(6,733) (INT(GEODAT2D(I,J)),I=1,nnxp,nnxp/15)
        enddo

  733   format(50(I2,1x))

	ENDIF

       ilndmsk=12  !land mask is 12 for both wrfsi and laps
       if(c10_grid_f(1:lf).eq.'wrfsi')then

	 IF (c6_maproj .ne. 'rotlat') THEN

c dominant category landuse
          ilnduse=11
          data(:,:,ilnduse)=GEODAT2D
c landmask for wrfsi
          data(:,:,ilndmsk) = 1.
          where(data(:,:,11) .eq. 16.)data(:,:,ilndmsk)=0.

c grids 15 thru 38 are percent distributions
          i=14
          do j=1,24
             data(:,:,i+j)=GEODAT3D(:,:,j)
          enddo
	ELSE ! rotlat case

        ilndmsk=7
        ilnduse=6

c landmask for wrfsi
           data(:,:,ilnduse)=GEODAT2D
           data(:,:,ilndmsk) = 1.
           where(data(:,:,ilnduse) .eq. 16.)data(:,:,ilndmsk)=0.

!!!! require a large % of water to make a water point
!!
	fraclk=0.7  !! require that 70% of points be water for output to be h2o

	DO JJ=1,nnyp
	 DO II=1,nnxp
	  IF (GEODAT3D(II,JJ,16) .gt. fraclk) THEN
	    data(II,JJ,ilndmsk)=0.
          ELSE
	    IF ( data(II,JJ,ilndmsk) .eq. 0.) THEN
!	      write(6,*) 'WATER WATER, BECOMING LAND:: ', II,JJ
	    ENDIF
	    data(II,JJ,ilndmsk)=1.
          ENDIF
         ENDDO
	ENDDO

!!	change dominant landuse  (otherwise will have land points with
!!	water as the dominant class)

	DO JJ=2,NNYP-1
        DO II=2,NNXP-1

	if (data(II,JJ,ilndmsk).eq.1.and.data(II,JJ,6).eq.16.) then
	onepos=amax1(data(II+1,JJ,6),data(II-1,JJ,6),
     +   data(II,JJ+1,6),data(II,JJ-1,6))
	twopos=amin1(data(II+1,JJ,6),data(II-1,JJ,6),
     +   data(II,JJ+1,6),data(II,JJ-1,6))

	if (twopos .ne. 16) then
	 	data(II,JJ,6)=twopos
	elseif (onepos .ne. 16) then
		data(II,JJ,6)=onepos
	else
!	leave data alone, revert the land mask to water
	write(6,*) 'switching back to water at II,JJ: ', II,JJ
		data(II,JJ,ilndmsk)=0.
	endif

	endif

	ENDDO
	ENDDO

! putting land-use percentages in 19-42 for NMM
        I=18
           do j=1,24
              data(:,:,i+j)=GEODAT3D(:,:,j)
           enddo

        ENDIF

       else
c dominant category landuse
          ilnduse=5
          data(:,:,ilnduse)=GEODAT2D   !landuse for laps
c landmask for laps
          data(:,:,ilndmsk)=1.
          where(data(:,:,5) .eq. 16.)data(:,:,ilndmsk)=0.
       endif
c
c **********************************************
c fix water areas to have 0.0m terrain elevation
c currently no fix going on!
c **********************************************
       lforce_ter_2_zero=.false.
       if(lforce_ter_2_zero)then
          where(data(:,:,ilndmsk).eq.0)topt_out=0.0
       endif

c for WRFSI we'll bilinearly interpolate to get the topo
c from the non-staggered topo

       if(c10_grid_f(1:lf).eq.'wrfsi')then

	IF (c6_maproj .ne. 'rotlat') THEN

          data(:,:,9)=topt_out
          call move(topt_out_s,data(1,1,48),nnxp,nnyp)
          call move(topt_out_ln,data(1,1,49),nnxp,nnyp)
          call move(topt_out_lt,data(1,1,50),nnxp,nnyp)

          indxter=51
          data(:,:,indxter)=r_missing_data
          do j=2,nnyp
          do i=2,nnxp
             call bilinear_interp(i,j,nnxp,nnyp,topt_out,result)
             if(data(i,j,51).gt.30000. .and.
     & data(i,j,51).lt.r_missing_data)then
                print*,'Here: Large topo: i/j/topo ',i,j,data(i,j,51)
             endif
             data(i-1,j-1,indxter)=result
          enddo
          enddo

          where((data(:,:,indxter).lt. 0.01).and.
     1          (data(:,:,indxter).gt.-0.01))data(:,:,indxter)=0.0

	ELSE

!       for rotlat, topt_out is already on mass point stagger.
!       No bilinear interp needed

!!	boundary smooth the topo data

	write(6,*) 'call smdhld'
	write(6,*) 'NNXP, NNYP: ', NNXP,NNYP
	call smdhld(nnxp,nnyp,topt_out,1.-data(:,:,ilndmsk),12,12)

c       write(6,*) 'topo after smoothing'
c       do J=nnyp,1,-nnyp/30
c       write(6,633) (topt_out(I,J),I=1,nnxp,nnxp/15)
c       enddo

!-------------4-point averaging of mountains along inner boundary-------

          do i=1,NNXP-1
      topt_out(i,2)=0.25*(topt_out(i,1)+topt_out(i+1,1)+
     +                    topt_out(i,3)+topt_out(i+1,3))
          enddo

          do i=1,NNXP-1
      topt_out(i,NNYP-1)=0.25*(topt_out(i,NNYP-2)+topt_out(i+1,NNYP-2)+
     +                         topt_out(i,NNYP)+topt_out(i+1,NNYP))
          enddo

          do j=4,NNYP-3,2
      topt_out(1,j)=0.25*(topt_out(1,j-1)+topt_out(2,j-1)+
     +                    topt_out(1,j+1)+topt_out(2,j+1))
          enddo

          do j=4,NNYP-3,2
      topt_out(NNXP,j)=0.25*(topt_out(NNXP-1,j-1)+topt_out(NNXP,j-1)+
     +                         topt_out(NNXP-1,j+1)+topt_out(NNXP,j+1))
          enddo

	indxter=18
        do J=1,nnyp
        do I=1,nnxp
          data(i,j,indxter)=topt_out(i,j)
        enddo
        enddo

          where((data(:,:,indxter).lt. 0.01).and.
     1          (data(:,:,indxter).gt.-0.01))data(:,:,indxter)=0.0

	ENDIF

       else  !must be LAPS

          indxter=3
          call move(topt_out,data(1,1,indxter),nnxp,nnyp)         ! KWD
          call move(topt_out_s,data(1,1,7),nnxp,nnyp)             ! JS
          call move(topt_out_ln,data(1,1,8),nnxp,nnyp)            ! JS
          call move(topt_out_lt,data(1,1,9),nnxp,nnyp)            ! JS

       endif

       print*
       print *,'terrain  = ',data(1,1,indxter),data(nnxp,nnyp,indxter)

c      print *,'# of grid pts using 30 sec terrain data =  ',icount_30
c      print *,'# of grid pts using blended terrain data = ',icount_ramp
       print*
C
C Output Section for topo, lat-lons and corners
C -----------------------------------------------
       if(c10_grid_f(1:lf).eq.'wrfsi')then
         if(c6_maproj.ne.'rotlat')then
          write(cnest,'(i2.2)')nest
          clatlon='latlon.d'//cnest//'.dat'            !binary non-staggered lat-lon file
          clatlons='latlon-mass.d'//cnest//'.dat'      !binary staggered lat-lon file
          ctopo='topo.d'//cnest//'.dat'                !binary non-staggered topo file
          ctopos='topo-mass.d'//cnest//'.dat'          !binary staggered topo file
          ctopog='topography.d'//cnest//'.dat'         !ascii non-staggered topo file
          ctopogs='topography-mass.d'//cnest//'.dat'   !ascii staggered topo file
          clatlon2d='latlon2d.d'//cnest//'.dat'        !ascii lat-lon non-staggered file
          clatlons2d='latlon2d-mass.d'//cnest//'.dat'  !ascii lat-lon staggered file
          corners='corners.d'//cnest//'.dat'           !non-staggered corners
          cornerss='corners-mass.d'//cnest//'.dat'     !staggered corners
         else
          clatlon='latlon-rotlat.dat'        !binary lat/lon file
          clatlon2d='latlon2d-rotlat.dat'    !ascii lat/lon file
          ctopo='topo-rotlat.dat'            !binary topo file
          ctopog='topography-rotlat.dat'     !ascii topo file
          corners='corners-rotlat.dat'       !domain corners file
         endif
       else
          clatlon='latlon.dat'        !binary lat/lon file
          clatlon2d='latlon2d.dat'    !ascii lat/lon file
          ctopo='topo.dat'            !binary topo file
          ctopog='topography.dat'     !ascii topo file
          corners='corners.dat'       !domain corners file
       endif

       call s_len(static_dir,len)
       if(c10_grid_f(1:lf).eq.'wrfsi')then

         if(c6_maproj.ne.'rotlat')then
           in=0
c write both non-stagger and stagger data for wrfsi, non-rotlat
           do k=1,2
            if(k.eq.1)then
             is=1
             idxt=9
             filename=static_dir(1:len)//clatlon
             open(10,file=filename,status='unknown',form='unformatted')
             filename=static_dir(1:len)//ctopo
             open(11,file=filename,status='unknown',form='unformatted')
             filename=static_dir(1:len)//corners
             open(15,file=filename,status='unknown')
             filename=static_dir(1:len)//ctopog
             open(666,file=filename)
             filename=static_dir(1:len)//clatlon2d
             open(667,file=filename)
            else
             is = 4
             in = 1
             idxt=51
             filename=static_dir(1:len)//clatlons
             open(10,file=filename,status='unknown',form='unformatted')
             filename=static_dir(1:len)//ctopos
             open(11,file=filename,status='unknown',form='unformatted')
             filename=static_dir(1:len)//cornerss
             open(15,file=filename,status='unknown')
             filename=static_dir(1:len)//ctopogs
             open(666,file=filename)
             filename=static_dir(1:len)//clatlons2d
             open(667,file=filename)
            endif

            write(10)lats(1:nnxp-in,1:nnyp-in,is)
     +,lons(1:nnxp-in,1:nnyp-in,is)
            close(10)
            write(11)data(1:nnxp-in,1:nnyp-in,idxt)
            close(11)
            write(15,*)lats(1,1,is),lons(1,1,is)
            write(15,*)lats(1,nnyp-in,is),lons(1,nnyp-in,is)
            write(15,*)lats(nnxp-in,1,is),lons(nnxp-in,1,is)
            write(15,*)lats(nnxp-in,nnyp-in,is),lons(nnxp-in,nnyp-in,is)
            close(15)

            do j=1,nnyp-in
            do i=1,nnxp-in
               write(666,*) data(i,j,idxt)
               write(667,*) lats(i,j,is),lons(i,j,is)
            enddo
            write(666,'()')
            write(667,'()')
            enddo
            close(666)
            close(667)

           enddo

         else   !if(c6_maproj.eq.'rotlat')then

           is=1
           idxt=18
           filename=static_dir(1:len)//clatlon
           open(10,file=filename,status='unknown',form='unformatted')
           filename=static_dir(1:len)//ctopo
           open(11,file=filename,status='unknown',form='unformatted')
           filename=static_dir(1:len)//corners
           open(15,file=filename,status='unknown')
           filename=static_dir(1:len)//ctopog
           open(666,file=filename)
           filename=static_dir(1:len)//clatlon2d
           open(667,file=filename)

           write(10)lats(1:nnxp,1:nnyp,is),lons(1:nnxp,1:nnyp,is)
           close(10)
           write(11)data(:,:,idxt)
           close(11)
           write(15,*)lats(1,1,is),lons(1,1,is)
           write(15,*)lats(1,nnyp,is),lons(1,nnyp,is)
           write(15,*)lats(nnxp,1,is),lons(nnxp,1,is)
           write(15,*)lats(nnxp,nnyp,is),lons(nnxp,nnyp,is)
           close(15)

           do j=1,nnyp
           do i=1,nnxp
              write(666,*) data(i,j,idxt)
              write(667,*) lats(i,j,is),lons(i,j,is)
           enddo
           write(666,'()')
           write(667,'()')
           enddo
           close(666)
           close(667)

         endif

       else !laps case

         is=1
         idxt=3
!        filename=static_dir(1:len)//clatlon
!        open(10,file=filename,status='unknown',form='unformatted')
!        filename=static_dir(1:len)//ctopo
!        open(11,file=filename,status='unknown',form='unformatted')
         filename=static_dir(1:len)//corners
         open(15,file=filename,status='unknown')
         filename=static_dir(1:len)//ctopog
         open(666,file=filename)
         filename=static_dir(1:len)//clatlon2d
         open(667,file=filename)

!        write(10)lats(1:nnxp,1:nnyp,is),lons(1:nnxp,1:nnyp,is)
!        close(10)
!        write(11)data(:,:,idxt)
!        close(11)
         write(15,*)lats(1,1,is),lons(1,1,is)
         write(15,*)lats(1,nnyp,is),lons(1,nnyp,is)
         write(15,*)lats(nnxp,1,is),lons(nnxp,1,is)
         write(15,*)lats(nnxp,nnyp,is),lons(nnxp,nnyp,is)
         close(15)

         do j=1,nnyp
         do i=1,nnxp
            write(666,*) data(i,j,idxt)
            write(667,*) lats(i,j,is),lons(i,j,is)
         enddo
         write(666,'()')
         write(667,'()')
         enddo
         close(666)
         close(667)

       endif
c ----------------------------------------------------------
c
c ------ Land Fraction from 30s Land Use -------------------
c type = L (no longer used ... no land fraction data base)
c
        print*
        print*,' Create from 30s land use fractional dist'
        print*

        GEODAT2D(:,:)=1.- GEODAT3D(:,:,16)
c
c this function allows variable amounts of smoothing of
c the land fraction data depending on grid spacing
c
!       iter=min(10,int(8000./deltax))
        iter = 0

	IF (c6_maproj .ne. 'rotlat') THEN

          print*,'Filter land fraction with 2dx ',iter
          do k=1,iter
             call filter_2dx(GEODAT2D,nnxp,nnyp,1, 0.5)
             call filter_2dx(GEODAT2D,nnxp,nnyp,1,-0.5)
             GEODAT2D = max(GEODAT2D,0.)
          enddo

	ENDIF

        if(c10_grid_f(1:lf).eq.'wrfsi')then

	    IF (c6_maproj .ne. 'rotlat') THEN
		idx=10
	    ELSE
		idx=5
	    ENDIF

	else ! not wrfsi
		idx=4
	endif

        data(:,:,idx)=GEODAT2D
        print *,'pctlandfrac=',GEODAT2D(1,1),GEODAT2D(nnxp,nnyp)
c
c -------------------------------------------------------------------
c potential fix of fictitous islands for certain resolution domains.
c Story here is that west coast terrain can be "funky" possibly due
c to steepness and method to compute avg terrain.

	write(6,*) 'idx, indxter: ', idx, indxter
        where(data(:,:,idx).le.0.1 .and. data(:,:,indxter).lt.5.0)
     &data(:,:,indxter)=0.0

c
c ------------------------------------------------------------
c -------------- Top Layer Soil Type -------------------------
c type = O
c
        itstatus=ishow_timer()
        print*
        print*,' Processing 30s soil type top layer data....'

        deallocate(GEODAT3D)
        allocate (GEODAT3D(nnxp,nnyp,16))

        IF (c6_maproj .ne. 'rotlat') THEN

        CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns),ytn(1,ns)
     +,deltax,deltay,GEODAT2D,GEODAT3D,adum2d,adum2d
     +,PATH_TO_SOILTYPE_TOP_30S,2.0,0.0,new_DEM,16   !maxdatacat
     +,istatus_soil)

        if(istatus_soil.ne.1)then
           print*
           print*,' Soil type data not processed completely'
           if(c10_grid_f(1:lf).eq.'wrfsi')then
             print*,' File(s) missing for soil type top layer data'
             print*,' Error:  Static file not created'
             print*
             return
           else
             print*,' File(s) missing for soil type top layer data'
             print*,'           *** WARNING ***'
             print*,' Soil Type Top data not added to static file'
             print*
           endif
           GEODAT2D=r_missing_data
           GEODAT3D=r_missing_data
        endif

	ELSE

	categorical=.true.
        call alt_10by10_all(nnxp,nnyp,1200,grid_spacing_m/1000.,
     &		hlat,hlon,PATH_TO_SOILTYPE_TOP_30S,GEODAT2D,
     &          GEODAT3D,16,adum2d,adum2d,adum2d,categorical,1,9999)

        write(6,*) 'soiltype_top'
        do J=nnyp,1,-nnyp/25
        write(6,733) (INT(GEODAT2D(I,J)),I=1,nnxp,nnxp/15)
        enddo

        ENDIF

        if(c10_grid_f(1:lf).eq.'wrfsi')then

         IF (c6_maproj .ne. 'rotlat') THEN
c
c make water points = 0 for adjust_geog. later we'll put it back to original
c
           where(GEODAT2D.eq.14)GEODAT2D=0.0

           call adjust_geog(nnxp-1,nnyp-1,1,'soiltype',istatus_soil
     &,lats(1:nnxp-1,1:nnyp-1,ns),data(1:nnxp-1,1:nnyp-1,indxter)
     &,data(1:nnxp-1,1:nnyp-1,ilndmsk),GEODAT2D(1:nnxp-1,1:nnyp-1)
     &,istatus)

           where(GEODAT2D.eq.0.0)GEODAT2D=14.
           idxstl=13
           data(:,:,idxstl)=GEODAT2D
           i=51
           do j=1,16
              data(:,:,i+j)=GEODAT3D(:,:,j)
           enddo
	 ELSE
           idxstl=8
           data(:,:,8)=GEODAT2D
           I=42
           do j=1,16
              data(:,:,i+j)=GEODAT3D(:,:,j)
           enddo
         ENDIF

        else
c
c make water points = 0 for adjust_geog. After put it back to original
c
           where(GEODAT2D.eq.14)GEODAT2D=0.0

           call adjust_geog(nnxp,nnyp,1,'soiltype',istatus_soil
     &,lats(1,1,ns),data(1,1,indxter),data(1,1,ilndmsk),GEODAT2D
     &,istatus)

           where(GEODAT2D.eq.0.0)GEODAT2D=14.
           idxstl=10
           data(:,:,idxstl)=GEODAT2D

        endif
c
c --------------------------------------------------------------
c -------------- Bottom Layer Soil Type ------------------------
c type = O
c
        itstatus=ishow_timer()
        print*
        print*,' Processing 30s soil type bottom layer data'

        IF (c6_maproj .ne. 'rotlat') THEN

        CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns)
     +,ytn(1,ns),deltax,deltay,GEODAT2D,GEODAT3D,adum2d,adum2d
     +,PATH_TO_SOILTYPE_BOT_30S,2.0,0.0,new_DEM,16   !maxdatacat
     +,istatus_soil)

        if(istatus_soil.ne.1)then
           print*
           print*,' Bottom layer soil data not processed completely'
           if(c10_grid_f(1:lf).eq.'wrfsi')then
             print*,' File(s) missing for soil type bot layer data'
             print*,' Error:  Static file not created'
             print*
             return
           else
             print*,' File(s) missing for soil type bot layer data'
             print*,'           *** WARNING ***'
             print*,' Soil Type Bot data not added to static file'
             print*
           endif
           GEODAT2D=r_missing_data
           GEODAT3D=r_missing_data
        endif

	ELSE

	categorical=.true.
        call alt_10by10_all(nnxp,nnyp,1200,grid_spacing_m/1000.,
     &		hlat,hlon,PATH_TO_SOILTYPE_BOT_30S,GEODAT2D,
     &          GEODAT3D,16,adum2d,adum2d,adum2d,categorical,1,9999)

        ENDIF
c
c soiltype bottom layer ... 68 thru 83
c
        if(c10_grid_f(1:lf).eq.'wrfsi')then

          IF (c6_maproj .ne. 'rotlat') THEN
c
c make water points = 0 for adjust_geog. After, put it back to original.
c
            where(GEODAT2D.eq.14)GEODAT2D=0.0

            call adjust_geog(nnxp-1,nnyp-1,1,'soiltype',istatus_soil
     &,lats(1:nnxp-1,1:nnyp-1,ns),data(1:nnxp-1,1:nnyp-1,indxter)
     &,data(1:nnxp-1,1:nnyp-1,ilndmsk),GEODAT2D(1:nnxp-1,1:nnyp-1)
     &,istatus)

            where(GEODAT2D.eq.0.)GEODAT2D=14.0
            idxsbl=14
            data(:,:,idxsbl)=GEODAT2D
            i=67
            do j=1,16
               data(:,:,i+j)=GEODAT3D(:,:,j)
            enddo
	  ELSE
            idxsbl=9
	    data(:,:,9)=GEODAT2D
            i=58
            do j=1,16
               data(:,:,i+j)=GEODAT3D(:,:,j)
            enddo
          ENDIF

        else
c
c make water points = 0 for adjust_geog. After, put it back to original.
c
            where(GEODAT2D.eq.14)GEODAT2D=0.0

            call adjust_geog(nnxp,nnyp,1,'soiltype',istatus_soil
     &,lats(1,1,ns),data(1,1,indxter),data(1,1,ilndmsk),GEODAT2D
     &,istatus)
            where(GEODAT2D.eq.0.)GEODAT2D=14.0
            idxsbl=11
            data(:,:,idxsbl)=GEODAT2D
        endif
c
c ---------------------------------------------------------
c ----------------- Greenness Fraction --------------------
c type = G
c
        itstatus=ishow_timer()
        print*
        print*,' Calling PROC_GEODAT: Process green frac data.'
        print*,' Re-allocate GEODAT3D ',nnxp,nnyp,' 12'
        print*

        deallocate(GEODAT3D)
        allocate  (GEODAT3D(nnxp,nnyp,12),stat=istat)
        if(istat.ne.0)then
           print*,'Error allocating data object GEODAT3D ',istat
           if(c10_grid_f(1:lf).eq.'wrfsi')then
              print*,'Error: Aborting process'
              print*,'Error: Maybe not enough memory'
              print*,'Error: Cannot allocate GEODAT3D'
              stop
           else
              print*,'Warning: Continue without green fraction'
           endif
        endif

	IF (c6_maproj .ne. 'rotlat') THEN

        print*,'Calling proc_geodat'
        call proc_geodat(nnxp,nnyp,12,path_to_green_frac
     +,lats(1,1,ns),lons(1,1,ns),data(1,1,ilndmsk)
     +,GEODAT3D,istatus_grn)

	ELSE

	categorical=.FALSE.
	useland=.TRUE.
        call alt_hemi_all(nnxp,nnyp,hlat,hlon,
     &                       PATH_TO_GREEN_FRAC,
     &                       type,12,grid_spacing_m/1000.,GEODAT3D,
     &                       GEODAT2D,categorical,data(1,1,ilndmsk),
     &			     useland)

!	convert from 0-100% to 0-1 fraction

	do N=1,12
	do J=1,NNYP
	do I=1,NNXP
	    GEODAT3D(I,J,N)=GEODAT3D(I,J,N)/100.
	enddo
	enddo
	enddo

        istatus_grn=1 !bogus value

        ENDIF

        if(istatus_grn.ne.1)then
         print*
         print*,'greenness fraction data not processed completely'
         if(c10_grid_f(1:lf).eq.'wrfsi')then
            print*,' Error: File(s) missing for green frac data'
            print*,' Error: Static file not created'
            print*
            istatus=0
            return
         else
            print*,' Warning: File(s) missing for green frac data'
            print*,' Warning: green frac not added to static file'
            print*
         endif
         GEODAT3D=r_missing_data
        endif

        if(c10_grid_f(1:lf).eq.'wrfsi')then

           IF (c6_maproj .ne. 'rotlat') THEN

              igrn=83

c             call adjust_geog(nnxp-1,nnyp-1,12,'greenfrac',istatus_grn
c    &,lats(1:nnxp-1,1:nnyp-1,ns),data(1:nnxp-1,1:nnyp-1,indxter)
c    &,data(1:nnxp-1,1:nnyp-1,ilndmsk),GEODAT3D(1:nnxp-1,1:nnyp-1,:)
c    &,istatus)

           ELSE
               igrn=74
           ENDIF

        else

           igrn=12

           call adjust_geog(nnxp,nnyp,igrn,'greenfrac',istatus_grn
     &,lats(:,:,ns),data(:,:,indxter),data(:,:,ilndmsk),GEODAT3D
     &,istatus)
        endif
        do j=1,12
           data(:,:,igrn+j)=GEODAT3D(:,:,j)
        enddo
c
c annual max/min greenness fraction in domain
c --------------------------------------------
        if(c10_grid_f(1:lf).eq.'wrfsi')then
           dommaxgf=0.0
           dommingf=99.0
           print*,' compute max/min greenness frac at grid points'
           IF (c6_maproj .ne. 'rotlat') THEN
               mxgiwrf=110
               mngiwrf=111
           ELSE
               mxgiwrf=101
               mngiwrf=102
           ENDIF
           do j=1,nnyp
           do i=1,nnxp
              data(i,j,mxgiwrf)=MAXVAL(GEODAT3D(i,j,:))
              data(i,j,mngiwrf)=MINVAL(GEODAT3D(i,j,:))
              if(data(i,j,mxgiwrf).gt.dommaxgf)
     &dommaxgf=data(i,j,mxgiwrf)
              if(data(i,j,mngiwrf).lt.dommingf)
     &dommingf=data(i,j,mngiwrf)
           enddo
           enddo
           write(6,*)'Domain Annual Max/Min Green Fraction',
     &dommaxgf,dommingf
        endif
c
c ---------------------------------------------------------
c --------------- Deep Soil Temperature -------------------
c type = T
c
	IF (c6_maproj .ne. 'rotlat') THEN
        itstatus=ishow_timer()
        print*
        print*,' Call PROC_GEODAT: Process 1 deg soiltemp data.'

        GEODAT2D=0.0

        call proc_geodat(nnxp,nnyp,1
     1,path_to_soiltemp_1deg,lats(1,1,ns),lons(1,1,ns)
     1,data(1,1,ilndmsk),GEODAT2D,istatus_tmp)

	ELSE

	write(6,*) 'call alt_10by10 for soiltemp'
        GEODAT2D=0.0
	categorical=.false.

	NCAT=1
        call alt_10by10_all(nnxp,nnyp,10,grid_spacing_m/1000.,
     &		hlat,hlon,PATH_TO_SOILTEMP_1DEG,GEODAT2D,
     &          ADUM3D,NCAT,adum2d,adum2d,adum2d,categorical,2,0)

!       convert to degrees

	write(6,*) 'convert to degrees ', nnxp,nnyp
        DO J=1,NNYP
        DO I=1,NNXP
                GEODAT2D(I,J)=GEODAT2D(I,J)/100.
        ENDDO
        ENDDO

        istatus_tmp=1 !bogus

        ENDIF

!       do J=NNYP,1,-NNYP/25
!         write(6,833) (GEODAT2D(I,J),I=1,NNXP,NNXP/15)
!       enddo

        if(istatus_tmp.ne.1)then
         print*
         print*,'soiltemp data not processed completely'
         if(c10_grid_f(1:lf).eq.'wrfsi')then
            print*,' Error:  File(s) missing for soiltemp data'
            print*,' Error:  Static file not created'
            print*
            istatus=0
            return
         else
            print*,' File(s) missing for soiltemp data'
            print*,'           *** WARNING ***'
            print*,' soiltemp data not added to static file'
            print*
         endif
         GEODAT2D=r_missing_data
        endif

        idst=25
        if(c10_grid_f(1:lf).eq.'wrfsi')then

            IF (c6_maproj .ne. 'rotlat') THEN

                idst=96

                call adjust_geog(nnxp-1,nnyp-1,1,'soiltemp',istatus_tmp
     &,lats(1:nnxp-1,1:nnyp-1,ns),data(1:nnxp-1,1:nnyp-1,indxter)
     &,data(1:nnxp-1,1:nnyp-1,ilndmsk),GEODAT2D(1:nnxp-1,1:nnyp-1)
     &,istatus)

            ELSE
                idst=87
            ENDIF

        else

            call adjust_geog(nnxp,nnyp,1,'soiltemp',istatus_grn
     &,lats(:,:,ns),data(:,:,indxter),data(:,:,ilndmsk),GEODAT2D
     &,istatus)

        endif
	write(6,*) 'idst= ', idst
        data(:,:,idst)=GEODAT2D
	write(6,*) 'idst written'
c
c -------------- Terrain slope index categories -------------------------
c type = I
c
        itstatus=ishow_timer()
        print*
        print*,' Processing 1 deg islope data'

        allocate (adum3d(nnxp,nnyp,9))  !temporarily holds %dist of islope cat

	IF (c6_maproj .ne. 'rotlat') THEN

        CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns)
     +,ytn(1,ns),deltax,deltay,GEODAT2D,adum3d,adum2d,adum2d
     +,PATH_TO_ISLOPE,2.0,0.0,new_DEM,9,istatus_slp)

	ELSE

	categorical=.false.
!!	should be done as n.n. interp by alt_interp code
	useland=.FALSE.
        call alt_hemi_all(nnxp,nnyp,hlat,hlon,
     &                       PATH_TO_ISLOPE,
     &                       type,1,grid_spacing_m/1000.,ADUM3D,
     &                       GEODAT2D,categorical,data(:,:,ilndmsk),
     &			     useland)

	istatus_slp=1 ! bogus

	ENDIF

        deallocate (adum3d)

        if(istatus_slp.ne.1)then
           print*
           print*,' ter slope category data not processed completely'
           if(c10_grid_f(1:lf).eq.'wrfsi')then
             print*,' Error: File(s) missing for islope data'
             print*,' Error: Static file not created'
             print*
             return
           else
             print*,' Warning: File(s) missing for islope data'
             print*,' Warning: islope not added to static file'
             print*
           endif
           GEODAT2D=r_missing_data
           GEODAT3D=r_missing_data
        endif

c put the categories back to the original raw data. if it is a land
c point but islope indicates water, force islope = 1.
        where(GEODAT2D .eq. 8)GEODAT2D=13.0
        where(GEODAT2D .eq. 9)GEODAT2D=0.0

        if(c10_grid_f(1:lf).eq.'wrfsi')then

	   IF (c6_maproj .ne. 'rotlat') THEN

               islp=109
               call adjust_geog(nnxp-1,nnyp-1,1,'islope',istatus_slp
     &,lats(1:nnxp-1,1:nnyp-1,ns),data(1:nnxp-1,1:nnyp-1,indxter)
     &,data(1:nnxp-1,1:nnyp-1,ilndmsk),GEODAT2D(1:nnxp-1,1:nnyp-1)
     &,istatus)

           ELSE
               islp=100 ! will need more adjustments down the road for this
           ENDIF

        else

           islp=38
           call adjust_geog(nnxp,nnyp,1,'islope',istatus_slp
     &,lats(1,1,ns),data(1,1,indxter),data(1,1,ilndmsk),GEODAT2D
     &,istatus)

        endif

        if(istatus.ne.1)then
           print*,' Warning: Processing incomplete: adjust_geog_data'
           if(c10_grid_f(1:lf).eq.'wrfsi')then
              print*,'Error: wrfsi static file not generated'
              return
           endif
        endif

        data(:,:,islp)=GEODAT2D

c force land points to have the correct (default) islope value.
        where(data(:,:,islp)   .eq. 0.0 .and.
     .        data(:,:,ilndmsk).eq. 1.0)data(:,:,islp)=1.0
c force water points to have the correct category for islope
        where(data(:,:,ilndmsk).eq. 0.0)data(:,:,islp)=0.0

!write(6,*) 'ISLOPE: '
!      do J=nnyp,1,-nnyp/30
!       write(6,733) (int(data(I,J,islp)),I=1,nnxp,nnxp/20)
!       enddo
c
c    subroutine adjust_geog_data now named adjust_geog and called
c    after each data type is processed (as opposed to all geog data
c    getting adjusted with one call).
c
c!JS: even though Matt Pyle suggests that adjust geog is fine for nmm
c     it currently isn't used.  If it should be used, then it would be
c     similar to the LAPS case which is not a staggered grid.
c
!mp	believe adjust_geog_data should be fine for nmm case
c
c -------------------------------------------------------------------
c ---------------Monthly Albedo Climatology--------------------------
c type = A
c
        itstatus=ishow_timer()
        print*
        print*,' Processing albedo climo data'

	IF (c6_maproj .ne. 'rotlat') THEN

        call proc_geodat(nnxp,nnyp,12,path_to_albedo
     +,lats(1,1,ns),lons(1,1,ns),data(1,1,ilndmsk)
     +,GEODAT3D,istatus_alb)

	ELSE

	categorical=.FALSE.
	useland=.TRUE.
        call alt_hemi_all(nnxp,nnyp,hlat,hlon,
     &                       PATH_TO_ALBEDO,
     &                       type,12,grid_spacing_m/1000.,GEODAT3D,
     &                       ADUM2D,categorical,data(1,1,ilndmsk),
     &			     useland)

        istatus_alb=1 !bogus

	do N=1,12
	 do J=1,NNYP
	  do I=1,NNXP
	   GEODAT3D(I,J,N)=GEODAT3D(I,J,N)*.01
	  enddo
	 enddo
	enddo

        ENDIF

        print*,'Done in proc_geodat: albedo'

        if(istatus_alb.ne.1)then
         print*
         print*,' Error: Albedo climo data not processed completely'
         print*,' Error: File(s) missing for albedo data'
         print*,' Error: Static file not created: albedo missing'
         print*
         istatus=0
         return
        endif

c force water points to 0.08 albedo

        do k=1,12
           where(data(:,:,ilndmsk).eq.0.0)GEODAT3D(:,:,k)=0.08
        enddo

        if(c10_grid_f(1:lf).eq.'wrfsi')then
        IF (c6_maproj .ne. 'rotlat') THEN
           ialb=96
	ELSE
	   ialb=87
	ENDIF
        else
           ialb=25
        endif
        do j=1,12
           data(:,:,ialb+j)=GEODAT3D(:,:,j)
        enddo
c
c ---------------- Max Snow Albedo ------------------
c type = M
c
	IF (c6_maproj .ne. 'rotlat') THEN

        call proc_geodat(nnxp,nnyp,1,path_to_maxsnoalb
     +,lats(1,1,ns),lons(1,1,ns),data(1,1,ilndmsk)
     +,GEODAT2D,istatus_alb)

	ELSE

	categorical=.false.

	useland=.TRUE.
	allocate(adum3d(nnxp,nnyp,1))
        call alt_hemi_all(nnxp,nnyp,hlat,hlon,
     &                       PATH_TO_MAXSNOALB,
     &                       type,1,grid_spacing_m/1000.,ADUM3D,
     &                       GEODAT2D,categorical,data(1,1,ilndmsk),
     &			     useland)

	do J=1,NNYP
	do I=1,NNXP
	GEODAT2D(I,J)=GEODAT2D(I,J)*.01
	enddo
	enddo

	ENDIF

        print*,'Done in proc_geodat: Max Snow Albedo'

        if(istatus_alb.ne.1)then
         print*
         if(c10_grid_f(1:lf).eq.'wrfsi')then
            print*,'--------------- WRFSI ------------------'
            print*,'Error: Max Snow Albedo not processed completely'
            print*,'Error: Static file not created '
            print*
            istatus=0
            return
         else
            print*,'Warning: Max Snow Albedo not processed completely'
         endif
        endif
c
c force max albedo = 0.08 over water. force max albedo = 0.7 over ice
c
        where(data(:,:,ilndmsk).eq.0.0)GEODAT2D=0.08
        where(data(:,:,ilnduse).eq.24.0)GEODAT2D=0.7

        if(c10_grid_f(1:lf).eq.'wrfsi')then

	IF (c6_maproj .ne. 'rotlat') THEN
           mxsnalb=47
           data(:,:,mxsnalb)=GEODAT2D
	ELSE
           mxsnalb=14
	   data(:,:,mxsnalb)=GEODAT2D
	ENDIF

        else
           mxsnalb=6
           data(:,:,mxsnalb)=GEODAT2D
        endif

        deallocate (GEODAT2D)
        deallocate (GEODAT3D)
c
c ---------------------------------------------------------------------------------
c Let's compare the grids to landmask to assess their consistency (or lack thereof)
c ---------------------------------------------------------------------------------

       if(c10_grid_f(1:lf).eq.'wrfsi')then
        print*,'Total number of gridpoints (nx*ny) = ',nnxp*nnyp
        print*
        do i=1,10
           if(i == 1)ii=indxter  !terrain
           if(i == 2)ii=idxstl   !soil top layer
           if(i == 3)ii=idxsbl   !soil bot layer
           if(i == 4)ii=igrn+6   !mxgiwrf  !max greenness
           if(i == 5)ii=mngiwrf  !min greenness
           if(i == 6)ii=idst     !deep soil temp
           if(i == 7)ii=islp     !terrain slope index
           if(i == 8)ii=ialb+6   !albedo; arbitrarily at month 6
           if(i == 9)ii=mxsnalb  !max snow albedo
           if(i == 10)ii=ilnduse !landuse dominant category
           call gridcompare(nnxp,nnyp,i,data(1,1,ii),data(1,1,ilndmsk)
     &,istatus)
        enddo
       endif
c
c ---------------------------------------------------------------------------------
c --------- This is were we prepare for writing the netCDF static file ------------
c ---------------------------------------------------------------------------------

        if(.not.allocated(var))allocate (var(ngrids)
     +,comment(ngrids))

        if(c10_grid_f(1:lf).eq.'wrfsi')then

	IF (c6_maproj .ne. 'rotlat') THEN

           in1=0
           in2=0
           do j=1,ns
              in1=in2+1
              in2=in1+1
              data(:,:,in1)=lats(:,:,j)
              data(:,:,in2)=lons(:,:,j)
           enddo

c compute map/grid information. Reallocate GEODAT3D array
           allocate (GEODAT3D(nnxp,nnyp,4))
           call get_projrot_grid(nnxp,nnyp,lats(1,1,ns)
     +,lons(1,1,ns),GEODAT3D,istatus)

           i=39

           data(:,:,i)  =GEODAT3D(:,:,1)
           data(:,:,i+1)=GEODAT3D(:,:,2)

           call get_map_factor_grid(nnxp,nnyp,n_staggers
     +,lats,lons ,GEODAT3D,istatus)
           if(istatus.ne.1)then
              print*,'Error returned: get_maps_factor_grid'
              return
           endif

           do j=1,ns
              data(:,:,i+j+1)=GEODAT3D(:,:,j)
           enddo
c
           call get_coriolis_components(nnxp,nnyp,lats(1,1,ns)
     +,GEODAT3D)
           data(:,:,i+6)=GEODAT3D(:,:,1)
           data(:,:,i+7)=GEODAT3D(:,:,2)

           deallocate(GEODAT3D)

c          call move(static_albedo,data(1,1,i+8),nnxp,nnyp)

           call move(topt_out_s,data(1,1,i+9),nnxp,nnyp)
           call move(topt_out_ln,data(1,1,i+10),nnxp,nnyp)
           call move(topt_out_lt,data(1,1,i+11),nnxp,nnyp)
c          call move(topt_stag_out,data(1,1,i+12),nnxp,nnyp) !51

           call get_gridgen_var(ngrids,ngrids,var,comment)

	ELSE ! rotlat case

!!!!!	these doloops possibly redundant.  Defined above????

!        DO J=1,NNYP
!        DO I=1,NNXP
!              data(I,J,1)=hlat(I,J)
!              data(I,J,2)=hlon(I,J)
!              data(I,J,3)=vlat(I,J)
!              data(I,J,4)=vlon(I,J)
!        ENDDO
!        ENDDO

           data(:,:,18)=topt_out

c compute map/grid information. Reallocate GEODAT3D array
           allocate (GEODAT3D(nnxp,nnyp,4))

           call vecrot_rotlat(nnxp,nnyp,mdlat,mdlon
     +,vlat,vlon,GEODAT3D(:,:,1),GEODAT3D(:,:,2))

!        write(6,*) 'cos term'
        do J=nnyp,1,-nnyp/40
!        write(6,833) (GEODAT3D(I,J,1),I=1,nnxp,nnxp/15)
        enddo
!        write(6,*) 'sin term'
        do J=nnyp,1,-nnyp/40
!        write(6,833) (GEODAT3D(I,J,2),I=1,nnxp,nnxp/15)
        enddo

  833	format(30(f5.1,1x))

           i=10

!       sin, then cosine?
           data(:,:,i)  =GEODAT3D(:,:,2)
           data(:,:,i+1)=GEODAT3D(:,:,1)
c

!!!!    coriolis just dependent on lat (hopefully)
           call get_coriolis_components(nnxp,nnyp,vlat
     +,GEODAT3D)
           data(:,:,i+2)=GEODAT3D(:,:,1)
           data(:,:,i+3)=GEODAT3D(:,:,2)

           deallocate(GEODAT3D)

c          call move(static_albedo,data(1,1,i+4),nnxp,nnyp)

           call move(topt_out_s,data(1,1,i+5),nnxp,nnyp)
           call move(topt_out_ln,data(1,1,i+6),nnxp,nnyp)
           call move(topt_out_lt,data(1,1,i+7),nnxp,nnyp)
c          call move(topt_stag_out,data(1,1,i+8),nnxp,nnyp) !18

!030513  Made modifications within the get_gridgen_var code
!	to handle NMM variable set

        write(6,*) 'calling get_gridgen_var with ngrids= ', ngrids
           call get_gridgen_var(ngrids,ngrids,var,comment)

        ENDIF

        else ! non-wrfsi

           call move(lats(1,1,1),data(1,1,1),nnxp,nnyp)            ! KWD
           call move(lons(1,1,1),data(1,1,2),nnxp,nnyp)            ! KWD
           call move(topt_out,data(1,1,3),nnxp,nnyp)               ! KWD
           call move(topt_out_s,data(1,1,7),nnxp,nnyp)             ! JS
           call move(topt_out_ln,data(1,1,8),nnxp,nnyp)            ! JS
           call move(topt_out_lt,data(1,1,9),nnxp,nnyp)            ! JS

           call get_gridgen_var(ngrids,ngrids,var,comment)

        endif

!mp
        if (c6_maproj .ne. 'rotlat') then
          if(c10_grid_f(1:lf).eq.'wrfsi')then
             filename=c10_grid_f(1:lf)//'.d'//cnest//'.cdl'
          else
             filename = c10_grid_f(1:lf)//'.cdl'
          endif
        else
          filename = c10_grid_f(1:lf)//'.rotlat.cdl'
        endif

	write(6,*) 'using filename: ', filename
        call s_len(filename,lfn)
        call get_directory('cdl',cdl_dir,lcdl)

	write(6,*) 'file= ', cdl_dir(1:lcdl)//filename(1:lfn)
        INQUIRE(FILE=cdl_dir(1:lcdl)//filename(1:lfn),EXIST=exist)

        if(.not.exist)then
           print*,'Error: Could not find file '
     +           ,cdl_dir(1:lcdl)//filename(1:lfn)
           print*,'c10_grid_f: ',c10_grid_f(1:lf)
           istatus = 0
           return
        endif

        IF (c6_maproj .ne. 'rotlat') THEN

        call check_domain(lats(1,1,ns),lons(1,1,ns)
     +,nnxp,nnyp,grid_spacing_m,1,istat_chk)

        write(6,*)'deltax = ',deltax

        if(istat_chk .eq. 1)then
            write(6,*)'check_domain: status = ',istat_chk
        else
            write(6,*)'ERROR in check_domain: status = ',istat_chk
        endif

	ENDIF

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	write(6,*) 'call put_laps_static ', ngrids

	write(6,*) 'model= ', model

        c10_grid_fname = 'nest7grid'
        call s_len(c10_grid_fname,len_fname)

!       Do zin calc (note this is last [kmax] element in data array)
        if(c10_grid_fname(1:len_fname).eq.'nest7grid')then
          avgelem=3
          zinelem=ngrids
          write(6,*)' Calculating zin field in array element ',zinelem
          do i = 1,nnxp
          do j = 1,nnyp
              psa = ztopsa(data(i,j,avgelem)) ! This is the AVG data (3rd element)
              data(i,j,zinelem) = (20.0 - ((psa - 100.0) * 0.02))
          enddo ! j
          enddo ! i
        endif

        call put_laps_static(grid_spacing_m,model,comment,var
     1,data,nnxp,nnyp,ngrids,ngrids,std_lat,std_lat2,std_lon
     1,c6_maproj,deltax,deltay,mdlat,mdlon,lli,llj,iuri,iurj
     1,iparent_id,iratio_to_parent,nest,c10_grid_fname)

        istatus = istat_chk

        itstatus=ishow_timer()
        print*,' -------------------------------'
        print*,' Total Elapsed time: ', itstatus
        print*,' -------------------------------'

        deallocate (data)

	return
	End
c
c ================================================================
c
      SUBROUTINE GEODAT(n2,n3,erad,rlat,wlon1,xt,yt,deltax,deltay
     1 ,DATR,DATS,DATLN,DATLT,OFN,WVLN,SILWT,which_data,maxdatacat
     1 ,istat_files)

      include 'trigd.inc'

      implicit none

      integer n2,n3
      integer lbf
      integer lb,mof,np,niq,njq,nx,ny,isbego,iwbego,
     1  iblksizo,no,iodim,istat_files,maxdatacat

      real erad,rlat,wlon1,deltax,deltay,wvln,silwt

c wvln = TOPTWVL_PARM_WRF fror wrfsi.nl section &hgridspec
c siwt = SILAVWT_PARM_WRF      "

      real DATR(N2,N3)
      real DATS(N2,N3,maxdatacat)
      real DATLN(N2,N3)
      real DATLT(N2,N3)

c this parameter can be increased if we ever read data with
c resolution finer than 30 sec, or if the tilesize for 30s
c data becomes greater than 10x10 deg.

      PARAMETER(IODIM=5800000)

      real xt(N2),YT(N3)
      real deltallo,deltaxq,deltayq,
     1  deltaxp,deltayp

      real RSOFF,RWOFF

      real std_lon
      integer istatus
      integer lcat

      CHARACTER*(*) OFN
      character*180 TITLE
      logical which_data
C
c *********************
      nx=n2-1
      ny=n3-1
c ****************************

      lcat = 1
      LBF=INDEX(OFN,' ')-1
      TITLE=OFN(1:LBF)//'HEADER'

      LB=INDEX(TITLE,' ')-1

      CALL JCLGET(29,TITLE(1:LB),'FORMATTED',1,istat_files)
      if(istat_files .ne. 1)then
         write(6,*)' Warning in gridgen_model opening HEADER: '
     1,' geog path = ', title(1:lb)
         return
      endif

      READ(29,*)IBLKSIZO,NO,ISBEGO,IWBEGO,RSOFF,RWOFF
c 2    FORMAT(4I5,2(F10.8))
      print *,'title=',title
      print *,'rsoff,rwoff = ',rsoff,rwoff
      print *,'isbego,iwbego =',isbego,iwbego
      print *,'iblksizo,no =',iblksizo,no
      CLOSE(29)

      if(NO .gt. 1)then
         DELTALLO=FLOAT(IBLKSIZO)/FLOAT(NO-1)
      elseif(NO .eq. 1)then
         DELTALLO=FLOAT(IBLKSIZO)/FLOAT(NO)
      else
         print*,'HEADER value NO = 0'
         return
         istat_files = 0
      endif

      MOF=IODIM/(NO*NO)

      if(ofn(lbf:lbf).eq.'G'.or.ofn(lbf:lbf).eq.'A')then
         lcat=12
         if(no.eq.1250)then
            MOF=1
         endif
      endif
c SG97 MOF determines the number of files held in buffer while reading
c SG97 DEM data; it saves some time when buffer data can be used instead
c SG97 of reading DEM file again. Originally MOF was 4.
      if (MOF.gt.10) MOF=5

      DELTAXQ=0.5*WVLN*DELTAX
      DELTAYQ=0.5*WVLN*DELTAY
      print *,'deltaxq,deltayq=',deltaxq,deltayq
      NP=MIN(10,MAX(1,INT(DELTAXQ/(DELTALLO*111000.))))
      print *,' np=',np
      DELTAXP=DELTAXQ/FLOAT(NP)
      DELTAYP=DELTAYQ/FLOAT(NP)
      NIQ=INT(FLOAT(NX)*DELTAX/DELTAXQ)+4
      NJQ=INT(FLOAT(NY)*DELTAY/DELTAYQ)+4
C
      call get_standard_longitude(std_lon,istatus)
      if(istatus .ne. 1)then
          write(6,*)' Error calling laps routine'
          stop
      endif

      CALL SFCOPQR(NO,MOF,NP,NIQ,NJQ,N2,N3,lcat
     +,XT,YT,90.,std_lon,ERAD,RWOFF,RSOFF
     +,DELTALLO,DELTAXP,DELTAYP,DELTAXQ,DELTAYQ
     +,IBLKSIZO,ISBEGO,IWBEGO,DATR,DATS,DATLN,DATLT
     +,OFN,WVLN,SILWT,which_data,maxdatacat,istat_files)
      RETURN
      END

C
C     ******************************************************************
C
      SUBROUTINE SFCOPQR(NO,MOF,NP,NIQ,NJQ,N2,N3,lcat
     +          ,XT,YT,RLAT,WLON1,ERAD,RWOFF,RSOFF
     +          ,DELTALLO,DELTAXP,DELTAYP,DELTAXQ,DELTAYQ
     +          ,IBLKSIZO,ISBEGO,IWBEGO,DATR,DATS,DATLN,DATLT
     +          ,OFN,WVLN,SILWT,dem_data,maxdatacat,istat_files)

c JS: removed dato array from subroutine argument list
c JS: added RWOFF/RSOFF - West and South offset of tile data

      real,  allocatable ::  dato(:,:,:,:)    !dato(no,no,mof,lcat)

      real,  allocatable ::  DATP(:,:,:)
      real,  allocatable ::  DATQ(:,:)
      real,  allocatable ::  DATQS(:,:,:)
      real,  allocatable ::  DATSM(:,:)
      real,  allocatable ::  DATSMX(:,:)
      real,  allocatable ::  DATSLN(:,:)
      real,  allocatable ::  DATSLT(:,:)

      real,  intent(inout) ::  DATLN(N2,N3)
      real,  intent(out)   ::  DATLT(N2,N3)
      real,  intent(out)   ::  DATR(N2,N3)
      real,  intent(out)   ::  DATS(N2,N3,maxdatacat)

      real ISO(MOF),IWO(MOF),XT(N2),YT(N3),rlat,wlon1,
     +     erad,deltallo,deltaxp,deltayp,deltaxq,deltayq,
     +     wvln,silwt,xq,yq,xp,yp,xcentr,ycentr,glatp,               ! pla,plo,
     +     glonp,rio,rjo,wio2,wio1,wjo2,wjo1,xq1,yq1,
     +     rwoff,rsoff

      real r_missing_data
      real xr,yr,rval,sh,sha,rh,rha,rhn,rht,shn,sht
      real shln,shlt,rhln,rhlt
      real delta_ln(np,np),delta_lt(np,np)

c     real xpmn,xpmx,ypmn,ypmx
      real xp1,xp2,yp1,yp2
      real xpcentr,ypcentr

      real  pctcat(maxdatacat)

      integer lp
      integer ixr,iyr
      integer lent

      CHARACTER*180 OFN
      CHARACTER*180 TITLE3,TITLE3_last_read,TITLE3_last_inquire
      CHARACTER*3   TITLE1
      CHARACTER*4   TITLE2
      CHARACTER*10  cdatatype

      LOGICAL L1,L2,dem_data,l_string_contains

      data icnt/0/
      save icnt
C
      print *,'no,mof,np,niq,njq=',no,mof,np,niq,njq

      istat_files = 1

      NONO=NO*NO
      XCENTR=0.5*(XT(1)+XT(N2))
      YCENTR=0.5*(YT(1)+YT(N3))
      print *,xt(1),xt(n2),xcentr
      print *,'deltaxp=',deltaxp
      NOFR=0
      DO 11 IOF=1,MOF
         ISO(IOF)=0
         IWO(IOF)=0
  11  continue

      TITLE3_last_read    = '/dev/null'
      TITLE3_last_inquire = '/dev/null'

      lcat=1
      len=index(ofn,' ')
      if(ofn(len-1:len-1).eq.'V'.or.
     &   ofn(len-1:len-1).eq.'I')then
         icnt = 0
         cdatatype='landuse'
         if(ofn(len-1:len-1).eq.'I')cdatatype='islope'
      elseif(ofn(len-1:len-1).eq.'O')then
         icnt = 0
         cdatatype='soiltype'
      elseif(ofn(len-1:len-1).eq.'U' .or.
     &       ofn(len-1:len-1).eq.'H')then
         cdatatype='topography'
      endif

      print*,'SFCOPQR: cdatatype = ',cdatatype

      call s_len(cdatatype,lent)

      allocate(dato(no,no,mof,lcat))

      allocate (DATP(NP,NP,lcat),
     &          DATQ(NIQ,NJQ),
     &          DATSM(NIQ,NJQ),
     &          DATSMX(NIQ,NJQ),
     &          DATSLN(NIQ,NJQ),
     &          DATSLT(NIQ,NJQ),
     &          DATQS(NIQ,NJQ,maxdatacat))

      call get_r_missing_data(r_missing_data,istatus)
      if(istatus.ne.1)then
         print*,'failed to get r_missing_data'
         return
      endif

      DO 15 JQ=1,NJQ
         print *,'jq,njq,niq,nofr=',jq,njq,niq,nofr
         DO 16 IQ=1,NIQ

            XQ=(FLOAT(IQ)-0.5*FLOAT(NIQ+1))*DELTAXQ+XCENTR
            YQ=(FLOAT(JQ)-0.5*FLOAT(NJQ+1))*DELTAYQ+YCENTR

            xpmn=1.0e30
            ypmn=1.0e30
c           xpmx=-1.0e30
c           ypmx=-1.0e30

            DO 17 JP=1,NP
               DO 18 IP=1,NP

                  XP=XQ+(FLOAT(IP)-0.5*FLOAT(NP+1))*DELTAXP
                  YP=YQ+(FLOAT(JP)-0.5*FLOAT(NP+1))*DELTAYP

                  call xy_to_latlon(XP,YP,erad,GLATP,GLONP)

                  glatp = max(-89.9999,min(89.9999,glatp - rsoff))
                  glonp = glonp - rwoff
                  if(glonp.ge.180.) glonp = glonp - 360.
                  if(glonp.le.-180.) glonp = glonp + 360.

c                 print *,'rlat,wlon1=',rlat,wlon1

                  ISOC=(INT((GLATP-FLOAT(ISBEGO))/FLOAT(IBLKSIZO)
     &          +200.)-200)*IBLKSIZO+ISBEGO
            IWOC=(INT((GLONP-FLOAT(IWBEGO))/FLOAT(IBLKSIZO)
     &          +400.)-400)*IBLKSIZO+IWBEGO

                  DO 19 IOFR=1,NOFR
                     JOFR=IOFR
                     IF(ISO(IOFR).EQ.ISOC.AND.IWO(IOFR).EQ.IWOC)GO TO 10
 19                 continue
                  ISOCPT=ABS(ISOC)/10
                  ISOCPO=ABS(ISOC)-ISOCPT*10
                  IWOCPH=ABS(IWOC)/100
                  IWOCPT=(ABS(IWOC)-IWOCPH*100)/10
                  IWOCPO=ABS(IWOC)-IWOCPH*100-IWOCPT*10
                  IF(ISOC.GE.0)THEN
                     WRITE(TITLE1,'(2I1,A1)')ISOCPT,ISOCPO,'N'
                  ELSE
                     WRITE(TITLE1,'(2I1,A1)')ISOCPT,ISOCPO,'S'
                  ENDIF

                  IF(IWOC.GE.0
     1               .and. IWOC .ne. 180                    ! 1998 Steve Albers
     1                                      )THEN
                     WRITE(TITLE2,'(3I1,A1)')IWOCPH,IWOCPT,IWOCPO,'E'
                  ELSE
                     WRITE(TITLE2,'(3I1,A1)')IWOCPH,IWOCPT,IWOCPO,'W'
                  ENDIF

                  LB=INDEX(OFN,' ')-1
                  TITLE3=OFN(1:LB)//TITLE1//TITLE2
                  LB=INDEX(TITLE3,' ')-1

                  if(TITLE3 .ne. TITLE3_last_inquire)then
                     INQUIRE(FILE=TITLE3(1:LB),EXIST=L1,OPENED=L2)
                     TITLE3_last_inquire = TITLE3
                  endif

                  IF(.NOT.L1)THEN
                     iwrite = 0

                     if(icnt .le. 100)then ! Reduce the output
                         iwrite=1

                     elseif(icnt .le. 1000)then
                         if(icnt .eq. (icnt/100)*100)iwrite=1

                     elseif(icnt .le. 10000)then
                         if(icnt .eq. (icnt/1000)*1000)iwrite=1

                     elseif(icnt .le. 100000)then
                         if(icnt .eq. (icnt/10000)*10000)iwrite=1

                     else
                         if(icnt .eq. (icnt/100000)*100000)iwrite=1

                     endif

                     if(iwrite .eq. 1)then
                        if(l_string_contains(TITLE3(1:LB),
     1                                       'world_topo_30s',
     1                                       istatus)             )then
                           PRINT*, ' ERROR: ',TITLE3(1:LB)
     1                            ,' DOES NOT EXIST ',icnt

                        elseif(l_string_contains(TITLE3(1:LB),
     1                                           'topo_30s',
     1                                       istatus)             )then
                           PRINT*, ' topo_30s file ',TITLE3(1:LB)
     1                            ,' does not exist, using topo_10m '
     1                            ,icnt

                        else ! Generic warning message
                           PRINT*, ' WARNING: ',TITLE3(1:LB)
     1                            ,' DOES NOT EXIST ',icnt

                        endif

                     endif ! iwrite

                     icnt = icnt + 1

c initialize these arrays as they may have some garbage in them
c if we don't actually read in any data.
c
                     DATP(IP,JP,:) = 0.
                     DELTA_LN(IP,JP) = 0.
                     DELTA_LT(IP,JP) = 0.
                     istat_files = 0
                     GO TO 20

                  ENDIF

                  IF(NOFR.GE.MOF)THEN
                     DO 21 IOF=1,MOF
                        ISO(IOF)=0
                        IWO(IOF)=0
21                    continue
                     NOFR=0
                  ENDIF
                  NOFR=NOFR+1
                  JOFR=NOFR

!                 Read the tile
                  if(TITLE3 .ne. TITLE3_last_read)then
                    if( (ofn(len-1:len).eq.'U').and.(no.eq.1200).or.
     .                   no.eq.1201 )then
                         if(no.eq.1201)then
                            print*,'Reading ', title3(1:lb)
                            CALL READ_DEM(29,TITLE3(1:LB),no,no,4,4, ! topo_3s experimental
     .                              DATO(1,1,NOFR,1),istat)
                         else
                            CALL READ_DEM(29,TITLE3(1:LB),no,no,2,2, ! world topo_30s
     .                              DATO(1,1,NOFR,1),istat)
                         endif
                      dem_data=.true.
                    elseif( (ofn(len-1:len).eq.'O') )then      ! soiltype top and bot layer
                      CALL READ_DEM(29,TITLE3(1:LB),no,no,1,4,
     .                              DATO(1,1,NOFR,1),istat)
                      dem_data=.true.
                    elseif( (ofn(len-1:len).eq.'V') )then      ! world USGS 30s landuse
                      CALL READ_DEM(29,TITLE3(1:LB),no,no,1,4,
     .                              DATO(1,1,NOFR,1),istat)
                      dem_data=.true.
                    elseif((ofn(len-1:len).eq.'I'))then      ! only islope in this code section
                      CALL READ_DEM_G(29,TITLE3(1:LB),no,no,1,lcat
     .                     ,nofr, 1,4, DATO,istat)
                      dem_data=.true.
                    elseif( (ofn(len-1:len).eq.'T') )then      ! soiltemp - obsolete in this code
                      CALL READ_DEM(29,TITLE3(1:LB),no,no,2,2,
     .                              DATO(1,1,NOFR,1),istat)
                      dem_data=.true.
                    else                                       ! other
                      CALL JCLGET(29,TITLE3(1:LB),'FORMATTED',0,istatus)
                      CALL VFIREC(29,DATO(1,1,NOFR,1),NONO,'LIN')
                      if ((ofn(len-1:len).eq.'U').and.(no.eq.121)) then
                        dem_data=.false.                       ! topo_30s
                      endif
                    endif

                    if(istat.ne.0)then
                       print*,'Error returned: SFCOPQR: READ_DEM'
                       return
                    endif

                    TITLE3_last_read = TITLE3

c                   print *,'nofr,dato=',nofr,dato(1,1,nofr)
                    CLOSE(29)

                  else
                    write(6,*)' We have made the code more efficient'

                  endif ! Is this a new file we haven't read yet?

                  ISO(NOFR)=ISOC
                  IWO(NOFR)=IWOC
10		  continue

                  RIO=(GLONP-FLOAT(IWOC))/DELTALLO+1.
                  RJO=(GLATP-FLOAT(ISOC))/DELTALLO+1.

!                 Prevent Bounds Error (Steve Albers)
                  if(RIO .lt. 1.0)then
                      if(RIO .gt. 0.98)then
                          write(6,*)' Reset RIO for Machine Epsilon'
                          RIO = 1.0
                      elseif(RIO .lt. 0.5)then
                          write(6,*)' ERROR: RIO out of bounds',RIO
                          stop
                      endif
                  endif

                  if(RJO .lt. 1.0)then
                      if(RJO .gt. 0.98)then
                          write(6,*)' Reset RJO for Machine Epsilon'
                          write(6,*)JQ,IQ,
     1                          IP,JP,IO1,JO1,JOFR,RIO,RJO,GLATP,ISOC
                          RJO = 1.0
                      elseif(RJO .lt. 0.5)then
                          write(6,*)' ERROR: RJO out of bounds',RJO
                          write(6,*)JQ,IQ,
     1                          IP,JP,IO1,JO1,JOFR,RIO,RJO,GLATP,ISOC
                          stop
                      endif
                  endif

C Interp OK for continuous data such as topography

                  if(cdatatype.eq.'topography')then

                   IO1=INT(RIO)
                   JO1=INT(RJO)
                   IO2=IO1+1
                   JO2=JO1+1
                   WIO2=RIO-FLOAT(IO1)
                   WJO2=RJO-FLOAT(JO1)
                   WIO1=1.0-WIO2
                   WJO1=1.0-WJO2

                   do LP = 1,lcat

                   DATP(IP,JP,LP)=WIO1*(WJO1*DATO(IO1,JO1,JOFR,LP)
     +                                 +WJO2*DATO(IO1,JO2,JOFR,LP))
     +                           +WIO2*(WJO1*DATO(IO2,JO1,JOFR,LP)
     +                                 +WJO2*DATO(IO2,JO2,JOFR,LP))

!S & W-facing slopes > 0.
                   DELTA_LN(IP,JP)=
     .           ((DATO(IO2,JO1,JOFR,LP)-DATO(IO1,JO1,JOFR,LP))+
     .            (DATO(IO2,JO2,JOFR,LP)-DATO(IO1,JO2,JOFR,LP)))*.5

                   DELTA_LT(IP,JP)=
     .           ((DATO(IO1,JO2,JOFR,LP)-DATO(IO1,JO1,JOFR,LP))+
     .            (DATO(IO2,JO2,JOFR,LP)-DATO(IO2,JO1,JOFR,LP)))*.5

                   enddo !LP = 1,lcat

                  else

C Nearest grid point for landuse and soiltype

                   IO1=NINT(RIO)
                   JO1=NINT(RJO)
                   do LP = 1,lcat
                    DATP(IP,JP,LP)= DATO(IO1,JO1,JOFR,LP)
                   enddo

                  endif ! cdatatype eq topography.

20               CONTINUE
18             continue ! IP
17           continue ! JP

!           print*,'xpmx/xpmn//ypmx/ypmn/ ',xpmx,xpmn,ypmx,ypmn

! Calculate average and silhouette terrain, then apply SILWT weight

            if(cdatatype(1:lent).eq.'topography')then

             SHA=0.
             RHA=0.
             RHLN=0.
             RHLT=0.
             shmax=0.

             DO 22 JP=1,NP
               SH=0.
               RH=0.
               RHN=0.
               RHT=0.
               DO 23 IP=1,NP
!                 Test for missing - then go to 16?
                  SH=max(SH,DATP(IP,JP,1))
                  RH=RH+DATP(IP,JP,1)
                  RHN=RHN+DELTA_LN(IP,JP)
                  RHT=RHT+DELTA_LT(IP,JP)
23             continue ! IP
               SHA=SHA+SH/(2.*FLOAT(NP))
               RHA=RHA+RH
               RHLN=RHLN+RHN
               RHLT=RHLT+RHT
               SHMAX=max(SHMAX,SH)
22           continue ! JP

             RHA=RHA/FLOAT(NP*NP)
             RMS=0.0
             DO 24 IP=1,NP
               SH=0.
               DO 25 JP=1,NP
                  SH=max(SH,DATP(IP,JP,1))
                  RMS=RMS+((DATP(IP,JP,1)-RHA)*(DATP(IP,JP,1)-RHA))
25             continue ! JP
               SHA=SHA+SH/(2.*FLOAT(NP))
24           continue ! IP

             DATQS(IQ,JQ,1)=SQRT(RMS/FLOAT(NP*NP))
             DATQ(IQ,JQ)=SHA*SILWT+RHA*(1.-SILWT)
             DATSM(IQ,JQ)=RHA                           !mean value of points used for IQ,JQ
             DATSMX(IQ,JQ)=SHMAX                        !max value from points used for IQ,JQ
             DATSLN(IQ,JQ)=RHLN/FLOAT(NP*NP)/DELTAXP
             DATSLT(IQ,JQ)=RHLT/FLOAT(NP*NP)/DELTAYP

c            print *,'datq=',datq(iq,jq)

            elseif(cdatatype(1:lent).eq.'islope'    .or.
     &             cdatatype(1:lent).eq.'landuse'   .or.
     &             cdatatype(1:lent).eq.'soiltype'    )then

             call compute_categories(cdatatype,np*np,DATP(1,1,1)
     &               ,maxdatacat,domcat,pctcat)
             datq(iq,jq)=domcat
             datqs(iq,jq,:)=pctcat(:)

            elseif(cdatatype(1:lent).eq.'greenfrac'   )then

c dominant greenness fraction for each month

             do lp=1,lcat
              call compute_categories(cdatatype,np*np,DATP(1,1,lp)
     &                               ,1,domcat,pctcat)
              datqs(iq,jq,lp)=domcat
             enddo

            endif

16       continue ! IQ
15    continue ! JQ

      print *,'after 15'

      XQ1=(1.-0.5*FLOAT(NIQ+1))*DELTAXQ+XCENTR
      YQ1=(1.-0.5*FLOAT(NJQ+1))*DELTAYQ+YCENTR

      if(cdatatype(1:lent).eq.'topography')then

        print*
        print*,'Before GDTOST2'
        print*,'--------------'
        print*,'datq(1,1)/(niq,njq)= ',datq(1,1),datq(niq,njq)
        print*,'datqs(1,1,1)/(niq,njq)= ',datqs(1,1,1),datqs(niq,njq,1)
        print*,'datsln(1,1)/(niq,njq)= ',datsln(1,1),datsln(niq,njq)
        print*,'datslt(1,1)/(niq,njq)= ',datslt(1,1),datslt(niq,njq)
        print*,'Mean/Max topo at IQ,JQ (1,1)/(niq,njq): '
     +,datsm(1,1),datsmx(1,1),datsm(niq,njq),datsmx(niq,njq)

        DO 28 JR=1,N3
         DO 29 IR=1,N2
           XR=(XT(IR)-XQ1)/DELTAXQ+1.
           YR=(YT(JR)-YQ1)/DELTAYQ+1.

           CALL GDTOST2(DATQ,NIQ,NJQ,XR,YR,RVAL)
           DATR(IR,JR)=max(0.,RVAL)
           if( DATR(IR,JR).gt.30000. )then
               print*,'Warning: value out of bounds'
           endif

           CALL GDTOST2(DATQS,NIQ,NJQ,XR,YR,RVAL)
           DATS(IR,JR,1)=max(0.,RVAL)
           CALL GDTOST2(DATSLN,NIQ,NJQ,XR,YR,RVAL)
           DATLN(IR,JR)=RVAL
           CALL GDTOST2(DATSLT,NIQ,NJQ,XR,YR,RVAL)
           DATLT(IR,JR)=RVAL

 29      CONTINUE
 28     CONTINUE

        print*,'After GDTOST2'
        print*,'-------------'
        print*,'datr(1,1)/(n2,n3)= ',datr(1,1),datr(N2,N3)
        print*,'dats(1,1,1)/(n2,n3)= ',dats(1,1,1),dats(n2,n3,1)
        print*,'datln(1,1)/(n2,n3)= ',datln(1,1),datln(n2,n3)
        print*,'datlt(1,1)/(n2,n3)= ',datlt(1,1),datlt(n2,n3)

      elseif(cdatatype(1:lent).eq.'landuse'.or.
     +       cdatatype(1:lent).eq.'islope' .or.
     +       cdatatype(1:lent).eq.'soiltype')then

        DO 38 JR=1,N3
         DO 39 IR=1,N2
            IXR=NINT((XT(IR)-XQ1)/DELTAXQ)+1.
            IYR=NINT((YT(JR)-YQ1)/DELTAYQ)+1.
            if(ixr.lt.1)ixr=1
            if(iyr.lt.1)iyr=1
            if(ixr.gt.n2)ixr=niq
            if(iyr.gt.n3)iyr=njq

            datr(ir,jr)=  datq(ixr,iyr)     !dominant category
            dats(ir,jr,:)=datqs(ixr,iyr,:)  !percent dist for ea category

 39      CONTINUE
 38     CONTINUE

      endif

      deallocate(dato)
      deallocate(DATP,
     &           DATQ,
     &           DATQS,
     &           DATSM,
     &           DATSMX,
     &           DATSLN,
     &           DATSLT)

      RETURN
      END

      subroutine vfirec(iunit,a,n,type)
      character*1 vc
      character*(*) type
      common/vform/vc(0:63)
      character line*80, cs*1
      dimension a(*)

      if(vc(0).ne.'0') call vfinit

      ich0=ichar('0')
      ich9=ichar('9')
      ichcz=ichar('Z')
      ichlz=ichar('z')
      ichca=ichar('A')
      ichla=ichar('a')

      read(iunit,10)nn,nbits,bias,fact
 10   format(2i8,2e20.10)
      if(nn.ne.n) then
         print*,' Word count mismatch on vfirec record '
         print*,' Words on record - ',nn
         print*,' Words expected  - ',n
         stop 'vfirec'
      endif

      nvalline=(78*6)/nbits
      nchs=nbits/6
      do 20 i=1,n,nvalline
         read(iunit,'(a78)') line
         ic=0
         do 30 ii=i,i+nvalline-1
            isval=0
            if(ii.gt.n) go to 20
            do 40 iii=1,nchs
               ic=ic+1
               cs=line(ic:ic)
               ics=ichar(cs)
               if(ics.le.ich9)then
                  nc=ics-ich0
               elseif(ics.le.ichcz) then
                  nc=ics-ichca+10
               else
                  nc=ics-ichla+36
               endif
               isval=ior(ishft(nc,6*(nchs-iii)),isval)
 40         continue
            a(ii)=isval
 30      continue
 20   continue

      facti=1./fact
      if(type.eq.'LIN') then
         do 48 i=1,n
            a(i)=a(i)*facti-bias
 48      continue
      elseif(type.eq.'LOG') then
         scfct=2.**(nbits-1)
         do 55 i=1,n
            a(i)=sign(1.,a(i)-scfct)
     +           *(10.**(abs(20.*(a(i)/scfct-1.))-10.))
 55      continue
      endif

      return
      end
c
cc ------------------------------------------------------------------
c
      subroutine vfinit
      character*1vc,vcscr(0:63)
      common/vform/vc(0:63)
      data vcscr/'0','1','2','3','4','5','6','7','8','9'
     +,'A','B','C','D','E','F','G','H','I','J'
     +,'K','L','M','N','O','P','Q','R','S','T'
     +,'U','V','W','X','Y','Z','a','b','c','d'
     +,'e','f','g','h','i','j','k','l','m','n'
     +,'o','p','q','r','s','t','u','v','w','x'
     +,'y','z','{','|'/

      do10n=0,63
      vc(n)=vcscr(n)
  10  continue

      return
      end
C +------------------------------------------------------------------+
      FUNCTION INTLSHFT(IWORD,NSHFT)
C
C       This function shifts IWORD to the left NSHFT bits in a
C         circular manner.
C
      INTLSHFT=ISHFT(IWORD,NSHFT)
      RETURN
      END
C +------------------------------------------------------------------+
      FUNCTION INTOR(IWORD1,IWORD2)
C
C       This function performs a bit-by-bit OR between IWORD1 and
C         IWORD2.
C
      INTOR=IOR(IWORD1,IWORD2)
      RETURN
      END

      SUBROUTINE BINOM2(X1,X2,X3,X4,Y1,Y2,Y3,Y4,XXX,YYY)
      implicit none
      real x1,x2,x3,x4,y1,y2,y3,y4,xxx,yyy,
     +   wt1,wt2,yz22,yz23,yz24,yz11,yz12,yz13,yoo
      integer istend
c      COMMON/BIN/ITYPP,I0X,I1X,I2X,YOO
       YYY=1E30
       IF(X2.GT.1.E19.OR.X3.GT.1.E19.OR.
     +   Y2.GT.1.E19.OR.Y3.GT.1.E19)RETURN
      WT1=(XXX-X3)/(X2-X3)
      WT2=1.0-WT1
      ISTEND=0
      IF(Y4.LT.1.E19.AND.X4.LT.1.E19) GO TO 410
      YZ22=WT1
      YZ23=WT2
      YZ24=0.0
      ISTEND= 1
410   IF(Y1.LT.1.E19.AND.X1.LT.1.E19) GO TO 430
      YZ11=0.0
      YZ12=WT1
      YZ13=WT2
      IF(ISTEND.EQ.1)GO TO 480
      GO TO 450
430   YZ11=(XXX-X2)*(XXX-X3)/((X1-X2)*(X1-X3))
      YZ12=(XXX-X1)*(XXX-X3)/((X2-X1)*(X2-X3))
      YZ13=(XXX-X1)*(XXX-X2)/((X3-X1)*(X3-X2))
      IF(ISTEND.EQ.  1    ) GO TO 470
450   YZ22=(XXX-X3)*(XXX-X4)/((X2-X3)*(X2-X4))
      YZ23=(XXX-X2)*(XXX-X4)/((X3-X2)*(X3-X4))
      YZ24=(XXX-X2)*(XXX-X3)/((X4-X2)*(X4-X3))
470   YYY=WT1*(YZ11*Y1+YZ12*Y2+YZ13*Y3)+WT2*(YZ22*Y2+YZ23*Y3+YZ24*Y4)
       GO TO 490
480      YYY=WT1*Y2+WT2*Y3
490   YOO=YYY
      RETURN
      END
c
c determine dominant category 1-05-01 JS
c

      subroutine compute_categories(ctype,nnp,data,nlcat,domcat
     +,pctcat)

      implicit none

      character*(*) ctype

      integer nnp
      integer igc
      integer i,j,k
      integer nlcat
      integer maxcat
      integer lcat(nlcat)
      real    pctcat(nlcat)
      real    data(nnp)
      real    domcat
      real    sum_g

c categorical data types
      if(ctype.eq.'landuse'   .or.
     +   ctype.eq.'soiltype'  .or.
     +   ctype.eq.'islope'  )then

         do k=1,nlcat
            lcat(k)=0
            pctcat(k)=0.
         enddo

         do i=1,nnp
         do k=1,nlcat
            if(nint(data(i)).eq.k)then
               lcat(k)=lcat(k)+1
            endif
         enddo
         enddo

         maxcat=-1
         do k=1,nlcat
            pctcat(k)=lcat(k)/float(nnp)
            if(lcat(k).gt.maxcat)then
               maxcat=lcat(k)
               domcat=float(k)
            endif
         enddo
         if(ctype.eq.'landuse')then
          if(pctcat(16).ge.0.5) then                            !!JBresch
	    domcat = 16                                         !!JBresch
          else if(domcat.eq.16.and.pctcat(16).lt.0.5)then       !!JBresch
            maxcat=-1
            do k=1,nlcat
               if(k.ne.16)then
                  if(lcat(k).gt.maxcat)then
                     maxcat=lcat(k)
                     domcat=float(k)
                  endif
               endif
            enddo
          endif
         endif

c quantitative data types
      elseif(ctype.eq.'greenfrac'.or.ctype.eq.'soiltemp')then

         sum_g=0.
         igc=0
         do i=1,nnp
          if(data(i).gt.0.0)then
             sum_g=sum_g+data(i)
             igc = igc+1
          endif
         enddo
         if(igc.gt.0)then
            domcat=sum_g/float(igc)
         else
            domcat=0.0
         endif

      endif

      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	subroutine zero2d(adum2d,nnxp,nnyp)

	integer I,J,NNXP,NNYP
	real adum2d(nnxp,nnyp)

	do J=1,nnyp
	do I=1,nnxp
	adum2d(I,J)=0.
	enddo
	enddo

	end subroutine zero2d

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	subroutine dumtodata(adum2d,nnxp,nnyp,IVAL,data)

	INTEGER NNXP,NNYP,IVAL
	real ADUM2D(:,:),data(:,:,:)

	do J=1,nnyp
	do I=1,nnxp
	data(I,J,IVAL)=adum2d(I,J)
	enddo
	enddo

	end subroutine dumtodata

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine smdhld(ime,jme,h,s,lines,nsmud)
!      parameter(ime=im+2,jme=jm+4)
      dimension ihw(jme),ihe(jme)
      dimension h(ime,jme),s(ime,jme)
     &         ,hbms(ime,jme),hne(ime,jme),hse(ime,jme)
!-----------------------------------------------------------------------
          do j=1,jme
      ihw(j)=-mod(j,2)
      ihe(j)=ihw(j)+1
          enddo
!-----------------------------------------------------------------------

              do j=1,jme
          do i=1,ime
      hbms(i,j)=1.-s(i,j)
          enddo
              enddo
!
      jmelin=jme-lines+1
      ibas=lines/2
      m2l=mod(lines,2)
!
              do j=lines,jmelin
          ihl=ibas+mod(j,2)+m2l*mod(j+1,2)
          ihh=ime-ibas-m2l*mod(j+1,2)
!
          do i=ihl,ihh
      hbms(i,j)=0.
          enddo
              enddo
!-----------------------------------------------------------------------
                  do ks=1,nsmud
!-----------------------------------------------------------------------
              do j=1,jme-1
          do i=1,ime-1
      hne(i,j)=h(i+ihe(j),j+1)-h(i,j)
          enddo
              enddo
              do j=2,jme
          do i=1,ime-1
      hse(i,j)=h(i+ihe(j),j-1)-h(i,j)
          enddo
              enddo
!
              do j=2,jme-1
          do i=1+mod(j,2),ime-1
      h(i,j)=(hne(i,j)-hne(i+ihw(j),j-1)
     &       +hse(i,j)-hse(i+ihw(j),j+1))*hbms(i,j)*0.125+h(i,j)
          enddo
              enddo
!-----------------------------------------------------------------------
                      enddo
!-----------------------------------------------------------------------
      return
      end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense