dir=  grid
get_tile_list.f gridgen_model.f gridgen_utils.f proc_geodat.f put_laps_static.f rd_static_attr_sub.f
COMPARE get_tile_list.f
14a15,18
>       integer   imxlat
>       integer   imnlat
>       integer   imxlon
>       integer   imnlon
19d22
<       integer   iwocdif
28c31,32
<       character*(*) ctile_name_list(maxtiles)
---
> !mp      character*(*) ctile_name_list(maxtiles)
>       character*8 ctile_name_list(maxtiles)
33a38,40
> 	write(6,*) 'min_lat, max_lat: ', min_lat, max_lat
> 	write(6,*) 'min_lon, max_lon: ', min_lon, max_lon
> 
45,54c52,54
<       if(IWOC1.lt.-180)then
<          if(itilesize.eq.180)then
<             iwocdif=iwoc2-iwoc1
<          else
<             iwocdif=iwoc2-iwoc1-itilesize
<          endif
<          IWOC1=IWOC1+IWOCDIF
<          IWOC2=IWOC2+abs(IWOCDIF)
<       endif
<       print*,'Noddy IWOC1, IWOC2 ',IWOC1,IWOC2
---
> 	write(6,*) 'ISOC1,IWOC1,ISOC2,IWOC2: ', ISOC1,IWOC1,ISOC2,IWOC2
> 
>       print*,'Noddy IWOC1, IWOC2 ',IWOC1,IWOC2, itilesize
64c64
<          IF(IWOC.GE.0.and.IWOC.LT.180) THEN
---
>          IF(IWOC.GE.0) THEN
69a70
> 	write(6,*) 'ewtitle= ', ewtitle
85a87,88
> 	write(6,*) 'nstitle= ', nstitle
> 
92,94c95,96
<       if(num_tiles_needed.eq.3.and.maxtiles.eq.2)then
<          num_tiles_needed = 2
<       endif
---
> 	write(6,*) 'num_tiles_needed,maxtiles: ', num_tiles_needed,
>      &			maxtiles
98a101,104
> 	write(6,*) 'ctilenamelist(itile_ns): ', 
>      & ctilenamelist(itile_ns)
> 	write(6,*) 'ctile_name_list(itile_ns): ', 
>      &			ctile_name_list(itile_ns)
done get_tile_list.f ======================
COMPARE gridgen_model.f
64,70d63
<         character cstaticdir*200
<         logical   localize
<         real,     allocatable  ::  moad_lat(:,:)
<         real,     allocatable  ::  moad_lon(:,:)
<         real,     allocatable  ::  dum2d(:,:)
<         real      La1,Lo1,LoV,Latin1,Latin2
<         real      Dx,Dy
72,73d64
<         include 'grid_fname.cmn'
< 
77,81c68,72
<         call get_config(istatus)
<         if(istatus.ne.1)then
<            print*,'Error returned: get_config'
<            stop
<         endif
---
> 	call get_grid_dim_xy(NX_L,NY_L,istatus)
> 	if (istatus .ne. 1) then
>            write (6,*) 'Error getting horizontal domain dimensions'
> 	   goto999
> 	endif
87a79
> 
90,97c82,86
<         if(c10_grid_fname(1:lf).eq.'wrfsi')then
<            call get_num_domains(num_domains,istatus)
<            if(istatus.ne.1)then
<               print*,'Error returned: get_num_domains'
<               stop
<            endif
<         else
<            num_domains=1
---
>         call get_n_staggers(n_staggers,istatus)
>         if(istatus.ne.1)then
>            print*,'Error getting n_staggers for',
>      &': ',c10_grid_fname(1:lf)
>            stop
100c89,93
<         do nest=1,num_domains
---
>         call get_stagger_index(istag,istatus)
>         if(istatus.ne.1)then
>            print*,'Error getting num_staggers_wrf from wrfsi.nl'
>            stop
>         endif
102c95
< c first detemine if this domain requires localizing or re-localizing
---
>         call Gridmap_sub(NX_L,NY_L,n_staggers,istag,istatus)
104c97,98
<            call get_directory('static',cstaticdir,len_dir)
---
>  999	write(6,*)' gridgen_model finish: istatus = ',istatus
>         write(6,*)
106,191d99
<            call rd_static_attr(cstaticdir,nest
<      .,c10_grid_fname,Nx_l,Ny_l,Dx,Dy,La1,Lo1,Latin1
<      .,Latin2,LoV,istatus)
< 
<            call eval_localization(cstaticdir,nest,localize
<      .,c10_grid_fname,La1,Lo1,istatus)
< 
<            if(localize)then
< 
<               call get_grid_dim_xy(NX_L,NY_L,istatus)
<               if (istatus .ne. 1) then
<                  print*, 'Error getting horizontal domain dimensions'
<                  goto999
<               endif
< 
<               call get_grid_spacing(Dx,istatus)
<               call get_domain_origin_parent(i_orig,j_orig,istatus)
< 
<               print*,'grid_spacing(m)  = ',Dx
<               print*,'grid dims        = ',nx_l,ny_l
<               print*,'grid origin i/j  = ',i_orig,j_orig
<               Dy=Dx
< 
<               if(nest.gt.1)then
<                  print*,' calculate center lat/lon '
<                  La1=moad_lat(i_orig,j_orig)
<                  Lo1=moad_lon(i_orig,j_orig)
<                  print*,' SW lat/lon = ',La1,Lo1
<                  call get_earth_radius(erad,istatus)
< c Get X/Y for grid center
<                  CALL POLAR_GP(La1,Lo1,Xcen,Ycen,Dx,Dy,NX_L,NY_L,-1.)
<                  call xy_to_latlon(Xcen,Ycen,erad,rlat,rlon)
<               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_fname(1:lf)
<                  stop
<               endif
< 
<               call Gridmap_sub(NX_L,NY_L,n_staggers,istag,nest,istatus)
< 
<               if(istatus.eq.1)then
<                  print*,' Gridmap_sub finish successfully'
<                  if(nest.eq.1)then
<                     allocate (moad_lat(nx_l,ny_l),
<      .                        moad_lon(nx_l,ny_l),
<      .                        dum2d(nx_l,ny_l))
<                     call get_laps_domain(nx_l,ny_l,c10_grid_fname
<      +,moad_lat,moad_lon,dum2d,istatus)
<                     if(istatus.lt.1)then
<                        print *,'Error reading MOAD lat, lon data.'
<                        stop
<                     endif
<                  endif
<               endif
< 
<               print*,' Gridmap_sub finished '
<               print*
< 
<            else
<               print*
<               print*,' No need to re-localize this domain '
<               print*
<               if(nest.eq.1)then
<                  allocate (moad_lat(nx_l,ny_l),
<      .                     moad_lon(nx_l,ny_l),
<      .                        dum2d(nx_l,ny_l))
<                  call get_laps_domain(nx_l,ny_l,c10_grid_fname
<      +,moad_lat,moad_lon,dum2d,istatus)
<                  if(istatus.lt.1)then
<                     print *,'Error reading MOAD lat, lon data.'
<                     stop
<                  endif
<               endif
<            endif
< 
<         enddo
< 
<  999	print*,' gridgen_model finish: istatus = ',istatus
<         print*
< 
<         stop
193,196d100
< c
< c --------------------------------------------------------------------
< c
<         subroutine Gridmap_sub(nnxp,nnyp,n_staggers,istag,nest,istatus)
197a102,105
> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
>        
>         subroutine Gridmap_sub(nnxp,nnyp,n_staggers,istag,istatus)
> 
201d108
<         logical lforce_ter_2_zero
206d112
<         integer nest
213a120
>         real  adum(nnxp,nnyp)
215,216c122,123
<         real, allocatable ::  topt_10   (:,:)
<         real, allocatable ::  topt_10_s (:,:)
---
>         real, allocatable ::  topt_10(:,:)
>         real, allocatable ::  topt_10_s(:,:)
219,220c126,127
<         real, allocatable ::  topt_30   (:,:)
<         real, allocatable ::  topt_30_s (:,:)
---
>         real, allocatable ::  topt_30(:,:)
>         real, allocatable ::  topt_30_s(:,:)
224,225c131,132
<         real  topt_out   (nnxp,nnyp)
<         real  topt_out_s (nnxp,nnyp)
---
>         real  topt_out(nnxp,nnyp)
>         real  topt_out_s(nnxp,nnyp)
238,240d144
<         real, allocatable ::  adum2d  (:,:)
<         real, allocatable ::  adum3d  (:,:,:)
< 
243a148,151
> !mp
> !	real test_hold(nnxp,nnyp)
> !mp
> 
246a155,156
> 	real, allocatable, dimension(:,:):: hlat,hlon,vlat,vlon
> 
263d172
<         character*200 path_to_islope
271c180
<         character*1   cdatatype
---
>         character*1   cdatatype,type
278c187,188
<  
---
> 
> 
281,282c191,192
<      &,istatgrn,istattmp,istatslp,lat,topt_out,path_to_soiltemp
<      &,landmask,soiltemp_1deg,greenfrac,islope,istatus)
---
>      &,istatgrn,istattmp,lat,topt_out,path_to_soiltemp
>      &,landmask,soiltemp_1deg,greenfrac,istatus)
288d197
<           integer istatslp
295d203
<           real    islope(nnxp,nnyp)
309a218,226
> !mp
>         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
> 
313d229
< c add 1 more for islope JS: 2-25-03
315,316c231,243
<            ngrids=112     !97+12   ! 2d grids (including %dist for landuse
<                                    ! and two soiltype categories, green frac, albedo and others).
---
> 
> !mp
> !mp	need to tweak this number for rotlat
> !mp
> 
> 	IF (c6_maproj .eq. 'rotlat') THEN
> 	write(6,*) 'ngrids set to 100'
> 	ngrids=100
> 	ELSE !normal
>            ngrids=97+12 ! 2d grids (including %dist for landuse
>                         ! and two soiltype categories, green frac and albedo).
> 	ENDIF
> 
318c245
<            ngrids=38      !26+12   doesn't include terrain slope index and other stuff.
---
>            ngrids=26+12
321,328c248,252
<         if(.not. allocated (data))then
<            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
---
>         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
343a268,269
> c   the 10m pctl covers the world
> cc        ipctlfn=static_dir(1:len)// 'model/land_10m/L'
356a283,288
>         call get_path_to_pctl_10m(path_to_pctl10m,istatus)
>         if(istatus .ne. 1)then
>            write(6,*) 'Error getting path_to_pctl10m'
>            return
>         endif
> 
402,407d333
<         call get_path_to_islope(path_to_islope,istatus)
<         if(istatus .ne. 1)then
<            print*, 'Error getting path_to_islope'
<            return
<         endif
< 
417a344,347
>         call s_len(path_to_pctl10m,len)
>         print*,'path to pctl_10m:       ',path_to_pctl10m(1:len)
> 	path_to_pctl10m(len+1:len+2)='/L'
> 
446,449d375
<         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'
< 
568a495
> 	IF (c6_maproj .ne. 'rotlat') THEN
575a503,536
> 	ELSE ! rotlat case
> 
>         fact=2.*3.14159265*erad/360.
>         dphd=grid_spacing_m*(1/( (2.04)**0.5 * fact))
>         dlmd=dphd*1.02
> 
>         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,*) 'hlat(1,1),hlon(1,1): ', hlat(1,1),hlon(1,1)
>         write(6,*) 'vlat(1,1),vlon(1,1): ', vlat(1,1),vlon(1,1)
> 
> 
> 	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
> 
> 	ENDIF
> 
> 
> 
> 	
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
608a570
> 	ENDIF
609a572
> 
639,642d601
< c ------------ Terrain USGS 30 sec ----------------
< c ------------ and 10m if necessary ---------------
< c type = U
< c
643a603
>        print*,'Elapsed time: ', itstatus
648,650c608,610
<      +           topt_30_s(nnxp,nnyp),
<      +           topt_30_ln(nnxp,nnyp),
<      +           topt_30_lt(nnxp,nnyp))
---
>      +              topt_30_s(nnxp,nnyp),
>      +              topt_30_ln(nnxp,nnyp),
>      +              topt_30_lt(nnxp,nnyp))
651a612,634
> 	IF (c6_maproj .eq. 'rotlat') THEN
> 
> 	type='U'
> 
> 	CALL alt_10by10(nnxp,nnyp,1200,grid_spacing_m/1000.,
>      &	hlat,hlon,type,PATH_TO_TOPT30S,topt_30,topt_30_lt,
>      &	topt_30_ln,topt_30_s)
> 
> 
> 	do J=nnyp,1,-nnyp/30
> 	write(6,633) (topt_30(I,J),I=1,nnxp,nnxp/15)
> 	enddo
> 
>   633	format(20(f5.0,1x))
> 
> 
> !mp
> 	istatus_30s=1
> !mp
> 
> 
> 	ELSE
> 
656a640,641
> 	ENDIF
> 
681c666
<              print*,'Error: Aborting --- No static file created'
---
>              print*,'Aborting --- No static file created'
705a691,692
> 	write(6,*) 'all tiles avail...defining topt_out here: '
> 
711c698
<           
---
> 
720,731d706
< c
< c -------------- 30s Landuse -------------------------
< c type = V
< c
<        itstatus=ishow_timer()
<        print*
<        print*,' Calling GEODAT: Processing 30s landuse data.'
<        print*,' Re-allocate GEODAT3D ',nnxp,nnyp,' 24'
<        print*
<        allocate  (adum2d(nnxp,nnyp))
<        allocate  (GEODAT2D(nnxp,nnyp))
<        allocate  (GEODAT3D(nnxp,nnyp,24))
733,779d707
<        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
< 
<        ilndmsk=12  !land mask is 12 for both wrfsi and laps
<        if(c10_grid_fname(1:lf).eq.'wrfsi')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 39 are percent distributions
<           i=14
<           do j=1,24
<              data(:,:,i+j)=GEODAT3D(:,:,j)
<           enddo
<        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
< 
< 
784,790c712
< 
<           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
---
> 	IF (c6_maproj .ne. 'rotlat') THEN
793c715,716
<              data(i,j,indxter) = r_missing_data
---
>              data(i,j,51) = r_missing_data
> c            topt_stag_out(i,j)=r_missing_data
795c718
<              data(i-1,j-1,indxter)=result
---
>              data(i-1,j-1,51)=result
799,800c722,723
<           where((data(:,:,indxter).lt. 0.01).and.
<      1          (data(:,:,indxter).gt.-0.01))data(:,:,indxter)=0.0
---
>           where((data(:,:,51).lt. 0.01).and.
>      1          (data(:,:,51).gt.-0.01))data(:,:,51)=0.0
802c725
<        else  !must be LAPS
---
> 	ELSE
804,808c727,728
<           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
---
> !	for rotlat, topt_out is already on mass point stagger.
> !	No bilinear interp needed
809a730,740
> 	do J=1,nnyp	
> 	do I=1,nnxp
> 	data(i,j,18)=topt_out(i,j)
> 	enddo
> 	enddo
> 
>           where((data(:,:,18).lt. 0.01).and.
>      1          (data(:,:,18).gt.-0.01))data(:,:,18)=0.0
> 
> 	ENDIF
> 
811a743,747
>        write(6,*)
>        print *,'topt_out   =',topt_out(1,1),topt_out(nnxp,nnyp)
>        print *,'# of grid pts using 30 sec data =  ',icount_30
>        print *,'# of grid pts using 10 min data =  ',icount_10
>        print *,'# of grid pts using blended data = ',icount_ramp
813,817d748
<        print *,'terrain  = ',data(1,1,indxter),data(nnxp,nnyp,indxter)
<        print *,'# of grid pts using 30 sec terrain data =  ',icount_30
<        print *,'# of grid pts using 10 min terrain data =  ',icount_10
<        print *,'# of grid pts using blended terrain data = ',icount_ramp
<        print*
819,820c750,751
<        call s_len(static_dir,len)
<        open(10,file=static_dir(1:len)//'latlon.dat'
---
>         call s_len(static_dir,len)
>         open(10,file=static_dir(1:len)//'latlon.dat'
822c753
<        open(11,file=static_dir(1:len)//'topo.dat'
---
>         open(11,file=static_dir(1:len)//'topo.dat'
824c755
<        open(15,file=static_dir(1:len)//'corners.dat'
---
>         open(15,file=static_dir(1:len)//'corners.dat'
826,834c757,775
<        write(10)lats(1:nnxp,1:nnyp,1),lons(1:nnxp,1:nnyp,1)
<        close(10)
<        write(11)data(:,:,indxter)
<        close(11)
<        write(15,*)lats(1,1,ns),lons(1,1,ns)
<        write(15,*)lats(1,nnyp,ns),lons(1,nnyp,ns)
<        write(15,*)lats(nnxp,1,ns),lons(nnxp,1,ns)
<        write(15,*)lats(nnxp,nnyp,ns),lons(nnxp,nnyp,ns)
<        close(15)
---
> 	IF (c6_maproj .ne. 'rotlat') THEN
>         write(10)lats(1:nnxp,1:nnyp,1),lons(1:nnxp,1:nnyp,1)
> 	ELSE
>         write(10)hlat(1:nnxp,1:nnyp),hlon(1:nnxp,1:nnyp)
> 	ENDIF
>         close(10)
>         write(11)topt_out
>         close(11)
> 	IF (c6_maproj .ne. 'rotlat') THEN
>         write(15,*)lats(1,1,1),lons(1,1,1)
>         write(15,*)lats(1,nnyp,1),lons(1,nnyp,1)
>         write(15,*)lats(nnxp,1,1),lons(nnxp,1,1)
>         write(15,*)lats(nnxp,nnyp,1),lons(nnxp,nnyp,1)
> 	ELSE
>         write(15,*)hlat(1,1),hlon(1,1)
>         write(15,*)hlat(1,nnyp),hlon(1,nnyp)
>         write(15,*)hlat(nnxp,1),hlon(nnxp,1)
>         write(15,*)hlat(nnxp,nnyp),hlon(nnxp,nnyp)
> 	ENDIF
835a777,779
> c        write(16,*)topt_pctlfn
>         close(15)
> 
846c790
< 	    write(666,*) data(i,j,indxter)
---
> 	    write(666,*) topt_out(i,j)
851,854c795,804
< c
< c ------ Land Fraction from 30s Land Use ----------------------------
< c type = L (no longer used ... no land fraction data base)
< c
---
> c -----------------------------------------------------------
> 
> 
>         itstatus=ishow_timer()
>         print*,'Elapsed time: ', itstatus
>         allocate  (GEODAT2D(nnxp,nnyp))
>         allocate  (GEODAT3D(nnxp,nnyp,24))
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
855a806,882
>         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
>      +,adum,adum,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
> 
> 	call alt_categorical(nnxp,nnyp,hlat,hlon,PATH_TO_LUSE_30S,
>      &			     type,24,grid_spacing_m/1000.,GEODAT3D,
>      &			     GEODAT2D)
> 
> 	do J=nnyp,1,-nnyp/40
> 	write(6,733) (INT(GEODAT2D(I,J)),I=1,nnxp,nnxp/30)
> 	enddo
> 
>   733	format(50(I2,1x))
> 
> 
> 
> 	ENDIF
> 
>         ilndmsk=12
>         if(c10_grid_fname(1:lf).eq.'wrfsi')then
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
> c landmask for wrfsi
>            data(:,:,11)=GEODAT2D
>            data(:,:,ilndmsk) = 1.
>            where(data(:,:,11) .eq. 16.)data(:,:,ilndmsk)=0.
> 
> c grids 15 thru 39 are percent distributions
>            i=14
>            do j=1,24
>               data(:,:,i+j)=GEODAT3D(:,:,j)
>            enddo
> 
> 	ELSE ! rotlat case
> 
> 	ilndmsk=7
> 
> c landmask for wrfsi
>            data(:,:,6)=GEODAT2D
>            data(:,:,ilndmsk) = 1.
>            where(data(:,:,6) .eq. 16.)data(:,:,ilndmsk)=0.
> 
> ! putting land-use percentages in 19-42 for NMM
> 	I=18
>            do j=1,24
>               data(:,:,i+j)=GEODAT3D(:,:,j)
>            enddo
> 
> 	ENDIF
> 
> 
> 
>         else
> c landmask for laps
>            data(:,:,5)=GEODAT2D   !landuse for laps
>            data(:,:,ilndmsk)=1.
>            where(data(:,:,5) .eq. 16.)data(:,:,ilndmsk)=0.
>         endif
> 
> c ----------------------------------------------------------------
>         print*
>         print*,' Processing 10m land fraction data ....'
858,859d884
<  
<         GEODAT2D(:,:)=1.- GEODAT3D(:,:,16)
861,863c886,890
< c this function allows variable amounts of smoothing of
< c the land fraction data depending on grid spacing
< c
---
> 
> !mp	GEODAT3D(:,:,16) is water.  Should GEODAT3D range from 0-1?
> 
>         GEODAT2D(:,:)=1.- GEODAT3D(:,:,16)
> 
865c892,901
<         
---
> 
> 
> !!!	what does this filtering achieve?
> 
> !mp
> !mp	tossing into this test_hold array for temp storage
> 
> !	test_hold=GEODAT2D
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
868,869c904,905
<            call filter_2dx(GEODAT2D,nnxp,nnyp,1, 0.5)
<            call filter_2dx(GEODAT2D,nnxp,nnyp,1,-0.5)
---
>            call filter_2dx(geodat2d,nnxp,nnyp,1, 0.5)
>            call filter_2dx(geodat2d,nnxp,nnyp,1,-0.5)
872,873c908
<         idx=4
<         if(c10_grid_fname(1:lf).eq.'wrfsi')idx=10
---
> 	ENDIF
874a910,944
> !	do J=nnyp,1,-nnyp/40
> !	write(6,833) ((GEODAT2D(I,J)-test_hold(I,J)),
> !     &				I=1,nnxp,nnxp/12)
> !	enddo
> 
>   833	format(30(f6.4,1x))
> 
> c JS: 1-10-03:
> c These calls might become obsolete and the land_10m geog
> c data base will also become obsolete.
> c
> c       CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns)
> c    +,ytn(1,ns),deltax,deltay,GEODAT2D,adum,adum,adum
> c    +,PATH_TO_PCTL10M,1.,0.,new_DEM,1,istatus)
> c       print*,'Calling proc_geodat - land fraction'
> c       call proc_geodat(nnxp,nnyp,1,path_to_pctl10m
> c    +,lats(1,1,ns),lons(1,1,ns),data(1,1,ilndmsk)
> c    +,GEODAT2D,istatus_ldf)
> c       if(istatus_ldf .ne. 1)then
> c          write(6,*)' File(s) missing for 10m land frac data'
> c          write(6,*)' Static file not created.......ERROR'
> c          return
> c       endif
> 
>         if(c10_grid_fname(1:lf).eq.'wrfsi')then
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
>            idx=10
> 	ELSE
> 	   idx=5
> 	ENDIF
> 
>         else
>            idx=4
>         endif
880,885d949
< c Story here is that west coast terrain can be "funky" possibly due
< c to steepness and method to compute avg terrain.
< 
<         where(data(:,:,idx).le.0.1 .and. data(:,:,indxter).lt.5.0)
<      &data(:,:,indxter)=0.0
< 
887,888c951,952
< c -------------- Top Layer Soil Type -------------------------
< c type = O
---
>         where(data(:,:,idx).le.0.1 .and. topt_out(:,:).lt.5.0)
>      &topt_out=0.0
889a954,955
> c -------------------------------------------------------------------
> c
890a957
>         print*,'Elapsed time: ', itstatus
896a964,965
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
898c967
<      +,deltax,deltay,GEODAT2D,GEODAT3D,adum2d,adum2d
---
>      +,deltax,deltay,GEODAT2D,GEODAT3D,adum,adum
919a989,1003
> 	ELSE
> 
> 	call alt_categorical(nnxp,nnyp,hlat,hlon,
>      &			     PATH_TO_SOILTYPE_TOP_30S,
>      &			     type,16,grid_spacing_m/1000.,GEODAT3D,
>      &			     GEODAT2D)
> 
> 	write(6,*) 'soiltype_top'
> 	do J=nnyp,1,-nnyp/40
> 	write(6,733) (INT(GEODAT2D(I,J)),I=1,nnxp,nnxp/30)
> 	enddo
> 
> 	ENDIF
> 
> 
922a1007,1008
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
927a1014,1021
> 	ELSE ! 43-58
>            data(:,:,8)=GEODAT2D
> 	   I=42
>            do j=1,16
>               data(:,:,i+j)=GEODAT3D(:,:,j)
>            enddo
> 	ENDIF
> 
931,934c1025,1026
< c
< c -------------- Bottom Layer Soil Type -------------------------
< c type = O
< c
---
> 
> c ---------------------------------------------------------------
935a1028
>         print*,'Elapsed time: ', itstatus
937c1030
<         print*,' Processing 30s soil type bottom layer data'
---
>         print*,' Processing 30s soil type bottom layer data....'
938a1032,1034
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
940c1036
<      +,ytn(1,ns),deltax,deltay,GEODAT2D,GEODAT3D,adum2d,adum2d
---
>      +,ytn(1,ns),deltax,deltay,GEODAT2D,GEODAT3D,adum,adum
961a1058,1071
> 	ELSE
> 
> 	call alt_categorical(nnxp,nnyp,hlat,hlon,
>      &			     PATH_TO_SOILTYPE_BOT_30S,
>      &			     type,16,grid_spacing_m/1000.,GEODAT3D,
>      &			     GEODAT2D)
> 
> 	write(6,*) 'soiltype_bot'
> 	do J=nnyp,1,-nnyp/40
> 	write(6,733) (INT(GEODAT2D(I,J)),I=1,nnxp,nnxp/30)
> 	enddo
> 
> 	ENDIF
> 
967a1078,1079
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
972a1085,1097
> 
> 	ELSE !59-74
>             data(:,:,9)=GEODAT2D
> 
> 	write(6,*) 'set sbl values to : ', data(5,5,9),data(20,20,9)
>             i=58
>             do j=1,16
>                data(:,:,i+j)=GEODAT3D(:,:,j)
>             enddo
> 
> 	ENDIF
> 	    
> 	
976,979c1101,1102
< c
< c ------------- Greenness Fraction --------------------
< c type = G
< c
---
> 
> c ----------------------------------------------------------------
980a1104
>         print*,' Elapsed time: ', itstatus
982c1106
<         print*,' Calling PROC_GEODAT: Process green frac data.'
---
>         print*,' Calling GEODAT: Processing green frac data.'
985a1110
> 
991,993c1116
<               print*,'Error: Aborting process'
<               print*,'Error: Maybe not enough memory'
<               print*,'Error: Cannot allocate GEODAT3D'
---
>               print*,'Terminating: Maybe not enough memory'
996c1119
<               print*,'Warning: Continue without green fraction'         
---
>               print*,'Continue without greenness fraction'         
999a1123,1131
> c       CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns)
> c    +,ytn(1,ns),deltax,deltay,GEODAT2D,GEODAT3D
> c    +,adum,adum,path_to_green_frac,2.0,0.0,new_DEM,12
> c    +,istatus_grn)
> 
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 	
> 
1004a1137,1149
> 	ELSE
> 
> 	call alt_categorical_hemi(nnxp,nnyp,hlat,hlon,
>      &			     PATH_TO_GREEN_FRAC,
>      &			     type,12,grid_spacing_m/1000.,GEODAT3D,
>      &			     GEODAT2D)
> 
> 	istatus_grn=1
> 
> 
> 	ENDIF
> 	
> 
1009,1010c1154,1155
<             print*,' Error: File(s) missing for green frac data'
<             print*,' Error: Static file not created'
---
>             print*,' File(s) missing for green frac data'
>             print*,' Error:   Static file not created'
1015,1016c1160,1162
<             print*,' Warning: File(s) missing for green frac data'
<             print*,' Warning: green frac not added to static file'
---
>             print*,' File(s) missing for green frac data'
>             print*,'           *** WARNING ***'
>             print*,' green fraction data not added to static file'
1023c1169,1175
<            igrn=83
---
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
>            i=83
> 	ELSE
> 	   i=74
> 	ENDIF
> 	
1025c1177
<            igrn=12
---
>            i=12
1028c1180
<            data(:,:,igrn+j)=GEODAT3D(:,:,j)
---
>            data(:,:,i+j)=GEODAT3D(:,:,j)
1030,1033c1182,1186
< c
< c --------------- Deep Soil Temperature -------------------
< c type = T
< c
---
> 
> 44      continue
> c ------------------------------------------------------------------
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
1034a1188
>         print*,' Elapsed time: ', itstatus
1036c1190
<         print*,' Call PROC_GEODAT: Process 1 deg soiltemp data.'
---
> c       print*,' Calling GEODAT: Processing 1 degree soiltemp data.'
1038c1192,1194
<         GEODAT2D=0.0
---
> c       CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns)
> c    +,ytn(1,ns),deltax,deltay,GEODAT2D,adum,data(1,1,ilndmsk),adum  !send in landmask to help 
> c    +,path_to_soiltemp_1deg,2.0,0.0,new_DEM,1,istatus_tmp)
1043a1200,1219
> 	ELSE
> 
> 	CALL alt_10by10(nnxp,nnyp,10,grid_spacing_m/1000.,
>      &	hlat,hlon,type,PATH_TO_SOILTEMP_1DEG,GEODAT2D,adum,
>      &	adum,adum)
> 	
> 
> !	convert to degrees
> 
> 	DO J=1,NNYP
> 	DO I=1,NNXP
> 		GEODAT2D(I,J)=GEODAT2D(I,J)/100.
> 	ENDDO
> 	ENDDO
> 
> 	istatus_tmp=1
> 
> 	ENDIF
> 	
> 
1048,1049c1224,1225
<             print*,' Error:  File(s) missing for soiltemp data' 
<             print*,' Error:  Static file not created'
---
>             print*,' File(s) missing for soiltemp data' 
>             print*,' Error:   Static file not created'
1062d1237
<         idst=25
1064c1239,1245
<            idst=96
---
> 	    IF (c6_maproj .ne. 'rotlat') THEN
>            	i=96
> 	    ELSE
> 	   	i=87
> 	    ENDIF
>         else
>            i=25
1066,1073d1246
<         data(:,:,idst)=GEODAT2D
< c
< c -------------- Terrain slope index categories -------------------------
< c type = I
< c
<         itstatus=ishow_timer()
<         print*
<         print*,' Processing 1 deg islope data'
1075,1110c1248
<         allocate (adum3d(nnxp,nnyp,9))  !temporarily holds %dist of islope cat
< 
<         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)
< 
<         deallocate (adum2d,adum3d)
< 
<         if(istatus_slp.ne.1)then
<            print*
<            print*,' ter slope category data not processed completely'
<            if(c10_grid_fname(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_fname(1:lf).eq.'wrfsi')then
<            islp=109
<         else
<            islp=39
<         endif
<         data(:,:,islp)=GEODAT2D
---
>         data(:,:,i)=GEODAT2D
1115a1254,1257
>         if(.true.)then
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
1118,1121c1260,1262
<         call adjust_geog_data(nnxp,nnyp,12,istatus_grn,
<      &istatus_tmp,istatus_slp,lats(1,1,ns),data(1,1,indxter),
<      &path_to_soiltemp_1deg, data(1,1,ilndmsk),data(1,1,idst),
<      &GEODAT3D,data(1,1,islp),istatus)
---
>         call adjust_geog_data(nnxp,nnyp,12,istatus_grn
>      &,istatus_tmp,lats(1,1,ns),topt_out,path_to_soiltemp_1deg
>      &,data(1,1,ilndmsk),data(1,1,i),GEODAT3D,istatus)
1123c1264
<            print*,' Warning: Processing incomplete: adjust_geog_data'
---
>            print*,'Processing incomplete in adjust_geog_data'
1125c1266
<               print*,'Error: wrfsi static file not generated'
---
>               print*,'Error: wrfsi static not generated'
1129,1151c1270
< 
< 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
< 
< 
< c rewrite adjusted greenfrac into output array
<         do j=1,12
<            data(:,:,igrn+j)=GEODAT3D(:,:,j)
<         enddo
< c
< c annual max/min greenness fraction in domain
< c
<         if(c10_grid_fname(1:lf).eq.'wrfsi')then
<            print*,' compute max/min greenness frac at grid points'
<            do j=1,nnyp
<            do i=1,nnxp
<               data(i,j,110)=MAXVAL(GEODAT3D(i,j,:))
<               data(i,j,111)=MINVAL(GEODAT3D(i,j,:))
<            enddo
<            enddo
---
> 	ENDIF
1153,1156c1272
< c
< c ---------------Monthly Albedo Climatology--------------------------------
< c type = A
< c
---
> c ------------------------------------------------------------------
1157a1274
>         print*,' Elapsed time: ', itstatus
1159c1276
<         print*,' Processing albedo climo data'
---
> c       print*,' Calling GEODAT: Processing albedo climo data.'
1160a1278,1283
> c       CALL GEODAT(nnxp,nnyp,erad,90.,std_lon,xtn(1,ns)
> c    +,ytn(1,ns),deltax,deltay,GEODAT2D,adum,data(1,1,ilndmsk),adum  !send in landmask to help
> c    +,path_to_soiltemp_1deg,2.0,0.0,new_DEM,1,istatus_tmp)
> 
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
1164a1288,1297
> 	ELSE
> 
> 	call alt_categorical_hemi(nnxp,nnyp,hlat,hlon,
>      &			     PATH_TO_ALBEDO,
>      &			     type,12,grid_spacing_m/1000.,GEODAT3D,
>      &			     GEODAT2D)
> 
> 	istatus_alb=1
> 	ENDIF
> 
1169,1170c1302,1303
<          print*,' Error: Albedo climo data not processed completely'
<          print*,' Error: File(s) missing for albedo data'
---
>          print*,'Albedo climo data not processed completely'
>          print*,' File(s) missing for albedo data'
1175a1309,1317
> c
> c We can force water points to have albedo = 0.08.
> c Land points with missing albedo are given the mean domain
> c albedo.
> c We can be more sophisticated with this by
> c determining the albedo as a function of landuse category and
> c then assigning missing albedo with the same landuse albedo.
> c
> c       print*,'force albedo = 0.08 for water points'
1177,1178d1318
< c force water points to 0.08 albedo
< 
1183a1324
> 	IF (c6_maproj .ne. 'rotlat') THEN
1184a1326,1328
> 	ELSE
> 	   i=87
> 	ENDIF
1191,1194c1335,1343
< c
< c ---------------- Max Snow Albedo ------------------
< c type = M
< c
---
> 
> c       call get_static_albedo(nnxp,nnyp,lats(1,1,ns),lons(1,1,ns)
> c    +,data(1,1,12),GEODAT2D,istatus)
> c       if(c10_grid_fname(1:lf).eq.'wrfsi')then
> c          data(:,:,47)=GEODAT2D
> c       else
> c          data(:,:,6)=GEODAT2D
> c       endif
> 
1203,1217c1352,1353
<          print*,'Error: Max Snow Albedo not processed completely'
<          if(c10_grid_fname(1:lf).eq.'wrfsi')then
<             print*,' Error: Static file not created '
<             print*
<             istatus=0
<             return
<          else
<             print*,'Warning: File(s) missing for max snow albedo'
<          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
---
>          print*,'Max Snow Albedo data not processed completely'
>          print*,'Warning: File(s) missing for max snow albedo data'
1219,1222c1355,1358
<         if(c10_grid_fname(1:lf).eq.'wrfsi')then
<            data(:,:,47)=GEODAT2D
<         else
<            data(:,:,6)=GEODAT2D
---
> c        print*,' Error: Static file not created '
> c        print*
> c        istatus=0
> c        return
1224a1361,1364
> c Adjust geog data to conform to landuse data (the most
> c accurate for defining land-water boundaries).
> c Adjust soil temps to terrain elevations.
> 
1231a1372,1374
> !mp
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
1240a1384
>            data(:,:,9)=topt_out
1243a1388,1389
> 
> 
1279c1425,1434
<         else
---
> 	ELSE ! rotlat case
> 	
> 	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
1280a1436,1507
>            data(:,:,18)=topt_out
> 
> c compute map/grid information. Reallocate GEODAT3D array
>            allocate (GEODAT3D(nnxp,nnyp,4))
> 
> 
> !!!!	figure this business out
> 
>            call vecrot_rotlat(nnxp,nnyp,mdlat,mdlon
>      +,vlat,vlon,GEODAT3D(:,:,1),GEODAT3D(:,:,2))
> 
> 	write(6,*) 'sample rot values: '
> 	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
> 
> 	
>            i=10
> 
> !	sin, then cosine?
>            data(:,:,i)  =GEODAT3D(:,:,2)
>            data(:,:,i+1)=GEODAT3D(:,:,1)
>  
> !          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           
> 
> !!!!	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
> 
> !mp	believe need alternate version to handle different 
> !	variable set
> 
> !030513  Actually just made modifications within the get_gridgen_var code
> 
> 	write(6,*) 'calling get_gridgen_var with ngrids= ', ngrids
>            call get_gridgen_var(ngrids,ngrids,var,comment)
> 
> !	do N=1,ngrids
> !	write(6,*) 'N, var(N): ', n,var(N)
> !	enddo
> 
> 	ENDIF
> 
> 
>         else ! non-wrfsi case
> 
1305a1533,1535
> !!mp	what does check_domain do?
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 
1315a1546,1547
> 	
> 	ENDIF
1316a1549,1550
> 	write(6,*) 'call put_laps_static'
> 
1320c1554
<      1,c6_maproj,deltax,deltay,nest)
---
>      1,c6_maproj,deltax,deltay)
1331,1333c1565,1567
< c
< c ================================================================
< c
---
> 
> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
> 
1530,1531c1764
<       if(ofn(len-1:len-1).eq.'V'.or.
<      &   ofn(len-1:len-1).eq.'I')then
---
>       if(ofn(len-1:len-1).eq.'V')then
1534c1767,1769
<          if(ofn(len-1:len-1).eq.'I')cdatatype='islope'
---
>       elseif(ofn(len-1:len-1).eq.'A')then
>          icnt=0
>          cdatatype='albedo'
1539c1774,1775
<      &       ofn(len-1:len-1).eq.'H')then
---
>      &       ofn(len-1:len-1).eq.'H' .or.
>      &       ofn(len-1:len-1).eq.'L')then
1540a1777,1784
>          if(ofn(len-1:len-1).eq.'L')cdatatype='landfrac'
>       elseif(ofn(len-1:len-1).eq.'G')then
>          icnt = 0
>          cdatatype='greenfrac'
>          lcat = 12
>       elseif(ofn(len-1:len-1).eq.'T')then
>          icnt = 0
>          cdatatype='soiltemp'
1698,1702c1942,1943
<                     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
---
>                     if( (ofn(len-1:len).eq.'U').and.(no.eq.1200) )then
>                       CALL READ_DEM(29,TITLE3(1:LB),no,no,2,2, ! world topo_30s
1704,1707d1944
<                          else
<                             CALL READ_DEM(29,TITLE3(1:LB),no,no,2,2, ! world topo_30s
<      .                              DATO(1,1,NOFR,1))
<                          endif
1718,1719c1955
<      .                      (ofn(len-1:len).eq.'I').or.
<      .                      (ofn(len-1:len).eq.'A') )then      ! greenfrac/albedo/islope
---
>      .                      (ofn(len-1:len).eq.'A') )then      ! greenfrac/albedo
1834a2071
>      &         cdatatype(1:lent).eq.'landfrac'.or.
1883,1884c2120
<             elseif(cdatatype(1:lent).eq.'islope'    .or.
<      &             cdatatype(1:lent).eq.'landuse'   .or.
---
>             elseif(cdatatype(1:lent).eq.'landuse'   .or.
1901a2138,2142
> c           elseif(cdatatype(1:lent).eq.'soiltemp'    )then
> c             call compute_categories(cdatatype,np*np,DATP(1,1,1)
> c    &                               ,1,domcat,pctcat)
> c             datq(iq,jq)=domcat
> 
1907a2149
> c     stop
1912c2154,2155
<       if(cdatatype(1:lent).eq.'topography')then
---
>       if(cdatatype(1:lent).eq.'topography'.or.
>      +   cdatatype(1:lent).eq.'landfrac')then
1941a2185,2194
> c let's constrain the landfrac data between 0.0 and 1.0
>            if(ofn(len-1:len).eq.'L')then
>               if(DATR(IR,JR).gt. 1.0)then
>                  DATR(IR,JR)=1.0
>               endif
>               if(DATS(IR,JR,1).gt.1.0)then
>                  DATS(IR,JR,1)=1.0
>               endif
>            endif
> 
1953d2205
<      +       cdatatype(1:lent).eq.'islope' .or.
1970a2223,2267
>       elseif(cdatatype(1:lent).eq.'greenfrac' )then
> 
>         DO 47 LP=1,lcat
>          DO 48 JR=1,N3
>           DO 49 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
>             dats(ir,jr,lp)=datqs(ixr,iyr,lp)  !monthly categories for greenness frac
> 
>  49       CONTINUE
>  48      CONTINUE
>  47     CONTINUE
> 
>       elseif(cdatatype(1:lent).eq.'soiltemp'  )then
> 
>         where(datq .eq. 0.0)datq=r_missing_data
>         datr=r_missing_data
>         DO 58 JR=1,N3
>          DO 59 IR=1,N2
> c           IXR=NINT((XT(IR)-XQ1)/DELTAXQ)+1.
> c           IYR=NINT((YT(JR)-YQ1)/DELTAYQ)+1.
> c           if(ixr.lt.1)ixr=1
> c           if(iyr.lt.1)iyr=1
> c           if(ixr.gt.n2)ixr=niq
> c           if(iyr.gt.n3)iyr=njq
> c           datr(ir,jr)=datq(ixr,iyr)
> 
>             XR=(XT(IR)-XQ1)/DELTAXQ+1.
>             YR=(YT(JR)-YQ1)/DELTAYQ+1.
> c           call bilinear_interp_extrap(xr,yr,niq,njq
> c    1               ,datq,rval,istatus)
> 
> c           CALL GDTOST2(DATQ,NIQ,NJQ,XR,YR,RVAL)
> 
>             DATR(IR,JR)=max(0.,RVAL)
> 
>  59      CONTINUE
>  58     CONTINUE
> 
> 
2027,2032c2324,2353
<       subroutine vfirec(iunit,a,n,type)
<       character*1 vc
<       character*(*) type
<       common/vform/vc(0:63)
<       character line*80, cs*1
<       dimension a(*)
---
>       subroutinevfirec(iunit,a,n,type)                                  
>       character*1vc                                                     
>       character*(*)type                                                 
>       common/vform/vc(0:63)                                             
>       characterline*80,cs*1                                             
>       dimensiona(*)                                                     
>                                                                         
>       if(vc(0).ne.'0')callvfinit                                        
>                                                                         
>       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                                                      
>       do20i=1,n,nvalline                                                
>       read(iunit,'(a78)', end=15)line                                          
>       go to 16
2034,2041c2355
<       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')
---
>  15   write(6,*)' Warning, incomplete terrain file detected'
2043,2050c2357,2391
<       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
---
>  16   ic=0                                                              
>       do30ii=i,i+nvalline-1                                             
>       isval=0                                                           
>       if(ii.gt.n)goto20                                                 
>       do40iii=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=intor(intlshft(nc,6*(nchs-iii)),isval)                      
>  40   continue                                                          
>       a(ii)=isval                                                       
>  30   continue                                                          
>  20   continue                                                          
>                                                                         
>       facti=1./fact                                                     
>       if(type.eq.'LIN')then                                             
>       do48i=1,n                                                         
>       a(i)=a(i)*facti-bias                                              
>  48   continue                                                          
>       elseif(type.eq.'LOG')then                                         
>       scfct=2.**(nbits-1)                                               
>       do55i=1,n                                                         
>       a(i)=sign(1.,a(i)-scfct)                                          
>      +*(10.**(abs(20.*(a(i)/scfct-1.))-10.))                            
>  55   continue                                                          
>       endif                                                             
>       return                                                  
>       end                            
2052,2094d2392
<       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
2191,2193c2489
<       if(ctype.eq.'landuse'   .or.
<      +   ctype.eq.'soiltype'  .or.
<      +   ctype.eq.'islope'  )then
---
>       if(ctype.eq.'landuse'.or.ctype.eq.'soiltype')then
2255,2256c2551,2552
<      &,istatgrn,istattmp,istatslp,lat,topt_out,path_to_soiltemp
<      &,landmask,soiltemp_1deg,greenfrac,islope,istatus)
---
>      &,istatgrn,istattmp,lat,topt_out,path_to_soiltemp
>      &,landmask,soiltemp_1deg,greenfrac,istatus)
2272d2567
<       integer istatslp
2274d2568
<       integer isum
2291d2584
<       real,intent(inout)  ::    islope(nnxp,nnyp)
2295d2587
<       real,   allocatable ::    islopetmp     (:,:)
2300d2591
<       real    avgslpcat
2305d2595
<       real    islp
2312,2313c2602,2603
<       if(istatgrn.eq.0.and.istattmp.eq.0.and.istatslp.eq.0)then
< 	 print*,'Unable to process greenfrac soiltemp or slp ',
---
>       if(istatgrn.eq.0 .and. istattmp.eq.0)then
> 	 print*,'Unable to process greenfrac or soiltemp ',
2322,2323c2612
<      &          soiltmp(nnxp,nnyp),
<      &          islopetmp(nnxp,nnyp))
---
>      &          soiltmp(nnxp,nnyp))
2329d2617
<       islopetmp=islope
2335d2622
<       where(islopetmp.eq.0.0)islopetmp=r_missing_data
2359,2363d2645
<          if(islopetmp(i,j).ne. r_missing_data)then
<             isum=islopetmp(i,j)+isum
<             islp=islp+1
<          endif
< 
2392,2408c2674,2681
<          if(tatlmx.ne.0.and.tatlmn.ne.0.and.
<      .      tatlmx.ne.tatlmn)then
<             tslp=(rlatmx-rlatmn)/(tatlmx-tatlmn)
<             do j=1,nnyp
<             do i=1,nnxp
<                soiltemp_1deg(i,j)=tatlmn+(rlatmx-lat(i,j))/tslp
<                sumt=soiltemp_1deg(i,j)+sumt
<             enddo
<             enddo
<             avgtmp=sumt/(nnyp*nnxp)
<          elseif(tatlmx.gt.0)then
<             avgtmp=tatlmx
<          elseif(tatlmn.gt.0)then
<             avgtmp=tatlmn
<          else
<             avgtmp=0.0
<          endif
---
>          tslp=(rlatmx-rlatmn)/(tatlmx-tatlmn)
>          do j=1,nnyp
>          do i=1,nnxp
>             soiltemp_1deg(i,j)=tatlmn+(rlatmx-lat(i,j))/tslp
>             sumt=soiltemp_1deg(i,j)+sumt
>          enddo
>          enddo
>          avgtmp=sumt/(nnyp*nnxp)
2419,2420d2691
<          else
<             avggrn(l)=0.0
2425,2433d2695
<       avgslpcat=r_missing_data
<       if(islp.gt.0)then
<          avgslpcat=float(isum)/islp
<          print*,'Average slope category: ',avgslpcat
<       else
<          print*,'Could not compute average slpe cat ?'
<          print*,'Average slope category: ',avgslpcat
<       endif
< 
2438c2700
<       ijsthresh = int(nnxp/2)
---
>       ijsthresh = int(nnxp/4)
2541,2544d2802
< c ------------------------------------------
< c this section for terrain slope categories
< c
<          if(islopetmp(i,j).eq.r_missing_data)then !an inconsistency
2546,2596d2803
<           endsearch = .false.
< 
<           sumt=0.0
<           isc=0
<           islp=-9
<           is=1
<           js=1
<           jj=j
< 
<           do while (.not.endsearch)
<            do ii=i-is,i+is
<            do jj=j-js,j+js
< 
<             if( (ii.ge.1) .and. (ii.le.nnxp)
<      &     .and.(jj.ge.1) .and. (jj.le.nnyp)) then
< 
<              if(landmask(ii,jj).eq.1.and.
<      &         islopetmp(ii,jj).lt.r_missing_data)then
<                islp=islopetmp(ii,jj)
<                sumt=sumt+islopetmp(ii,jj)
<                isc=isc+1
<              endif
< 
<             endif
< 
<            enddo
<            enddo
< 
<            if(isc.gt.0)then
<               if(islp.gt.7.0)islp=13.0    !force mean to glacial ice as necessary
<               if(sumt.lt.isc)then
<                  islp=1.0
<               else
<                  islp=float(nint(sumt/float(isc)))
<               endif
<               islope(i,j)=islp
<               endsearch=.true.
<            else
<               is=is+1
<               js=js+1
<               if(is.gt.ijsthresh.or.js.gt.ijsthresh)endsearch=.true.
<            endif
< 
<           enddo
< 
<          else
< 
<           islope(i,j)=islopetmp(i,j)
< 
<          endif
< 
2610,2613d2816
<          if(islope(i,j).ne.r_missing_data)then
<             islope(i,j)=r_missing_data
<          endif
< 
2619c2822
<       deallocate (grnfrctmp,soiltmp,islopetmp)
---
>       deallocate (grnfrctmp,soiltmp)
2627c2830
<          if(landmask(i,j).eq.1)then  !a land point
---
>          if(landmask(i,j).eq.1)then
2639,2642d2841
<             if(islope(i,j).eq.r_missing_data)then
<                islope(i,j)=float(nint(avgslpcat))
<             endif
< 
2648,2649d2846
<       where(islope.eq.r_missing_data)islope=0.0  !put water category back to original value
< 
2653a2851,2920
> 
> CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
>        subroutine vecrot_rotlat(IM,JM,TPH0D,TLM0D,VLAT,VLON,DRAY,CRAY)
> 
> C  SUBPROGRAM: ROTLLE        ROTATE WINDS ON LAT/LONG GRID TO E-GRID
> C   PRGMMR: T.BLACK         ORG: W/NP22     DATE: ??-??-??
> C
> C ABSTRACT: ROTATES WINDS ON THE LAT/LONG GRID TO THE ETA MODEL GRID
> C               (same angles can be used to do reverse function)
> C
> C PROGRAM HISTORY LOG:
> C   ??-??-??  T.BLACK
> C   98-06-08  M.BALDWIN - CONVERT TO 2-D CODE
> C   01-01-03  T BLACK - MODIFIED FOR MPI
> C
> C   INPUT ARGUMENT LIST:
> C     VLAT     - LATITUDE OF E-GRID V POINTS (DEGREES)
> C     VLON     - LONGITUDE OF E-GRID V POINTS (DEGREES)
> C
> C   OUTPUT ARGUMENT LIST:
> C     DRAY     - ROTATION COSINE
> C     CRAY     - ROTATION SINE
> C
> C ATTRIBUTES:
> C   LANGUAGE: FORTRAN 90
> C   MACHINE:  IBM RS/6000 SP
> C
> C***
> C*** ROTATE THE LAT-LON WINDS TO/FROM THE E-GRID
> C***
> C
> C    N O T E : INPUT LAT/LONG MUST BE CONVERTED TO RADIANS !!!
> C
> C----------------------------------------------------------------------
>         parameter(D2RAD=1.745329E-2)
> 
>                          D I M E N S I O N
>      1 VLAT(IM,JM),VLON(IM,JM)
>      2,CRAY(IM,JM),DRAY(IM,JM)
> C----------------------------------------------------------------------
> 
>         ERPHI0=TPH0D*D2RAD
>         if (TLM0D .lt. 0) then
>         ERLAM0=(TLM0D+360.)*D2RAD
>         else
>         ERLAM0=TLM0D*D2RAD
>         endif
> 
>       SPHI0 = SIN(ERPHI0)
>       CPHI0 = COS(ERPHI0)
> C
> 
>       DO J = 1, JM
>       DO I = 1, IM
>         TLAT = VLAT(I,J) * D2RAD
>         TLON = VLON(I,J) * D2RAD
>         RELM = TLON - ERLAM0
>         SRLM = SIN(RELM)
>         CRLM = COS(RELM)
>         SPH = SIN(TLAT)
>         CPH = COS(TLAT)
>         CC = CPH * CRLM
>         TPH = ASIN(CPHI0 * SPH - SPHI0 * CC)
>         RCTPH = 1.0 / COS(TPH)
>         CRAY(I,J) = SPHI0 * SRLM * RCTPH
>         DRAY(I,J) = (CPHI0 * CPH + SPHI0 * SPH * CRLM) * RCTPH
>       ENDDO
>       ENDDO
> C----------------------------------------------------------------------
> 	end subroutine vecrot_rotlat
done gridgen_model.f ======================
COMPARE gridgen_utils.f
39a40
> 
58c59
<      +  NNXP,NNYP,1.)
---
>      +  NNXP,NNYP)
72a74
>             call get_grid_center(mdlat,mdlon,istatus)
81a84
>             call get_grid_center(mdlat,mdlon,istatus)
340c343
<       if(ngrids.eq.38)then  !this is the LAPS analysis section
---
>       if(ngrids.eq.38)then
347c350
<          var(6)    = 'ALB'  !now used for max snow alb 2-20-03 JS.
---
>          var(6)    = 'ALB'  !possibly not used now.
373,374c376,377
<          comment(1) = 'Lat: From MODEL by J. Smart/ S. Albers 2-03\0'
<          comment(2) = 'Lon: From MODEL by J. Smart/ S. Albers 2-03\0'
---
>          comment(1) = 'Lat: From MODEL by J. Snook/ S. Albers 1-95\0'
>          comment(2) = 'Lon: From MODEL by J. Snook/ S. Albers 1-95\0'
376c379
<          comment(4) = 'Land Fraction: derived from USGS land use \0'
---
>          comment(4) = 'Land Fraction \0'
378c381
<          comment(6) = 'Maximum Snow Albedo; defined over land only \0'
---
>          comment(6) = 'Clear Sky Albedo - fixed at .08 over water\0'
406c409
<       elseif(ngrids.eq.112)then   !this is the WRFSI section
---
>       elseif(ngrids.eq.109)then
414a418
> 	write(6,*) 'setting things here.....var(7) is LAC'
443c447
<          var(i+8)   = 'ALB'  ! Maximum snow albedo defined over land only
---
>          var(i+8)   = 'ALB'  ! Static (climatological) albedo
478c482
<             var(i+j)= 'A'//cat ! monthly albedo
---
>             var(i+j)= 'A'//cat   ! monthly albedo
481,484d484
<          var(109) = 'SLP'      !terrain slope index, dominant category
<          var(110) = 'GNX'      !max greenness fraction
<          var(111) = 'GNN'      !min greenness fraction
< 
487,488c487,488
<          comment(1) = 'Made from MODEL by J. Smart/ S. Albers 2-03\0'
<          comment(2) = 'Made from MODEL by J. Smart/ S. Albers 2-03\0'
---
>          comment(1) = 'Made from MODEL by J. Snook/ S. Albers 1-95\0'
>          comment(2) = 'Made from MODEL by J. Snook/ S. Albers 1-95\0'
496c496
<          comment(10)= 'Land Fraction computed from USGS landuse\0'
---
>          comment(10)= 'Land Fraction A-grid \0'
521c521
<          comment(i+8)= 'Maximum Snow Albedo (%) over land only \0'
---
>          comment(i+8)= 'Static Albedo (%) valid only over water atm \0'
560,562c560
<          comment(109)=  'terrain slope index'
<          comment(110)=  'annual max greenness fraction'
<          comment(111)=  'annual min greenness fraction'
---
>          comment(ngrids)= '\0'
563a562,691
> 	ELSEIF (ngrids .eq. 100) THEN ! NMM (rotlat) case
> 
>          var(1)   = 'LAH'  ! Mass(H)-point lats
>          var(2)   = 'LOH'  ! Mass(H)-point lons
>          var(3)   = 'LAV'  ! Mass(V)-point lats
>          var(4)   = 'LOV'  ! Mass(V)-point lons
>          var(5)   = 'LDF'  ! Land Fraction
>          var(6)   = 'USE'  ! landuse dominant category
>          var(7)   = 'LND'  ! Land-Water Mask based on USGS landuse
>          var(8)   = 'STL'  ! Soiltype top layer dominant category
>          var(9)   = 'SBL'  ! Soiltype bot layer dominant category
> 
>          i=10
> 
>          var(i)     = 'SPR'  ! Sin(projection rotation) from true
>          var(i+1)   = 'CPR'  ! Cos(projection rotation) from true
>          var(i+2)   = 'CPH'  ! Horizontal component of coriolis parameter
>          var(i+3)   = 'CPV'  ! Vertical component of coriolis parameter
>          var(i+4)   = 'ALB'  ! Static (climatological) albedo
>          var(i+5)   = 'STD'  ! Standard Deviation of Elevation Data (m)
>          var(i+6)  = 'SLN'  ! Terrain Slope; Longitudinal Component (m/m)
>          var(i+7)  = 'SLT'  ! Terrain Slope; Latitudinal Component (m/m)
>          var(i+8)  = 'AVC'  ! Topo (m) on c-stagger grid
> 
> 	i=18
>          do j=1,24
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           var(i+j)= 'U'//cat
>          enddo
> 
>          i=42
>          do j=1,16
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           var(i+j)= 'T'//cat   !top layer (0-30cm) soiltype (% dist)
>          enddo
>          i=58
>          do j=1,16
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           var(i+j)= 'B'//cat   ! bot layer (30-90cm) soiltype (% dist)
>          enddo
>          i=74
>          do j=1,12
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           var(i+j)= 'G'//cat   ! vegetation greenness fraction
>          enddo
> 
>          var(87)='TMP'
> 
>          i=87
>          do j=1,12
>             write(cat,'(i2.2)')j
>             if(cat(1:1).eq.' ')cat(1:1)='0'
>             if(cat(2:2).eq.' ')cat(2:2)='0'
>             var(i+j)= 'A'//cat   ! monthly albedo
>          enddo
> 
>          var(ngrids)   = 'ZIN'
> 
>          comment(1) = 'H-point (mass) latitudes \0'
>          comment(2) = 'H-point (mass) longitudes \0'
>          comment(3) = 'V-point (mass) latitudes \0'
>          comment(4) = 'V-point (mass) longitudes \0'
>          comment(5)= 'Land Fraction  \0'
>          comment(6)= 'Land Use Dominant category \0'
>          comment(7)= 'Land-Water Mask (0=water; 1 otherwise) \0'
>          comment(8)= 'Soiltype Top Layer dominant category \0'
>          comment(9)= 'Soiltype Bot Layer dominant category \0'
> 
>          i=10
>          comment(i)= 'Sin of projection rotation (rad) \0'
>          comment(i+1)= 'Cosine of projection rotation (rad) \0'
>          comment(i+2)= 'Horizontal component coriolis parameter \0'
>          comment(i+3)= 'Vertical component coriolis parameter \0'
>          comment(i+4)= 'Static Albedo (%) valid only over water atm \0'
>          comment(i+5)= 'Standard Deviation of Elevation data (m)\0'
>          comment(i+6)= 'Mean longitudinal terrain slope (m/m)\0'
>          comment(i+7)= 'Mean latitudinal terrain slope (m/m)\0'
>          comment(i+8)= 'Average terrain elevation (c-stagger) (m)\0'
> 
>          i=18
>          do j=1,24
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           var(i+j)= 'U'//cat
>           comment(i+j)= '% Dist Land Use Category '//cat//' \0'
>          enddo
> 
> 
>          i=42
>          do j=1,16
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           comment(i+j)='% Dist Top Layer Soiltype Category '//cat
>          enddo
> 
>          i=58
>          do j=1,16
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           comment(i+j)= '% Dist Bot Layer Soiltype Category '//cat
>          enddo
> 
>          i=74
>          do j=1,12
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           comment(i+j)= 'vegetation greenness fraction: mon = '//cat
>          enddo
> 
>          comment(87)='1 degree mean annual soiltemp (deg K)'
>          i=87
>          do j=1,12
>           write(cat,'(i2.2)')j
>           if(cat(1:1).eq.' ')cat(1:1)='0'
>           if(cat(2:2).eq.' ')cat(2:2)='0'
>           comment(i+j)= 'climatological albedo: mon = '//cat
>          enddo
> 
565a694
> 
598,612c727,733
<            if(nn1.eq.1201)then
<               open(unit_no,file=unit_name,status='old',
<      .form='unformatted')
<               read(unit_no)data
<               close(unit_no)
<            else
<               call read_binary_field(cdata,i1,i2,nn1*nn2,unit_name,len)
< 
<               do county=1,nn2
<               do countx=1,nn1
<                idata(countx,county) = ia (cdata(countx,county),2,0)
<               enddo
<               enddo
<               if(ctiletype.eq.'T')multiplier=.01  !(for T data these are temps * 100)
<            endif
---
>            call read_binary_field(cdata,i1,i2,nn1*nn2,unit_name,len)
>            do county=1,nn2
>            do countx=1,nn1
>               idata(countx,county) = ia (cdata(countx,county),2,0)
>            enddo
>            enddo
>            if(ctiletype.eq.'T')multiplier=.01  !(for T data these are temps * 100)
617d737
<         if(nn1.le.1200)then
631d750
<         endif
644c763
<         integer  len, lend, i1, i2, i
---
>         integer  len, i1, i2
653d771
<         character*1   ctype
667,673d784
< 
<         call s_len(unit_name,len)
<         call get_directory_length(unit_name,lend)
<         ctype=unit_name(lend+1:lend+1)
< 
<         print*,'read_dem_g: tile type = ',ctype
< 
676,680c787
< c       if(nn1.ne.1250 .and. nn2.ne.1250)then
<         if(ctype.ne.'A' .and. 
<      &     ctype.ne.'G' .and.
<      &     ctype.ne.'I' .and.
<      &     ctype.ne.'M')then
---
>         if(nn1.ne.1250 .and. nn2.ne.1250)then
715,722d821
< c we'll resurrect the actually categories when done
< c processing but for now we want categories 1-9 for
< c terrain slope index.
<         if(ctype == 'I')then
<            where(data .eq. 13)data = 8
<            where(data .eq. 0)data = 9
<         endif
< 
805c904
<        SUBROUTINE POLAR_GP(LAT,LON,X,Y,DX,DY,NX,NY,DIR)
---
>        SUBROUTINE POLAR_GP(LAT,LON,X,Y,DX,DY,NX,NY)
812,814c911
<        INTEGER   NX,NY
<        INTEGER   IDIR  !positive going from center Lat/Lon to SW X/Y;
< C                       negative going from SW Lat/Lon to center X/Y.
---
>        INTEGER*4 NX,NY
823a921,932
>        TLAT=90.0
>        call get_standard_longitude(std_lon,istatus)
>        if(istatus .ne. 1)then
>            write(6,*)' Error calling laps routine'
>            stop 
>        endif
>        TLON=std_lon
> C
> C      CALL GETOPS(PLAT,PLON,LAT,LON,TLAT,TLON)
> C      CALL PSTOXY(XDIF,YDIF,PLAT,PLON,ERAD)
> 
> C      call latlon_to_xy(LAT,LON,TLAT,TLON,ERAD,XDIF,YDIF)
826,828c935,938
<        X=XDIF+(1.-DIR*FLOAT(NX)/2.)*DX
<        Y=YDIF+(1.-DIR*FLOAT(NY)/2.)*DY
<  
---
> C
>        X=XDIF+(1.-FLOAT(NX)/2.)*DX
>        Y=YDIF+(1.-FLOAT(NY)/2.)*DY
> C
830c940
<  
---
> C
1135,1139d1244
< c
< c ---------------------------------------------------------
< c
<       subroutine eval_localization(cstaticdir,nest,localize
<      .,cgrid_fname,La1_dom,Lo1_dom,istatus)
1141,1253d1245
< c routine performs the following:
< c 1. tests if the static file has been generated for this domain or nest?
< c 2. tests if the static variables consistent with the current namelist specs
< c
< c if 1 or 2 is false then we need to either localize or re-localize
< c the domain and "localize" is returned indicating such (true or false).
< c
<       implicit  none
< 
<       integer   nest
<       integer   Nx,Ny
<       integer   nx_dom,ny_dom
<       integer   lf
<       integer   istatus
< 
<       character cstaticdir*200
<       character cgrid_fname*10
<       character cnest*2
< 
<       logical   localize
<       logical   static_exists
< 
<       real      Dx,Dy,La1,Lo1,LoV,Latin1,Latin2
<       real      grid_spacing_dom_m,grid_spacing_m
<       real      La1_dom,Lo1_dom
< 
<       localize = .false.
< 
<       call rd_static_attr(cstaticdir,nest,cgrid_fname
<      .,Nx, Ny, Dx, Dy, La1, Lo1, Latin1, Latin2, LoV
<      .,istatus)
< 
<       IF (istatus .NE. 1) THEN
<         print '(A,I5)', ' Error reading WRF static file: ',istatus
<         localize = .true.
<         return
<       END IF 
< 
<       call get_grid_spacing(grid_spacing_dom_m,istatus)
<       call get_grid_dim_xy(nx_dom,ny_dom,istatus)
<       if(Nx.ne.nx_dom)localize=.true.
<       if(Ny.ne.ny_dom)localize=.true.
<       if(La1.ne.La1_dom)localize=.true.
<       if(Lo1.ne.Lo1_dom)localize=.true.
<       if(dx.ne.grid_spacing_dom_m)localize=.true.
< 
<       RETURN
<       end
< c
< c --------------------------------------------------------
< c
<       subroutine rd_static_attr(cstaticdir,nest
<      .,cgrid_fname,nx,ny,dx,dy,La1,Lo1,Latin1
<      .,Latin2,LoV,istatus)
< 
< c routine performs the following:
< c 1. tests if the static file has been generated for this domain or nest?
< c 2. tests if the static variables consistent with the current namelist specs
< c
< c if 1 or 2 is false then we need to either localize or re-localize
< c the domain and "localize" is returned indicating such (true or false).
< c
<       implicit  none
< 
<       integer   nest
<       integer   Nx,Ny
<       integer   istatus
< 
<       character staticfile*200
<       character cstaticdir*200
<       character cgrid_fname*10
<       character cnest*2
< 
<       logical   localize
<       logical   static_exists
< 
<       real      Dx,Dy,La1,Lo1,LoV,Latin1,Latin2
<       real      grid_spacing_dom_m,grid_spacing_m
< 
<       istatus = 1
< 
<       staticfile=TRIM(cstaticdir)//'static.'//cgrid_fname
<       if(nest.gt.1)then
<          write(cnest,'(i2.2)')nest
<          staticfile=TRIM(staticfile)//'_d'//cnest
<       endif
< 
<       INQUIRE(FILE=staticfile, EXIST=static_exists)
<       IF (static_exists) THEN
< 
<         call rd_static_attr_sub(staticfile, Nx, Ny
<      .,La1, Latin1, Latin2, Lo1, LoV, Dx, Dy
<      .,istatus)
< 
<         IF (istatus .NE. 1) THEN
<            print '(A,I5)', ' Error reading WRF static file: ',istatus
<            return
<         END IF
< 
<         IF (LoV .LT. -180.) LoV = LoV + 360.
<         IF (LoV .GT. 180.) LoV = LoV - 360.
<         IF (Lo1 .LT. -180.) Lo1 = Lo1 + 360.
<         IF (Lo1 .GT. 180.) Lo1 = Lo1 - 360.
< 
<       ELSE
< 
<         PRINT '(A)', 'Static file not found: ', staticfile
<         istatus = 0
< 
<       ENDIF
< 
<       return
<       end
done gridgen_utils.f ======================
COMPARE proc_geodat.f
0a1,62
> !1.     given lat/lon points of analysis domain
> !2.     determine the domain subset corner points (lat0 and lon0
> !       lat0=lat(1,1)
> !         where real*4 lat0,lon0,dlat,dlon     !SW corner lat, lon, lat, lon spacing
> !3.     adjust for staggering is done via lat/lon 
> !
> !4.     use latlon_2_llij(nx*ny,lat,lon,gri,grj)
> !         with common /llgrid/ni,nj,nk,lat0,lon0,dlat,dlon
> !         set appropriately for the raw (30s or whatever) data
> !         ni/nj = num raw data gridpoints covering the domain
> !5.     use read_usgs_geog (formerly read_usgs_veg.F)
> !         to extract that part of the full-sized raw data
> !         relevant to the domain.  This routine will need
> !         to be written (starting with read_usgs_veg.F)
> !6.     need function that determines the starting/ending i/j
> !         values for the raw data surrounding the domain lat/lon point.
> !
> !         grid_ratio=(resolution input)/(resolution output)
> 
> !         function istartg(r_grid_ratio,gri)
> !            is = nint(gri-(1./r_grid_ratio) * 0.5)
> !            return
> !            end
> 
> !         function iendg(r_grid_ratio,gri)
> !            ie = int(gri+(1./r_grid_ratio) * 0.5)
> !            return
> !            end
> 
> !         function jstartg(r_grid_ratio,grj)
> !            js = nint(grj-(1./r_grid_ratio) * 0.5)
> !            return
> !            end
> 
> !         function jendg(r_grid_ratio,grj)
> !            je = int(grj+(1./r_grid_ratio) * 0.5)
> !            return
> !            end
> 
> !7.     now do this:
> !
> !       do j=1,ny
> !          do i=1,nx
> !             
> !             is=istartg(r_grid_ratio,gri(i,j)) 
> !             ie=iendg(r_grid_ratio,(gri(i,j))
> !             js=jstartg(r_grid_ratio,gri(i,j))
> !             je=jendg(r_grid-ratio,gri(i,j))
> !             do jd=js,je
> !                do id=is,ie
> !                   dindex=rawdata(id,jd)
> !                   dcnt(i,j,dindex)=dcnt(i,j,dindex)+1
> !                   totnum(i,j)=totnum(i,j)+1 
> !                enddo
> !             enddo
> 
> !          enddo
> !       enddo
> 
> c-----------------------------------------------
> 
> 
11,14d72
< ! J. Smart (NOAA/FSL) : 2003
< !                    Further refined processing of tiles to accomodate
< !                    dateline, Greenwich mean, and Poles. Added smoothing
< !                    of temp field (subroutine one_two_one).
23,24c81
<       integer ntn,nt
<       integer bgmodel
---
>       integer ntn
35a93
>       integer itilelon
42a101
>       integer icat
81d139
<       real    sw(2),ne(2)
87,90d144
<       real    minlat
<       real    maxlat
<       real    minlon
<       real    maxlon
103,104c157
<       logical  lgot0
<       logical  lgotE,lgotW
---
>       logical  lgotW,lgotE
107d159
<       logical  lfndtile
173,176c225,228
<       if(ctiletype.eq.'T'.or.
<      1   ctiletype.eq.'A'.or.
<      1   ctiletype.eq.'G'.or.
<      1   ctiletype.eq.'M' )then
---
>       if(ctiletype.eq.'V'.or.ctiletype.eq.'O'.or.
>      1   ctiletype.eq.'U'.or.ctiletype.eq.'T'.or.
>      1   ctiletype.eq.'A'.or.ctiletype.eq.'G'.or.
>      1   ctiletype.eq.'L'.or.ctiletype.eq.'M' )then
184,185d235
<       bgmodel=6
< 
202c252
<       maxtiles=(360/iblksizo)*(180/iblksizo)
---
>       maxtiles=(360/iblksizo)+(180/iblksizo)
229,232c279,282
<       minlat = MINVAL(dom_lats)
<       maxlat = MAXVAL(dom_lats)
<       minlon = MINVAL(dom_lons)
<       maxlon = MAXVAL(dom_lons)
---
>       min_lat = MINVAL(dom_lats)
>       max_lat = MAXVAL(dom_lats)
>       min_lon = MINVAL(dom_lons)
>       max_lon = MAXVAL(dom_lons)
237,240c287,290
<       min_lat = max(-89.9999, min(89.9999, minlat - abs(rsoff)))
<       max_lat = max(-89.9999, min(89.9999, maxlat + abs(rsoff)))
<       min_lon = max(-359.9999,min(359.9999,minlon - abs(rwoff)))
<       max_lon = max(-359.9999,min(359.9999,maxlon + abs(rwoff)))
---
>       min_lat = max(-89.9999,min(89.9999,min_lat - abs(rsoff)))
>       max_lat = max(-89.9999,min(89.9999,max_lat + abs(rsoff)))
>       min_lon = max(-359.9999,min(359.9999,min_lon - abs(rwoff)))
>       max_lon = max(-359.9999,min(359.9999,max_lon + abs(rwoff)))
268a319,325
> c preliminary ... change sign of lat/lon as necessary and
> c find SW anchor point.
> 
>       min_lon=360
>       max_lon=-360
>       min_lat=90
>       max_lat=0
271d327
<       lgot0=.false.
276,285d331
<          if(itile_lon(itile).eq.0)lgot0=.true.
<          if(ctile_name_list(itile)(7:7).eq.'W')then
<             lgotW=.true.
<             if(itile_lon(itile).lt.180)then
<                itile_lon(itile)=360-itile_lon(itile)
<             endif
<          endif
<          if(ctile_name_list(itile)(7:7).eq.'E')then
<             lgotE=.true.
<          endif
290,297d335
<       enddo
< 
<       min_lon=360
<       max_lon=-360
<       min_lat=90
<       max_lat=0
< 
<       do itile=1,ntn
300,302d337
< c        if(itile_lon(itile).lt.min_lon)min_lon=itile_lon(itile)
< c        if(itile_lon(itile).gt.max_lon)max_lon=itile_lon(itile)
<       enddo
304,321c339,341
<       tile_s_lat = min_lat
<       tile_n_lat = max_lat
< 
<       if(itilesize_d.lt.180)then
<          if(lgotW.and.lgotE)then
<             if(lgot0)then
<                tile_w_lon = float(maxval(itile_lon(1:ntn)))
<                tile_e_lon = float(minval(itile_lon(1:ntn)))
<                if(abs(tile_w_lon-tile_e_lon).lt.180)then
<                   tile_w_lon=tile_w_lon-360.
<                endif
<             else
<                tile_w_lon = float(minval(itile_lon(1:ntn)))
<                tile_e_lon = float(maxval(itile_lon(1:ntn)))
<             endif
<          else
<             tile_w_lon = float(minval(itile_lon(1:ntn)))
<             tile_e_lon = float(maxval(itile_lon(1:ntn)))
---
>          if(ctile_name_list(itile)(7:7).eq.'W')then
> c           lgotW=.true.
>             itile_lon(itile)=-1.0*itile_lon(itile)
323,328c343,344
<       else
<          tile_s_lat = float(minval(itile_lat(1:ntn)))  !+rsoff
<          tile_n_lat = float(maxval(itile_lat(1:ntn)))  !+rsoff
<          tile_w_lon = float(minval(itile_lon(1:ntn)))  !+rwoff
<          tile_e_lon = float(maxval(itile_lon(1:ntn)))  !+rwoff
<       endif
---
>          if(itile_lon(itile).lt.min_lon)min_lon=itile_lon(itile)
>          if(itile_lon(itile).gt.max_lon)max_lon=itile_lon(itile)
329a346,356
> c        if(ctile_name_list(itile)(7:7).eq.'W')then
> c           itile_lon(itile)=360-itile_lon(itile)
> c        endif
> c        if(ctile_name_list(itile)(7:7).eq.'E')lgotE=.true.
>       enddo
> 
>       tile_s_lat = float(minval(itile_lat(1:ntn)))+rsoff
>       tile_n_lat = float(maxval(itile_lat(1:ntn)))+rsoff
>       tile_w_lon = float(minval(itile_lon(1:ntn)))+rwoff
>       tile_e_lon = float(maxval(itile_lon(1:ntn)))+rwoff
> 
335c362,363
<       rlondif=abs(tile_e_lon-tile_w_lon)
---
>       rlondif=tile_e_lon-tile_w_lon
>       if(rlondif.lt.0.0)rlondif=360.+rlondif
351,360c379,381
<       lat0=tile_s_lat+rsoff
<       if(tile_w_lon.gt.180)then
<          lon0=tile_w_lon-360+rwoff 
<       else
<          lon0=tile_w_lon+rwoff
<       endif
<       sw(1)=lat0
<       sw(2)=lon0
<       ne(1)=tile_n_lat+itilesize_d+rsoff
<       ne(2)=tile_e_lon+itilesize_d+rwoff
---
>       lat0=min_lat+rsoff 
> c     if(lgotE .and. lgotW)min_lon=360+min_lon+rwoff
>       lon0=min_lon+rwoff 
363c384
<       print*,'generate supertile for domain'
---
>       print*,'generating supertile for domain'
372,375c393,400
<        IF (icurEW > 180)THEN    ! .and. icurEW < 360) THEN
<           icurEW = 360 - icurEW
<        ENDIF
<        IF (icurEW == 180.and.curEW.ne.'W') THEN
---
>        IF (icurEW > 180) THEN
>           icurEW = icurEW - 180
>           IF (curEW == 'W') THEN
>              curEW = 'E'
>           ELSE
>              curEW = 'W'
>           ENDIF
>        ELSE IF (icurEW == 180) THEN
378d402
< 
416c440
<        elseif( ctiletype.eq.'T')then ! .or. ctiletype.eq.'M')then      ! soiltemp
---
>        elseif( ctiletype.eq.'T' )then      ! soiltemp
437,443c461,462
< c
<        rlondif=float(itile_lon(itile))-tile_w_lon
<        if(rlondif.lt.0)then
<           rlondif=rlondif+360.
<        elseif(rlondif.ge.360)then
<           rlondif=rlondif-360
<        endif
---
>        if(curEW .eq. 'W') icurEW=-1.0*icurEW
>        itilex=nint(1.0+(float(icurEW)-tile_w_lon)/float(itilesize_d))
445,446d463
<        itilex=nint(1.0+(rlondif/float(itilesize_d)))
< 
450,451c467,468
<        print*,'itilex/itiley ',itilex,itiley
<        print*,'itx/ity ',itx,ity
---
> c      print*,'itilex/itiley ',itilex,itiley
> c      print*,'itx/ity ',itx,ity
469,475d485
<       print*,'nxst= ',nxst
<       print*,'nyst= ',nyst
<       print*,'dlat= ',dlat
<       print*,'dlon= ',dlon
<       print*,'lat0= ',lat0
<       print*,'lon0= ',lon0
< 
477c487
<      .dom_lats_in,dom_lons_in,grx,gry,bgmodel,'     ')
---
>      .dom_lats_in,dom_lons_in,grx,gry,1,'     ')
503,507c513
<          if(icnta.gt.0)then
<             amean(k)=asum(k)/icnta
<          else
<             amean(k)=r_missing_data
<          endif
---
>          amean(k)=asum(k)/icnta
509a516,520
> c     print*,'is/ie,js/je ',is,ie,js,je
> c     if(ctiletype.eq.'A')then
> c        write(10)data_proc(is:ie,js:je,1)
> c     endif
> 
511a523
>       method  = 2
513d524
<       nt=1
519,520d529
<          method  = 1
<          geodat(:,:,:)=r_missing_data
526,527d534
<          method  = 2
<          where(amean .eq. r_missing_data)amean=18.
533,534d539
<          method  = 2
<          where(amean .eq. r_missing_data)amean=65.
535a541,545
>       elseif(ctiletype.eq.'L')then
> 
>          min_val = 0.0
>          max_val = 1.0
> 
538c548
<          min_val = 1.0
---
>          min_val = 0.0
540,541d549
<          method  = 1
<          geodat(:,:,:)=65.
549c557,572
<          call interpolate_masked_val(nxst, nyst
---
>          if(.false.)then
>             call hinterp_field(nxst,nyst,nx_dom,ny_dom,1
>      .,grx,gry,data_proc(1,1,ii),geodat(1,1,ii),1)
>             geodat(:,:,ii)=geodat(:,:,ii)/100.
> 
>          elseif(ctiletype.eq.'L')then
> 
>             call hinterp_field(nxst,nyst,nx_dom,ny_dom,1
>      .,grx,gry,data_proc(1,1,ii),geodat(1,1,ii),1)
> 
>             call filter_2dx(geodat,nx_dom,ny_dom,1, 0.5)
>             call filter_2dx(geodat,nx_dom,ny_dom,1,-0.5)
> 
>          else
> 
>             call interpolate_masked_val(nxst, nyst
552a576,578
>             if(ctiletype.eq.'A'.or.ctiletype.eq.'M')then
>                geodat(:,:,ii)=geodat(:,:,ii)/100.
>             endif
554,555d579
<          if(ctiletype.eq.'A'.or.ctiletype.eq.'M')then
<             geodat(:,:,ii)=geodat(:,:,ii)/100.
557d580
<          if(ctiletype.eq.'T')then  !.or.ctiletype.eq.'M')then
559,562d581
<             call one_two_one(nx_dom,ny_dom,nt,geodat(1,1,ii))
< 
<          endif
< 
637,692d655
<       return
< 
<       end
< c
< c-----------------------------------------------------------
< c
<       subroutine one_two_one(nx,ny,nt,data)
< 
<       implicit none
<       integer i,j,ii,jj,icnt,n
<       integer nx,ny,nt
<       integer istatus
<       real  r_missing_data
<       real  fact,sum
<       real  data(nx,ny)
<       real, allocatable :: temp(:,:)
< 
<       if(.not. allocated(temp))allocate(temp(nx,ny))
< 
<       call get_r_missing_data(r_missing_data,istatus)
< 
<       temp=data
< 
<       do n=1,nt
< 
<          do j=2,ny-1
<          do i=2,nx-1
< 
<             sum=0.0
<             icnt=0
<             do jj=j-1,j+1
<             do ii=i-1,i+1
< 
<                fact=1
<                if(jj==j.and.ii==i)fact=2
<                if(data(ii,jj).lt.r_missing_data)then
<                   icnt=icnt+fact
<                   sum=sum+data(ii,jj)*fact
<                endif
< 
<             enddo
<             enddo
< 
<             if(icnt.gt.0)then
<                temp(i,j)=sum/float(icnt)
<             else
<                temp(i,j)=data(i,j)
<             endif
< 
<          enddo
<          enddo
< 
<          data=temp
< 
<       enddo
< 
done proc_geodat.f ======================
COMPARE put_laps_static.f
5c5
<      1                          ,c6_maproj,deltax,deltay,nest)
---
>      1                          ,c6_maproj,deltax,deltay)
32d31
<       character*2       cnest           !nest domain number
57a57
>       laps_dom_file = c10_grid_fname(1:len_fname)
62,69c62,67
<       laps_dom_file = c10_grid_fname(1:len_fname)
< 
<       if(c10_grid_fname(1:lf).eq.'wrfsi')then
<          avgelem=9
<          if(nest.gt.1)then
<             write(cnest,'(i2.2)')nest
<             laps_dom_file = laps_dom_file//'_d'//cnest
<          endif
---
>       if(c10_grid_fname(1:lf).eq.'wrfsi') then
> 	IF (c6_maproj .ne. 'rotlat') THEN
> 		avgelem=9
> 	ELSE
> 		avgelem=18 ! actually AVC, not AVG.
> 	ENDIF
101a100
> 	write(6,*) 'call wrt_laps_static'
done put_laps_static.f ======================
COMPARE rd_static_attr_sub.f
diff: /home/lapb/mpyle/wrfsi_20030217/src/grid/rd_static_attr_sub.f: A file or directory in the path name does not exist.
done rd_static_attr_sub.f ======================
dir=  lib
constant.f cv_asc_i4time.f cv_i4tim_asc_lp.f downcase.f fflxc.f filter_2dx.f filter_filenames.f fname_conversions.f get_dir_length.f get_directory.f get_domain_perimeter.f get_file_names.f get_systime.f get_wrfsi_config.f gridconv.f i4time_asc_gg.f i4time_fname_lp.f lapsgrid.f latlon_to_rlaps.f make_fnam_lp.f map_proj.f move.f rd_laps_static.f read_wrfsi_nl.f timer_fnc.f wfo_fname_conversions.f wrt_laps_static.f ztopsa.f
COMPARE constant.f
done constant.f ======================
COMPARE cv_asc_i4time.f
done cv_asc_i4time.f ======================
COMPARE cv_i4tim_asc_lp.f
done cv_i4tim_asc_lp.f ======================
COMPARE downcase.f
done downcase.f ======================
COMPARE fflxc.f
241a242,245
> 	elseif(c6_maproj .eq. 'rotlat') then
> 
> 	sigma=1.0
> 
done fflxc.f ======================
COMPARE filter_2dx.f
done filter_2dx.f ======================
COMPARE filter_filenames.f
done filter_filenames.f ======================
COMPARE fname_conversions.f
364a365,366
>       elseif(map_proj_name(1:len).eq.'rotlat')then
>          c6_maproj='rotlat'
done fname_conversions.f ======================
COMPARE get_dir_length.f
done get_dir_length.f ======================
COMPARE get_directory.f
done get_directory.f ======================
COMPARE get_domain_perimeter.f
done get_domain_perimeter.f ======================
COMPARE get_file_names.f
done get_file_names.f ======================
COMPARE get_systime.f
done get_systime.f ======================
COMPARE get_wrfsi_config.f
1c1
<       subroutine get_wrfsi_config(nest,istatus)
---
>       subroutine get_wrfsi_config(istatus)
63d62
< 
65,72d63
<            grid_cen_lat_cmn   =moad_known_lat(nest)
<            grid_cen_lon_cmn   =moad_known_lon(nest)
<            nx_l_cmn           =xdim(nest)
<            ny_l_cmn           =ydim(nest)
<            grid_spacing_m_cmn =grid_spacing_wrf_m(nest)
<            i_orig_cmn         =domain_origin_parent_x(nest)
<            j_orig_cmn         =domain_origin_parent_y(nest)
< 
134,137d124
<         if(nest.eq.0)then
<            nest=1    !this would be the MOAD
<         endif
< 
141,146c128,129
<         grid_cen_lat_cmn   =moad_known_lat(nest)
<         grid_cen_lon_cmn   =moad_known_lon(nest)
<         num_domains_cmn    =num_domains
<         i_orig_cmn         =domain_origin_parent_x(nest)
<         j_orig_cmn         =domain_origin_parent_y(nest)
< 
---
>         grid_cen_lat_cmn   =moad_known_lat
>         grid_cen_lon_cmn   =moad_known_lon
157d139
<         path_to_islope = islope
178a161
>         nest = 1
186a170
> 
188a173
> 
done get_wrfsi_config.f ======================
COMPARE gridconv.f
697d696
< c     smart/pw li (04/1/03)
698a698,699
> ccc      implicit none
> c
700,701c701
<       integer nx,ny,nz             !No. of LL domain grid points
<  
---
> c
703,708c703,707
<       real   lli(np),llj(np)       !Lat-lon grid i,j
<       real   diff
< c     real   sw(2),ne(2)
< c     common /llgrid/nx,ny,sw,ne,cgrddef
< 
<       real   lat0,lon0,dlat,dlon     !SW corner lat, lon, lat, lon spacing
---
>       real   lli(np),llj(np),      !Lat-lon grid i,j
>      .       diff
> c
>       integer nx,ny,nz             !No. of LL domain grid points
>       real*4 lat0,lon0,dlat,dlon     !SW corner lat, lon, lat, lon spacing
717,725d715
< 
< c pw li code for integration later on.
< c     dlat=(ne(2)-sw(2))/(nx-1)
< c     dlon=(ne(1)-sw(1))/(ny-1)
< c     do i=1,n
< c        lli(i)=(glon(i)-sw(2))/dlon + 1.
< c        llj(i)=(glat(i)-sw(1))/dlat + 1.
< c     enddo
<  
825,828c815,816
<       diff=xmax-xmin
<       if(diff.lt.0.000001)then
<          xmax=dabs(xmin)
<       endif
---
>       if(xmax.eq.xmin)xmax=abs(xmin)
>       if(ymax.eq.ymin)ymax=abs(ymin)
868,871c856,857
<       diff=xmax-xmin
<       if(diff.lt.0.000001)then
<          xmax=dabs(xmin)
<       endif
---
>       if(xmax.eq.xmin)xmax=abs(xmin)
>       if(ymax.eq.ymin)ymax=abs(ymin)
917,918c903,904
<       ymax=R*(log(tan(PI/4.+0.5*(ne(1)-rlatc)*deg2rad)))
<       ymin=R*(log(tan(PI/4.+0.5*(sw(1)-rlatc)*deg2rad)))
---
> c     ymax=R*(1./tanh(sind(ne(1))))
>       ymax=R*(log(tan(PI/4.+(ne(1)-rlatc)*deg2rad)))
920a907,910
> 
>       ymin=R*(1./tanh(sind(sw(1)-rlatc)))
>       ymin=R*(log(tan(PI/4.+(sw(1)-rlatc)*deg2rad)))
> 
932c922,923
<          y=R*(log(tan(PI/4.+0.5*(dlat*deg2rad))))
---
>          y=R*(1./tanh(sind(rlat(i))))
>          y=R*(log(tan(45.+0.5*rlat(i)*deg2rad)))
942,970d932
<       subroutine latlon_2_npij(n,rlat,rlon,ri,rj)
< c
< c this is conversion of lat/lon to ri/rj in a
< c lat-lon grid. Very similar to latlon_2_llij.
< c
<       integer n,i
<       integer nx,ny
<       real*4 rlat(n),rlon(n),ri(n),rj(n)
<       real*4 sw(2),ne(2)
<       real*4 dx,dy
<       common /npgrid/nx,ny,sw,ne
< 
< c     print *, ' Inside latlon_2_npij'
< c     print *, nx,ny,dx,dy
< c     print *, sw(1), sw(2),ne(1),ne(2)
< 
<       dx=(ne(2)-sw(2))/(nx-1)
<       dy=(ne(1)-sw(1))/(ny-1)
<  
<       do i=1,n
<         ri(i)=(rlon(i)-sw(2))/dx + 1.
<         rj(i)=(rlat(i)-sw(1))/dy + 1.
<       enddo
<  
<       return
<       end
< 
< c -------------------------------------------------------
< 
1002c964
<       real*4    lat1_lc,lat2_lc,lon0_lc,sw_lc(2),ne_lc(2)
---
>       real*4    lat1_lc,lat2_lc,lon0_lc,sw(2),ne(2)
1004c966
<      &,lon0_lc,sw_lc,ne_lc
---
>      &,lon0_lc,sw,ne
1011,1018d972
< c
< c *** Common block variables for polar stereographic grid.
< c
<       integer nx_ps,ny_ps,nz_ps    !No. of PS domain grid points
<       real*4 lat0_ps,lon0_ps,rota  !Pol ste. std lat, lon and rotation
<      .      ,sw_ps(2),ne_ps(2)     !SW lat, lon, NE lat, lon
<       common /psgrid/nx_ps,ny_ps,nz_ps,lat0_ps,lon0_ps
<      .              ,rota,sw_ps,ne_ps
1020,1022d973
<       integer nx_np,ny_np
<       real    sw_np(2),ne_np(2)
<       common /npgrid/nx_np,ny_np,sw_np,ne_np
1024d974
< 
1032,1035c982,986
<          sw_lc(1)=sw1
<          sw_lc(2)=sw2
<          ne_lc(1)=ne1
<          ne_lc(2)=ne2
---
>          sw(1)=sw1
>          sw(2)=sw2
>          ne(1)=ne1
>          ne(2)=ne2
>          return
1042,1043c993,994
<          lat0_ll=Lat0
<          lon0_ll=Lon0
---
>          lat0_ll=lat0
>          lon0_ll=lon0
1046a998
>          return
1058a1011,1013
> c        dx=dlon*111100.
> c        dy=dlat*111100.
>          return 
1061,1081d1015
<       if(gproj.eq.'PS')then
<          nx_ps=nxbg
<          ny_ps=nybg
<          nz_ps=nzbg
<          lat0_ps=Lat0
<          lon0_ps=Lon0
<          sw_ps(1)=sw1
<          sw_ps(2)=sw2
<          ne_ps(1)=ne1
<          ne_ps(2)=ne2
<       endif
< 
<       if(gproj.eq.'NP')then
<          nx_np=nxbg
<          ny_np=nybg
<          sw_np(1)=sw1
<          sw_np(2)=sw2
<          ne_np(1)=ne1
<          ne_np(2)=ne2
<       endif
< 
1101,1102c1035,1036
<       character*132 cmodel
<       character*2   gproj
---
>       character*(*) cmodel
>       character*2 gproj
1107d1040
<       integer nx,ny
1109d1041
<       real nw(2),se(2),rlatc,rlonc
1116c1048
< c________________________________________________________________________________
---
> c_______________________________________________________________________________
1139,1140d1070
<       elseif (gproj .eq. 'NP') then
<          call latlon_2_npij(nx_laps*ny_laps,lat,lon,grx,gry)
1160,1161c1090,1093
<       tolx=(abs(grxdifsum1)/ny_laps+abs(grxdifsum2)/ny_laps)*0.5
<       toly=(abs(grydifsum1)/nx_laps+abs(grydifsum2)/nx_laps)*0.5
---
>       tolx=(abs(grxdifsum1)/ny_laps
>      .+abs(grxdifsum2)/ny_laps)*0.5
>       toly=(abs(grydifsum1)/nx_laps
>      .+abs(grydifsum2)/nx_laps)*0.5
1163a1096
> 
1173c1106
<                if (grx(i,j) .lt. 1.) grx(i,j)=grx(i,j)+float(nx_bg)
---
>                if (grx(i,j) .lt. 1) grx(i,j)=grx(i,j)+float(nx_bg)
1175c1108
<                if (gry(i,j) .lt. 1.) then
---
>                if (gry(i,j) .lt. 1) then
1178c1111
<                   if (grx(i,j) .lt. 1.) grx(i,j)=grx(i,j)+float(nx_bg)
---
>                   if (grx(i,j) .lt. 1) grx(i,j)=grx(i,j)+float(nx_bg)
1184c1117
<                   if (grx(i,j) .lt. 1.) grx(i,j)=grx(i,j)+float(nx_bg)
---
>                   if (grx(i,j) .lt. 1) grx(i,j)=grx(i,j)+float(nx_bg)
1187,1189d1119
<                if (grx(i,j) .lt. 1) then
<                   grx(i,j)=grx(i,j)+1.
<                endif
1195c1125
<                if (grx(i,j) .lt. 1.) grx(i,j)=grx(i,j)+float(nx_bg)
---
>                if (grx(i,j) .lt. 1) grx(i,j)=grx(i,j)+float(nx_bg)
1197c1127
<                if (gry(i,j) .lt. 1.) then
---
>                if (gry(i,j) .lt. 1) then
1200c1130
<                   if (grx(i,j) .lt. 1.) grx(i,j)=grx(i,j)+float(nx_bg)
---
>                   if (grx(i,j) .lt. 1) grx(i,j)=grx(i,j)+float(nx_bg)
1206c1136
<                   if (grx(i,j) .lt. 1.) grx(i,j)=grx(i,j)+float(nx_bg)
---
>                   if (grx(i,j) .lt. 1) grx(i,j)=grx(i,j)+float(nx_bg)
done gridconv.f ======================
COMPARE i4time_asc_gg.f
done i4time_asc_gg.f ======================
COMPARE i4time_fname_lp.f
done i4time_fname_lp.f ======================
COMPARE lapsgrid.f
368c368
<            call get_wrfsi_config(nest,istatus)
---
>            call get_wrfsi_config(istatus)
410c410
<      1  ,l_compress_radar
---
>      1  ,l_highres
425,426c425,427
<      1  ,c8_project_common,c_raddat_type,c80_description
<      1  ,path_to_topt30s ,path_to_topt10m
---
>      1  ,c8_project_common
>      1  ,c_raddat_type, c80_description, path_to_topt30s
>      1  ,path_to_topt10m, path_to_pctl10m
430c431
<      1  ,path_to_islope,path_to_sst,fdda_model_source_cmn
---
>      1  ,path_to_sst,fdda_model_source_cmn
902d902
<       subroutine get_num_domains(num_domains,istatus)
904,948d903
<       include 'lapsparms.cmn'              ! num_domains_cmn
<       include 'grid_fname.cmn'             ! grid_fnam_common
< 
< !     character*80 grid_fnam
< 
< !     This routine accesses the num_domains variable from the
< !     namelist file via the common block. Note the variable names in the
< !     argument list may be different in the calling routine
< 
<       call get_laps_config(grid_fnam_common,istatus)
<       if(istatus .ne. 1)then
<           write(6,*)' ERROR, get_laps_config not successfully called'
<           return
<       endif
< 
<       num_domains = num_domains_cmn
< 
<       istatus = 1
<       return
<       end
< 
<       subroutine get_domain_origin_parent(i_orig,j_orig,istatus)
< 
<       include 'lapsparms.cmn'              ! i_orig_cmn,j_orig_cmn
<       include 'grid_fname.cmn'             ! grid_fnam_common
< 
< !     character*80 grid_fnam
< 
< !     This routine accesses the I&J_DOMAIN_ORIGIN_PARENT variables from the
< !     namelist file via the common block. Note the variable names in the
< !     argument list may be different in the calling routine
< 
<       call get_laps_config(grid_fnam_common,istatus)
<       if(istatus .ne. 1)then
<           write(6,*)' ERROR, get_laps_config not successfully called'
<           return
<       endif
< 
<       i_orig=i_orig_cmn
<       j_orig=j_orig_cmn
< 
<       istatus = 1
<       return
<       end
< 
1206,1230d1160
<       subroutine get_l_compress_radar(l_compress_radar_ret, istatus)
< 
<       include 'lapsparms.cmn' ! l_compress_radar
<       include 'grid_fname.cmn'! grid_fnam_common
< 
<       logical l_compress_radar_ret
< 
< !     This routine accesses the 'l_compress_radar' variable from the
< !     .parms file via the common block. Note the variable names in the
< !     argument list may be different in the calling routine
< 
<       call get_laps_config(grid_fnam_common,istatus)
<       if(istatus .ne. 1)then
<          write(6,*)' ERROR, get_laps_config not successfully called'       
<          return
<       endif
< 
<       l_compress_radar_ret = l_compress_radar
< 
<       istatus = 1
<       return
<       end
< c
< c-----------------------------------------------------------------------
< c
1384a1315,1336
>       subroutine get_path_to_pctl_10m(path_to_pctl_10m,istatus)
> 
>       include 'lapsparms.cmn' ! path_to_pctl_10m
>       include 'grid_fname.cmn'! grid_fnam_common
> 
> !     This routine accesses the path_to_pctl_10m  variable from the
> !     .parms file via the common block. 
> 
>       character*200 path_to_pctl_10m
> 
>       call get_laps_config(grid_fnam_common,istatus)
> 
>       if(istatus .ne. 1)then
>           write(6,*)' ERROR, get_laps_config not successfully called'
>           return
>       endif
> 
>       path_to_pctl_10m =  path_to_pctl10m
> 
>       return
>       end
> 
1542c1494
<       pathtomaxsnoalb = path_to_maxsnoalb
---
>       pathtomaxsnoalb =  path_to_maxsnoalb
1546,1568d1497
< c
< c---------------------------------------------------------------
< 	subroutine get_path_to_islope(pathtoislope,istatus)
< 
< 	include 'lapsparms.cmn' ! path_to_islope
< 	include 'grid_fname.cmn'! grid_fnam_common
< 
< !	This routine accesses the path_to_islope variable from the
< !	.parms file via the common block.
< 
< 	character*200 pathtoislope
< 
< 	call get_laps_config(grid_fnam_common,istatus)
< 
< 	if(istatus .ne. 1)then
< 	  write(6,*)' ERROR, get_laps_config not successfully called'
< 	  return
< 	endif
< 
<         pathtoislope = path_to_islope
< 
< 	return
< 	end
done lapsgrid.f ======================
COMPARE latlon_to_rlaps.f
41d40
< !       2002            Dan Birkenheuer - added input test for valid lats
56,61d54
<         if (abs(rlat) > 90.000) then
<            write(6,*) 'rejecting invalid latitude ',rlat
<            istatus = -1
<            return
<         endif
< 
204c197
<             call get_grid_center(cenlat_dum,cenlon,istatus)
---
>             call get_standard_longitude(cenlon,istatus)
206c199
<                print*,'Error returned from get_grid_center'
---
>                print*,'Error returned from get_standard_longitude'
213a207,216
> !mp
> 	elseif(c6_maproj .eq. 'rotlat')then
> 
> !!	do nothing for now...test
> !!	do nothing for now...test
> !!	do nothing for now...test
> !!	do nothing for now...test
> 
> 	
> 
362c365
<             call get_grid_center(cenlat_dum,cenlon,istatus)
---
>             call get_standard_longitude(cenlon,istatus)
364c367
<                print*,'Error returned from get_grid_center'
---
>                print*,'Error returned from get_standard_longitude'
584a588,596
> 
> 
> !mp	
> 	elseif(c6_maproj .eq. 'rotlat') then ! rotated lat-lon
> 
> 	write(6,*) 'shouldnt be here for rotlat'
> 	projrot_laps=0.
> 
> 
917c929
<               rconst =       tand(colat1)    / tand(colat1/2.)**n
---
>               rconst =  s *  tand(colat1)    / tand(colat1/2.)**n
967a980
> 	if (c6_maproj .ne. 'rotlat') then
968a982,984
> 	else
> 	sigma=1.0
> 	endif
done latlon_to_rlaps.f ======================
COMPARE make_fnam_lp.f
done make_fnam_lp.f ======================
COMPARE map_proj.f
done map_proj.f ======================
COMPARE move.f
done move.f ======================
COMPARE rd_laps_static.f
done rd_laps_static.f ======================
COMPARE read_wrfsi_nl.f
192c192
<      &,islope, sstemp
---
>      &,islope_index, sstemp
413d412
<       include 'grid_fname.cmn'
415a415
>       integer  nest
425c425
<       grid_center_lon=moad_known_lon(nest)
---
>       grid_center_lon=moad_known_lon
done read_wrfsi_nl.f ======================
COMPARE timer_fnc.f
done timer_fnc.f ======================
COMPARE wfo_fname_conversions.f
done wfo_fname_conversions.f ======================
COMPARE wrt_laps_static.f
40c40
<       integer        f_len
---
>       integer        f_len, n
86a87,92
> !mp
>       if (map_proj .eq. 'rotlat') then
> 	 map_projection = 'rotated lat-lon'
> 	write(6,*) 'set map_proj to rotlat'
> 	endif
> !mp
91a98,100
> 	write(6,*) 'n_grids= ', n_grids
> 
> 	write(6,*) 'call write_cdf_static '
99a109,110
> 
> 	write(6,*) 'returned from write_cdf_static with status: ', status
done wrt_laps_static.f ======================
COMPARE ztopsa.f
done ztopsa.f ======================