diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90
index 70d60ac9e..a5ba74ffc 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90
@@ -29,7 +29,7 @@ program mkOverlaySimple
 
   REAL_,   parameter     :: PI        = RASTER_PI
 
-  integer                :: IARGC
+  integer                :: command_argument_count
   integer                :: nxt, argl, fill
   integer                :: i, j, k, l, ip
   integer                :: STATUS, i1, i2, nvars, rvars
@@ -75,7 +75,7 @@ program mkOverlaySimple
     rstdir='rst/'   ! Write in current dir
     maxtiles=4000000
 
-    I = iargc()
+    I = command_argument_count()
 
     if(I < 2 .or. I > 11) then
        print *, trim(Usage)
@@ -83,7 +83,7 @@ program mkOverlaySimple
     end if
 
     nxt = 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
 
     do while(arg(1:1)=='-')
        opt=arg(2:2)
@@ -92,7 +92,7 @@ program mkOverlaySimple
        if(argl==2) then
           if(scan(opt,'zvh')==0) then
              nxt = nxt + 1
-             call getarg(nxt,arg)
+             call get_command_argument(nxt,arg)
           end if
        else
           arg = arg(3:)
@@ -118,13 +118,13 @@ program mkOverlaySimple
        end select
 
        nxt = nxt + 1
-       call getarg(nxt,arg)
+       call get_command_argument(nxt,arg)
     end do
 
     Grid1 = ARG
 
     nxt = nxt + 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
 
     Grid2 = ARG
 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/ConvertAlb.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/ConvertAlb.F90
index b82b1a14b..360179290 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/ConvertAlb.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/ConvertAlb.F90
@@ -11,7 +11,7 @@
 allocate(albo(nx,ny))
 allocate(alb_out(nx,ny))
 
-call getarg(1,ifile)
+call get_command_argument(1,ifile)
 
 ofile = "data/Attic/foo"!trim(ifile)//'int1'
 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/FillMomGrid.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/FillMomGrid.F90
index 7f3502fbf..3b596641f 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/FillMomGrid.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/FillMomGrid.F90
@@ -27,7 +27,7 @@ program FillMomGrid
 
   REAL(KIND=REAL64),   parameter     :: PI        = MAPL_PI
 
-  integer                :: IARGC
+  integer                :: command_argument_count
   integer                :: nxt, argl, fill
   integer                :: i, j, k, l, ip
   integer                :: STATUS, i1, i2, nvars, rvars
@@ -94,7 +94,6 @@ program FillMomGrid
     end if
 
     nxt = 1
-    !call getarg(nxt,arg)
     call get_command_argument(nxt,arg)
 
     do while(arg(1:1)=='-')
@@ -104,7 +103,6 @@ program FillMomGrid
        if(argl==2) then
           if(scan(opt,'zvh')==0) then
              nxt = nxt + 1
-             !call getarg(nxt,arg)
              call get_command_argument(nxt,arg)
           end if
        else
@@ -133,20 +131,17 @@ program FillMomGrid
        end select
 
        nxt = nxt + 1
-       !call getarg(nxt,arg)
        call get_command_argument(nxt,arg)
     end do
 
     Grid1 = ARG
 
     nxt = nxt + 1
-    !call getarg(nxt,arg)
     call get_command_argument(nxt,arg)
 
     Grid2 = ARG
 
     nxt = nxt + 1
-    !call getarg(nxt,arg)
     call get_command_argument(nxt,arg)
 
     GridFile = arg
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/chk_clsm_params.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/chk_clsm_params.F90
index 9eb0f6eea..9a64355ac 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/chk_clsm_params.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/chk_clsm_params.F90
@@ -42,11 +42,11 @@ PROGRAM chk_clsm_params
 character*300 :: input_data
 
 read_file = .true.
-call system ('mkdir -p idl_out')
-call system ('cd bin/ ; /bin/cp ../src/plot_curves.csh . ; cd .. ; /bin/cp bin/plot_curves.csh idl_out/. ; chmod 755 idl_out/plot_curves.csh')
+call execute_command_line ('mkdir -p idl_out')
+call execute_command_line ('cd bin/ ; /bin/cp ../src/plot_curves.csh . ; cd .. ; /bin/cp bin/plot_curves.csh idl_out/. ; chmod 755 idl_out/plot_curves.csh')
 
 
-n = iargc()
+n = command_argument_count()
 
 if(n < 3) then
        print *, "Usage : chk_clsm_params -s Y(N) -m MaskFile"
@@ -56,13 +56,13 @@ PROGRAM chk_clsm_params
 end if
 
 n = 1
-    call getarg(n,arg)
+    call get_command_argument(n,arg)
     do while(arg(1:1)=='-')
        opt=arg(2:2)
        if(len(trim(arg))==2) then
           if(scan(opt,'zvh')==0) then
              n = n + 1
-             call getarg(n,arg)
+             call get_command_argument(n,arg)
           endif
        else
           arg = arg(3:)
@@ -74,7 +74,7 @@ PROGRAM chk_clsm_params
          MaskFile  = trim(arg)
        end select
        n = n + 1
-       call getarg(n,arg)
+       call get_command_argument(n,arg)
     end do
 
 if((single_tile == 'Y').or.(single_tile == 'y')) read_file = .false.
@@ -195,6 +195,6 @@ PROGRAM chk_clsm_params
 close (10,status = 'keep')
 close (11,status = 'keep')
 
-call system ('cd idl_out/ ; ./plot_curves.csh ; cd ..')
+call execute_command_line ('cd idl_out/ ; ./plot_curves.csh ; cd ..')
 
 END PROGRAM chk_clsm_params
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/create_vegdyn_ndvi.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/create_vegdyn_ndvi.F90
index e68138f70..71f0d88c6 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/create_vegdyn_ndvi.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/create_vegdyn_ndvi.F90
@@ -21,7 +21,7 @@ PROGRAM create_vegdyn_ndvi
 
   implicit none
   
-  integer            :: NTILES, N, I, k, iargc, JPLH
+  integer            :: NTILES, N, I, k, command_argument_count, JPLH
   integer, parameter :: nx = 8640, ny = 4320
 
   character*400 :: BCSDIR, OUTDIR, GFILE, IMxJM, arg(5)
@@ -30,18 +30,18 @@ PROGRAM create_vegdyn_ndvi
   real, pointer, dimension (:)  :: z2, z0, ityp
   include 'netcdf.inc'	
 
-  I = iargc()
+  I = command_argument_count()
   IMxJM =''
   GFILE =''
 
-  if(iargc() /= 5) then
-     print *, "Wrong Number of arguments: ", iargc()
+  if(command_argument_count() /= 5) then
+     print *, "Wrong Number of arguments: ", command_argument_count()
      print *, "Usage : ./create_vegdyn_ndvi BCSDIR GFILE, IMxJM JPLH OUTDIR"
      stop
   endif
   
   do n=1,5
-     call getarg(n,arg(n))
+     call get_command_argument(n,arg(n))
   enddo
   
   read(arg(1),'(a)') BCSDIR
@@ -53,8 +53,8 @@ PROGRAM create_vegdyn_ndvi
   ! create dirs/links
   ! -----------------
 
-  call system('mkdir -p data ; cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')
-  call system('mkdir -p '//trim(OUTDIR)) 
+  call execute_command_line('mkdir -p data ; cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')
+  call execute_command_line('mkdir -p '//trim(OUTDIR)) 
 
   open (10,file = trim(BCSDIR)//'/clsm/catchment.def', &
        form= 'formatted', action = 'read', status = 'old')
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_during_day.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_during_day.f90
index b24cb32c7..904693112 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_during_day.f90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_during_day.f90
@@ -90,7 +90,7 @@ PROGRAM loss_during_day
     mult_jobs = .false.
     jseg = ntyps
     job=1
-    I = iargc()
+    I = command_argument_count()
 
     if(I < 1 ) then
        print *, "Job Segment is not specified: ", i
@@ -99,9 +99,9 @@ PROGRAM loss_during_day
     end if
 
     if (i>0) then 
-       call getarg(1,arg)
+       call get_command_argument(1,arg)
        if ( trim(arg) == '-js' ) then
-          call getarg(2,arg)
+          call get_command_argument(2,arg)
           read(arg,'(i2)') job
           mult_jobs = .true.
           jseg =32
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs
index a174aff89..bf9065677 100755
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs
@@ -51,11 +51,11 @@ if ( "$1" == "-noland" ) set SKIPLAND = YES
 
 if ( $HELPMODE != YES ) then
 
-   setenv NCPUS `/usr/bin/lscpu | grep '^CPU(s)' | cut -d ':' -f2 | head -1 `
-   @ NCPUS = $NCPUS / 4
-   @ NCPUS = $NCPUS * 3
+   #setenv NCPUS `/usr/bin/lscpu | grep '^CPU(s)' | cut -d ':' -f2 | head -1 `
+   #@ NCPUS = $NCPUS / 4
+   #@ NCPUS = $NCPUS * 3
    set NCPUS = 20
-   setenv OMP_NUM_THREADS 1
+   setenv OMP_NUM_THREADS  1
 
 endif
    
@@ -93,29 +93,33 @@ set lbcsv = NL3
 
 if ( $HELPMODE != YES ) then
    echo "${CR}"
-   echo "Boundary Conditions (BCs) Package:"    
+   echo "----------------------------------------------------------------------------------------------"
+   echo "                  Boundary Conditions (BCs) Package:"    
+   echo "----------------------------------------------------------------------------------------------"
+   echo " "
+   echo "Enter 3-character alphanumeric code for land BCs version:"
+   echo "(Select only one.)"
    echo " "
-   echo "Enter 3-character alphanumeric code for land BCs version (select only one):"
 else
    echo "Options for land BCs version:"
 endif
-echo "   ${C2}F25 : Fortuna-2_5"
-echo "   ${C2}GM4 : Ganymed-4_0  (/discover/nobackup/ltakacs/bcs/Ganymed-4_0/)"
-echo "   ${C2}ICA : Icarus       (/discover/nobackup/ltakacs/bcs/Icarus/)"
-echo "   ${C2}NL3 : Icarus-NLv3  (/discover/nobackup/ltakacs/bcs/Icarus-NLv3/)"
-echo "   ${C2}NL4 : NLv4 [SMAP]  (/discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)"
-echo "   ${C2}NL5 : NLv5 [SMAP]"
+echo "   ${C2}F25 : Fortuna-2_5   (archived${CR}${C1}*${CR}${C2}: n/a)${CR}"
+echo "   ${C2}GM4 : Ganymed-4_0   (archived${CR}${C1}*${CR}${C2}: /discover/nobackup/ltakacs/bcs/Ganymed-4_0/)${CR}"
+echo "   ${C2}ICA : Icarus        (archived${CR}${C1}*${CR}${C2}: /discover/nobackup/ltakacs/bcs/Icarus/)${CR}"
+echo "   ${C2}NL3 : Icarus-NLv3   (archived${CR}${C1}*${CR}${C2}: /discover/nobackup/ltakacs/bcs/Icarus-NLv3/)${CR}"
+echo "   ${C2}NL4 : NLv4 [SMAPL4] (archived${CR}${C1}*${CR}${C2}: /discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)${CR}"
+echo "   ${C2}NL5 : NLv5 [SMAPL4] (archived${CR}${C1}*${CR}${C2}: /discover/nobackup/projects/gmao/smap/SMAP_L4/L4_SM/bcs/CLSM_params/Icarus-NLv5_EASE/)${CR}"
 echo "   ${C2}DEV : Development version${CR}"
-echo "                                            "
-if ( $HELPMODE != YES ) then
-   echo "              NOTE: Due to compiler differences, code improvements and bug fixes that"
-   echo "                    have taken place since the above archived BCs were created, some parameter"
-   echo "                    files produced by current source code may differ from those in the archived BCs."
-   echo "                    Nevertheless, the impact of these differences on science is insignificant and"
-   echo "                    the parameter files produced by current source code is considered to be"
-   echo "                    scientifically equivalent to the corresponding archived BCs"
-   echo "                                      "
-   echo "              OR press ENTER to select $lbcsv (current default).${CR}"
+echo "  "
+if ( $HELPMODE != YES ) then        
+   echo "                         ${C1}       *BCs produced by this code will differ from BCs in archived directories\!\!\! ${CR}"
+   echo "                                 These differences are caused by compiler changes, code improvements and bug"
+   echo "                                 fixes that were implemented since the archived BCs in the above-mentioned"
+   echo "                                 directories were originally created.  The impact of these differences on"
+   echo "                                 science is insignificant, and the parameter files produced by current"
+   echo "                                 code are scientifically equivalent to the corresponding archived BCs."
+   echo "                                 "
+   echo " OR press ENTER to select $lbcsv (current default).${CR}"
    echo " "
 endif
 
@@ -133,7 +137,7 @@ if ( $HELPMODE != YES ) then
        $dummy == 'NL4'  | \
        $dummy == 'NL5'  | \
        $dummy == 'DEV') then
-      set lbcsv = $dummy
+      set lbcsv        = $dummy
    else if ( $dummy == '' ) then
       echo $lbcsv
    else
@@ -141,6 +145,7 @@ if ( $HELPMODE != YES ) then
       echo " ${C1} Invalid choice, try again:${CR}"
       goto LBCSV
    endif
+
 endif
 
 #######################################################################
@@ -312,17 +317,34 @@ endif
 # Experiment directory (for BCS output)
 
 # define default
-setenv EXPDIR /discover/nobackup/$USER/BCS_PACKAGE/$lbcsv/
+set EXPDIR = /discover/nobackup/$USER/BCS_PACKAGE/$lbcsv
+
+# append minor version string to default bcs output dir name if output from current code differs from archived bcs
+
+if( $lbcsv == 'F25'  | \
+    $lbcsv == 'GM4'  | \
+    $lbcsv == 'ICA'  | \
+    $lbcsv == 'NL3'  | \
+    $lbcsv == 'NL4'  | \
+    $lbcsv == 'NL5') then
+   set EXPDIR = $EXPDIR'_2/'
+   echo "----------------------------------------------------------------------------------------------"
+   echo " "
+   echo " ${C1} NOTE: $lbcsv BCs produced by this code will differ from BCs in archived directories\!\!\! ${CR}"
+   echo " ${C1}       For clarity, '_[v]' was appended to the default output directory shown below.${CR}"
+else
+   # just append trailing slash
+   set EXPDIR = $EXPDIR'/'
+endif
 
 if ( $HELPMODE != YES ) then
 
    echo " "
-   echo "Enter desired BCS output directory (incl. full path)"
-   echo "   or press ENTER to use the default:"
+   echo "Enter desired BCS output directory (incl. full path) or press ENTER to use the default:"
    echo "   [${C2}${EXPDIR}${CR}]"
    set NEWEXPDIR = $<
 
-   if( $NEWEXPDIR != '' ) setenv EXPDIR $NEWEXPDIR
+   if( $NEWEXPDIR != '' ) set EXPDIR = $NEWEXPDIR
 
    mkdir -p $EXPDIR
    if ($status>0) then
@@ -338,7 +360,7 @@ echo ""
 echo "${C1} Land BCs version:${CR}       ${C2}$lbcsv${CR}"
 echo "${C1} Atmospheric resolution:${CR} ${C2}$HRCODES${CR}"
 echo "${C1} Ocean resolution:${CR}       ${C2}$orslvs${CR}"
-echo "${C1} Experiment directory:${CR}   ${C2}$EXPDIR${CR}"
+echo "${C1} Output directory:${CR}       ${C2}$EXPDIR${CR}"
 
 #######################################################################
 #######################################################################
@@ -1290,13 +1312,14 @@ cd ../
 limit stacksize unlimited
 if ( $EVERSION == EASEv2 ) then
 setenv MASKFILE ${MASKFILE}
-if(${MGRID} == M09 | ${MGRID} == M36) then 
-    bin/mkLandRaster.x -x ${NX} -y ${NY} -v -t ${NT}
-    bin/mkSMAPTilesPara_v2.x -smap_grid ${MGRID} -pfaf_til T
-    bin/CombineRasters.x -f 0 -t 232000000 ${THISGRID} Pfafstetter > /dev/null
-    bin/CombineRasters.x -t 232000000 ${THISGRID} ${THISGRID}-Pfafstetter
-    /bin/mv til/${THISGRID}_${THISGRID}-Pfafstetter.til  til/${THISGRID}_${THISGRID}-Pfafstetter.ind 
-endif
+## This section was used to make Irrigated Tiles
+##if(${MGRID} == M09 | ${MGRID} == M36) then 
+##    bin/mkLandRaster.x -x ${NX} -y ${NY} -v -t ${NT}
+##    bin/mkSMAPTilesPara_v2.x -smap_grid ${MGRID} -pfaf_til T
+##    bin/CombineRasters.x -f 0 -t 232000000 ${THISGRID} Pfafstetter > /dev/null
+##    bin/CombineRasters.x -t 232000000 ${THISGRID} ${THISGRID}-Pfafstetter
+##    /bin/mv til/${THISGRID}_${THISGRID}-Pfafstetter.til  til/${THISGRID}_${THISGRID}-Pfafstetter.ind 
+##endif
 setenv OMP_NUM_THREADS 1
 bin/mkSMAPTilesPara_v2.x -smap_grid ${MGRID} -v $lbcsv 
 setenv OMP_NUM_THREADS ${NCPUS}
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90
index fe035a9e5..c53cab991 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90
@@ -43,7 +43,7 @@ PROGRAM mkCatchParam
   character*4          :: EASE ='    '
   character*2          :: DL ='DC'    
   integer              :: II, JJ, Type
-  integer              :: I, J, iargc, nxt
+  integer              :: I, J, command_argument_count, nxt
   real*8               :: dx, dy, lon0
   logical              :: regrid
   character(len=400), dimension (8) ::  Usage 
@@ -99,8 +99,8 @@ PROGRAM mkCatchParam
 !
 !$OMP ENDPARALLEL
 
-!   call system('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')
-!   call system('cd ..')
+!   call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')
+!   call execute_command_line('cd ..')
 
     USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE                                      "
     USAGE(2) ="     -x: Size of longitude dimension of input raster. DEFAULT: 8640                                     "
@@ -124,7 +124,7 @@ PROGRAM mkCatchParam
        write (log_file,'(a)')' '
     endif
 
-    I = iargc()
+    I = command_argument_count()
     if(I < 1 .or. I > 10) then
        write (log_file,'(a)') "Wrong Number of arguments: ", i
        do j = 1,size(usage)
@@ -134,13 +134,13 @@ PROGRAM mkCatchParam
     end if
 
     nxt = 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
     do while(arg(1:1)=='-')
        opt=arg(2:2)
        if(len(trim(arg))==2) then
           if(scan(opt,'zh')==0) then
              nxt = nxt + 1
-             call getarg(nxt,arg)
+             call get_command_argument(nxt,arg)
           endif
        else
           arg = arg(3:)
@@ -168,10 +168,10 @@ PROGRAM mkCatchParam
           call exit(1)
        end select
        nxt = nxt + 1
-       call getarg(nxt,arg)
+       call get_command_argument(nxt,arg)
     end do
 
-   call getenv ("MASKFILE"        ,MaskFile        )
+   call get_environment_variable ("MASKFILE"        ,MaskFile        )
  
   if(trim(Gridname) == '') then
       write (log_file,'(a)')'Unable to create parameters without til/rst files.... !'
@@ -663,7 +663,7 @@ PROGRAM mkCatchParam
        write (log_file,'(a)')'============================================================'
        write (log_file,'(a)')' '
        
-       !       call system ('chmod 755 bin/create_README.csh ; bin/create_README.csh')
+       !       call execute_command_line ('chmod 755 bin/create_README.csh ; bin/create_README.csh')
     endif
 
     close (log_file,status='keep') 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCubeFVRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCubeFVRaster.F90
index 6a73ac382..02c7c3c90 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCubeFVRaster.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCubeFVRaster.F90
@@ -27,7 +27,7 @@ program mkCubeFVRaster
     character*128        :: ARG, GridName=''
     character*1          :: opt
     integer              :: ncells !  Cells on edges of cubed faces
-    integer              :: I, J, N, status, iargc, nxt
+    integer              :: I, J, N, status, command_argument_count, nxt
     integer              :: js,jv
     real*8               :: dx, dy
     real*8,  allocatable :: xs(:,:), ys(:,:), xv(:,:,:), yv(:,:,:)
@@ -40,7 +40,7 @@ program mkCubeFVRaster
 ! Process Arguments
 !------------------
 
-    I = iargc()
+    I = command_argument_count()
 
     if(I < 1 .or. I > 10) then
        print *, Usage
@@ -48,14 +48,14 @@ program mkCubeFVRaster
     end if
 
     nxt = 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
     do while(arg(1:1)=='-')
 
        opt=arg(2:2)
        if(len(trim(arg))==2) then
           if(scan(opt,'zvh')==0) then
              nxt = nxt + 1
-             call getarg(nxt,arg)
+             call get_command_argument(nxt,arg)
           endif
        else
           arg = arg(3:)
@@ -80,7 +80,7 @@ program mkCubeFVRaster
        end select
 
        nxt = nxt + 1
-       call getarg(nxt,arg)
+       call get_command_argument(nxt,arg)
     end do
 
     read(arg,'(i6)') ncells
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLISTilesPara.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLISTilesPara.F90
index d224fb777..99ad7b2a9 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLISTilesPara.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLISTilesPara.F90
@@ -13,7 +13,7 @@ PROGRAM mkLISTilesPara
   integer, parameter   :: nc_gswp2 = 360, nr_gswp2 = 180, n_gswp2 =15238 
   integer, parameter :: max_pfaf_smap = 100
   character(40) :: arg
-  integer       ::  i, N_args, iargc, status
+  integer       ::  i, N_args, command_argument_count, status
   character*300 :: latlon_vector_file
   integer       :: nc, nr
   character*200 :: gfile
@@ -23,7 +23,7 @@ PROGRAM mkLISTilesPara
   real :: dx, dy
   integer :: ncells, dateline, nc_domain,nr_domain,i_offset,j_offset
   
-  N_args = iargc()
+  N_args = command_argument_count()
 
   if(N_args /= 2) then
      print *,'USAGE : bin/mkLISTilesPara -vfile filename'
@@ -36,11 +36,11 @@ PROGRAM mkLISTilesPara
      
      i = i+1
      
-     call getarg(i,arg)
+     call get_command_argument(i,arg)
      
      if     ( trim(arg) == '-vfile' ) then
         i = i+1
-        call getarg(i,latlon_vector_file)
+        call get_command_argument(i,latlon_vector_file)
 	
      else ! stop for any other arguments
 
@@ -50,14 +50,14 @@ PROGRAM mkLISTilesPara
      endif     
 end do
 
-call system('mkdir -p data/ ; mkdir -p til/ ; mkdir -p rst/ ; mkdir -p clsm/plots')
-call system('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')  
-call system('cd ..')
+call execute_command_line('mkdir -p data/ ; mkdir -p til/ ; mkdir -p rst/ ; mkdir -p clsm/plots')
+call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')  
+call execute_command_line('cd ..')
      
 !   Check for the 10 arc-sec MaskFile
 ! -----------------------------------
 
-call getenv ("MASKFILE"        ,MaskFile        )
+call get_environment_variable ("MASKFILE"        ,MaskFile        )
 
 if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then       
    ! Use new ESA based MaskFile
@@ -87,7 +87,7 @@ PROGRAM mkLISTilesPara
 tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1)
 print *,trim(tmpstring)
 
-call system(tmpstring)
+call execute_command_line(tmpstring)
 
 contains
 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLandRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLandRaster.F90
index 3af0ff1a6..f42a5adc0 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLandRaster.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLandRaster.F90
@@ -22,7 +22,7 @@ Program MakeLandRaster
     integer, parameter     :: nx0 = 8640
     integer, parameter     :: ny0 = 4320
 
-    integer                :: IARGC
+    integer                :: command_argument_count
     integer                :: i,j,k,n, status,ncid, ncid2
     integer                :: ip, nxt
     integer                :: type, maxtiles, nx, ny
@@ -73,8 +73,8 @@ Program MakeLandRaster
     character*128          :: &
     Usage = "mkLandRaster -x nx -y ny -v -h -z -t maxtiles -l LandFile -g GridName"
     include 'netcdf.inc'
-    call system('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')
-    call system('cd ..')
+    call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')
+    call execute_command_line('cd ..')
 
 ! Process Arguments
 !------------------
@@ -90,7 +90,7 @@ Program MakeLandRaster
     InputFile = &
           "data/CATCH/global.cat_id.catch.DL"
      
-    I = iargc()
+    I = command_argument_count()
 
     if(I > 13) then
        print *, "Wrong Number of arguments: ", i
@@ -101,7 +101,7 @@ Program MakeLandRaster
     nxt = 1
 
     do while(nxt<=I)
-       call getarg(nxt,arg)
+       call get_command_argument(nxt,arg)
        if(arg(1:1)/='-') then
           print *, trim(Usage)
           call exit(1)
@@ -110,7 +110,7 @@ Program MakeLandRaster
        if(len(trim(arg))==2) then
           if(scan(opt,'zvh')==0) then
              nxt = nxt + 1
-             call getarg(nxt,arg)
+             call get_command_argument(nxt,arg)
           end if
        else
           arg = arg(3:)
@@ -143,7 +143,7 @@ Program MakeLandRaster
     
 !   Check for the 10 arc-sec MaskFile (SM)
 
-    call getenv ("MASKFILE"        ,MaskFile        )
+    call get_environment_variable ("MASKFILE"        ,MaskFile        )
 
     print *,'Using Mask file : ',trim(MaskFile)
 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLatLonRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLatLonRaster.F90
index c7bb5bce2..6812a6bfe 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLatLonRaster.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLatLonRaster.F90
@@ -59,7 +59,7 @@ program mkLatLonRaster
     integer              :: II, JJ, Type
     logical              :: UseType = .false., DoZip=.false., Verb=.false.
     logical              :: Here = .false.
-    integer              :: I, J, status, iargc, nxt
+    integer              :: I, J, status, command_argument_count, nxt
     real*8               :: dx, dy, lon0
     real*8,  allocatable :: xs(:), ys(:), xv(:,:,:), yv(:,:,:)
     character*128        :: &
@@ -68,7 +68,7 @@ program mkLatLonRaster
 ! Process Arguments
 !------------------
 
-    I = iargc()
+    I = command_argument_count()
 
     if(I < 2 .or. I > 17) then
        print *, "Wrong Number of arguments: ", i
@@ -77,13 +77,13 @@ program mkLatLonRaster
     end if
 
     nxt = 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
     do while(arg(1:1)=='-')
        opt=arg(2:2)
        if(len(trim(arg))==2) then
           if(scan(opt,'zvh')==0) then
              nxt = nxt + 1
-             call getarg(nxt,arg)
+             call get_command_argument(nxt,arg)
           endif
        else
           arg = arg(3:)
@@ -113,13 +113,13 @@ program mkLatLonRaster
           call exit(1)
        end select
        nxt = nxt + 1
-       call getarg(nxt,arg)
+       call get_command_argument(nxt,arg)
     end do
 
     read(arg,'(i5)') ii
 
     nxt = nxt + 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
 
     read(arg,'(i5)') jj
 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMITAquaRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMITAquaRaster.F90
index 18a5a6889..5dd4732e8 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMITAquaRaster.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMITAquaRaster.F90
@@ -12,7 +12,7 @@ program MAIN
       integer, parameter      :: RKIND  = 8
       integer, parameter      :: LL = 1, LR = 2, UR = 3, UL = 4
 
-      integer                 :: iargc
+      integer                 :: command_argument_count
       integer                 :: Nc    = 8640, NR   = 4320
 
       type Ptr2
@@ -44,7 +44,7 @@ program MAIN
 ! Get source grid directory and destination raster file names
 !------------------------------------------------------------
 
-   i = iargc()
+   i = command_argument_count()
 
    if(I < 2 .or. i > 7) then
       print *, "Wrong Number of arguments: ", i
@@ -53,13 +53,13 @@ program MAIN
    end if
    
    nxt = 1
-   call getarg(nxt,arg)
+   call get_command_argument(nxt,arg)
    do while(arg(1:1)=='-')
       opt=arg(2:2)
       if(len(trim(arg))==2) then
          if(scan(opt,'zv')==0) then
             nxt = nxt + 1
-            call getarg(nxt,arg)
+            call get_command_argument(nxt,arg)
          endif
       else
          arg = arg(3:)
@@ -79,7 +79,7 @@ program MAIN
          call exit(1)
       end select
       nxt = nxt + 1
-      call getarg(nxt,arg)
+      call get_command_argument(nxt,arg)
    end do
 
    GridDir = arg
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMOMAquaRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMOMAquaRaster.F90
index fcb9a96e9..7e6fbe4ea 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMOMAquaRaster.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMOMAquaRaster.F90
@@ -19,7 +19,7 @@ program MOMraster
   REAL_,     pointer     :: yvert(:,:,:)          ! Lats of MOM's vertices
   REAL_                  :: xmin, xmax
   integer                :: i, j, nxt,k
-  integer                :: status, iargc
+  integer                :: status, command_argument_count
   character*(128)        :: GridFile
   character*(128)        :: GridName=''
   character*(128)        :: arg
@@ -40,7 +40,7 @@ program MOMraster
 ! Process Arguments
 !------------------
 
-    I = iargc()
+    I = command_argument_count()
 
     if(I < 1 .or. I > 8) then
        print *, "Wrong Number of arguments: ", i
@@ -49,13 +49,13 @@ program MOMraster
     end if
 
     nxt = 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
     do while(arg(1:1)=='-')
        opt=arg(2:2)
        if(len(trim(arg))==2) then
           if(scan(opt,'zvh')==0) then
              nxt = nxt + 1
-             call getarg(nxt,arg)
+             call get_command_argument(nxt,arg)
           endif
        else
           arg = arg(3:)
@@ -78,7 +78,7 @@ program MOMraster
           call exit(1)
        end select
        nxt = nxt + 1
-       call getarg(nxt,arg)
+       call get_command_argument(nxt,arg)
     end do
 
     GridFile = arg
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara.F90
index 2f5f4f7c6..fc62e7718 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara.F90
@@ -52,7 +52,7 @@ PROGRAM mkSMAPTilesPara_v1
       character*100 :: veg_class (12)
       character*5 :: MGRID
       character*100 :: gfile,gtopo30
-      integer :: nc_smap,nr_smap, N_args, iargc 
+      integer :: nc_smap,nr_smap, N_args, command_argument_count 
       real :: EASE_grid_area, CELL_km
       REAL :: dx,dy,d2r,lats,mnx,mxx,mny,mxy,sum1,sum2,jgv, VDUM,pix_area
       character(40) :: arg
@@ -61,7 +61,7 @@ PROGRAM mkSMAPTilesPara_v1
       character*128          :: MaskFile
       include 'netcdf.inc'
 
-      N_args = iargc()
+      N_args = command_argument_count()
 
       if(N_args < 1) then
         print *,'USAGE : bin/mkSMAPTiles_v1 -smap_grid MXX'
@@ -75,11 +75,11 @@ PROGRAM mkSMAPTilesPara_v1
 
          i = i+1
          
-         call getarg(i,arg)
+         call get_command_argument(i,arg)
          
          if     ( trim(arg) == '-smap_grid' ) then
             i = i+1
-            call getarg(i,MGRID)
+            call get_command_argument(i,MGRID)
             
          else ! stop for any other arguments
             
@@ -91,8 +91,8 @@ PROGRAM mkSMAPTilesPara_v1
          
       end do
       
-      call system('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')  
-      call system('cd ..')
+      call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')  
+      call execute_command_line('cd ..')
       
       
       ! Setting SMAP Grid specifications
@@ -127,7 +127,7 @@ PROGRAM mkSMAPTilesPara_v1
       !   Check for the 10 arc-sec MaskFile
       ! -----------------------------------
       
-      call getenv ("MASKFILE"        ,MaskFile        )
+      call get_environment_variable ("MASKFILE"        ,MaskFile        )
 
       print *, 'Using MaskFile ', trim(MaskFile)
 
@@ -690,7 +690,7 @@ PROGRAM mkSMAPTilesPara_v1
       tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1)
       print *,trim(tmpstring)
       
-      call system (tmpstring)   
+      call execute_command_line (tmpstring)   
 
    END PROGRAM mkSMAPTilesPara_v1
 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90
index 79ec28562..70025ec11 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90
@@ -56,7 +56,7 @@ PROGRAM mkSMAPTilesPara_v2
       character*100 :: veg_class (12)
       character*5 :: MGRID
       character*100 :: gfile,gtopo30
-      integer :: nc_smap,nr_smap, N_args, iargc 
+      integer :: nc_smap,nr_smap, N_args, command_argument_count 
       real :: EASE_grid_area, CELL_km
       REAL :: dx,dy,d2r,lats,mnx,mxx,mny,mxy,sum1,sum2,jgv, VDUM,pix_area
       character(40) :: arg, EASElabel 
@@ -68,7 +68,7 @@ PROGRAM mkSMAPTilesPara_v2
       
       include 'netcdf.inc'
 
-      N_args = iargc()
+      N_args = command_argument_count()
 
       if(N_args < 1) then
         print *,'USAGE : bin/mkSMAPTiles -smap_grid MXX'
@@ -82,20 +82,20 @@ PROGRAM mkSMAPTilesPara_v2
 
          i = i+1
          
-         call getarg(i,arg)
+         call get_command_argument(i,arg)
          
          if     ( trim(arg) == '-smap_grid' ) then
             i = i+1
-            call getarg(i,MGRID)
+            call get_command_argument(i,MGRID)
 
          elseif ( trim(arg) == '-pfaf_til' ) then
             i = i+1
-            call getarg(i,PF)
+            call get_command_argument(i,PF)
             if (PF == 'T') pfaf_til = .true.
 
          elseif ( trim(arg) == '-v' ) then
             i = i+1
-            call getarg(i,LBSV)
+            call get_command_argument(i,LBSV)
                         
          else ! stop for any other arguments
             
@@ -107,8 +107,8 @@ PROGRAM mkSMAPTilesPara_v2
          
       end do
       
-      call system('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')  
-      call system('cd ..')
+      call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH')  
+      call execute_command_line('cd ..')
       
       
       ! Setting SMAP Grid specifications
@@ -184,7 +184,7 @@ PROGRAM mkSMAPTilesPara_v2
       !   Check for the 10 arc-sec MaskFile
       ! -----------------------------------
       
-      call getenv ("MASKFILE"        ,MaskFile        )
+      call get_environment_variable ("MASKFILE"        ,MaskFile        )
 
       print *, 'Using MaskFile ', trim(MaskFile)
       
@@ -194,8 +194,9 @@ PROGRAM mkSMAPTilesPara_v2
          nr = 21600
          call mkEASEv2Raster
          
-      else
-         if((trim(MGRID) == 'M09').or.(trim(MGRID) == 'M36'))call write_tilfile 
+      !else
+      !   This section was used to make Irrigated Tiles 
+      !   if((trim(MGRID) == 'M09').or.(trim(MGRID) == 'M36'))call write_tilfile 
       endif
       
       if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then         
@@ -749,8 +750,6 @@ PROGRAM mkSMAPTilesPara_v2
 
       ! CALL CREATE_ROUT_PARA_FILE (NC, NR, trim(gfile), MGRID=MGRID)  
       
-      call system (tmpstring)
-
       ! now run mkCatchParam
       ! --------------------
 
@@ -759,7 +758,7 @@ PROGRAM mkSMAPTilesPara_v2
       tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1)
       print *,trim(tmpstring)
       
-      call system (tmpstring)
+      call execute_command_line (tmpstring)
 
     contains
 
@@ -776,8 +775,8 @@ SUBROUTINE mkEASEv2Raster
         allocate (xs ( nc_smap+1, nr_smap+1))
         allocate (ys ( nc_smap+1, nr_smap+1))
         
-        do  j = 1, nr_smap
-           do i = 1, nc_smap
+        do  j = 1, nr_smap+1
+           do i = 1, nc_smap+1
               x = real(i-1)        -0.5
               y = real(nr_smap - j)+0.5
               call easeV2_inverse(MGRID, x, y, yout, xout)
@@ -785,31 +784,9 @@ SUBROUTINE mkEASEv2Raster
               xs (i,j) = dble(xout)
            end do
         end do
-        
-        do  j = nr_smap + 1, nr_smap + 1
-           do i = nc_smap + 1, nc_smap + 1
-              x = real(i-1)         -0.5
-              y =  -0.5
-              call easeV2_inverse(MGRID, x, y, yout, xout)
-              ys (i,j) = dble(yout)
-              xs (i,j) = dble(xout)        
-           end do
-        end do
+	
+        call  LRRasterize(EASElabel,xs,ys,nc=nc,nr=nr,xmn = xs(1,1), xmx= xs(nc_smap+1, nr_smap+1),  ymn=ys(1,1), ymx = ys(nc_smap+1, nr_smap+1), Here=.false., Verb=.false.)       
 
-        where (ys > 90.)
-           ys = 90.D0
-        endwhere
-        where (ys < -90.)
-           ys = -90.D0
-        endwhere
-        where (xs > 180.)
-           xs = 180.D0
-        endwhere
-        where (xs < -180.)
-           xs = -180.D0
-        endwhere
-        
-        call  LRRasterize(EASElabel,xs,ys,nc=nc,nr=nr,Here=.false., Verb=.false.)       
         stop
       end SUBROUTINE mkEASEv2Raster
 
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mk_runofftbl.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mk_runofftbl.F90
index f69829259..c097f69bc 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mk_runofftbl.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mk_runofftbl.F90
@@ -18,7 +18,7 @@ program Runoff
   character*100          :: fileLL="data/CATCH/Outlet_latlon."
   character*5            :: C_NX, C_NY
 
-  call getarg(1,file)
+  call get_command_argument(1,file)
 
   fileT = "til/"//trim(file)//".til"
   fileR = "rst/"//trim(file)//".rst"
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90
index a75036635..0522879c4 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90
@@ -1139,7 +1139,7 @@ SUBROUTINE create_mapping (nc,nr,nc_data,nr_data,rmap, gfile)
       allocate(iraster(nc_data,nr_data),stat=STATUS); VERIFY_(STATUS)
       call RegridRaster(tile_id,iraster)
       NPLUS = count(iraster >= 1 .and. iraster <= ncatch)
-      allocate (rmap%ij_index(1:nc_data, 1:nr_data))
+      allocate (rmap%ij_index(1:nc_data, 1:nr_data), source = 0)
       allocate (rmap%map (1:NPLUS))     
       rmap%map%NT = 0
       pix_count = 1
@@ -1158,7 +1158,7 @@ SUBROUTINE create_mapping (nc,nr,nc_data,nr_data,rmap, gfile)
 
    else
       NPLUS = count(tile_id >= 1 .and. tile_id <= ncatch)
-      allocate (rmap%ij_index(1:nc_data, 1:nr_data))
+      allocate (rmap%ij_index(1:nc_data, 1:nr_data), source = 0)
       allocate (rmap%map (1:NPLUS))     
       rmap%map%NT = 0
       pix_count   = 1
@@ -1393,15 +1393,15 @@ SUBROUTINE modis_scale_para_high (ease_grid,MA,gfile)
     character*6 :: MA
     CHARACTER*20 :: version,resoln,continent
     integer :: nc_gcm,nr_gcm,nc_ocean,nr_ocean
-    REAL :: latt,lont,fr_gcm,fr_cat,tsteps,zth, slr,tarea
-    INTEGER :: typ,pfs,ig,jg,j_dum,ierr,indx_dum,indr1,indr2,indr3 ,ip2
+    REAL :: tsteps,zth, slr,tarea
+    INTEGER :: typ,j_dum,ierr,indr1,ip2
     character*100 :: path,fname,fout,metpath
     character (*) :: gfile
     integer :: n,maxcat,ip
     integer :: yy,j,month
     integer, allocatable, dimension (:) :: vegcls 
     real, allocatable, dimension (:) :: &
-         modisvf, modisnf,albvf,albnf, lat,lon, &
+         modisvf, modisnf,albvf,albnf, &
          green,lai,lai_before,lai_after,grn_before,grn_after
     real, allocatable, dimension (:) :: &
          calbvf,calbnf, zero_array, one_array, albvr,albnr
@@ -1435,8 +1435,8 @@ SUBROUTINE modis_scale_para_high (ease_grid,MA,gfile)
     close (10,status='keep')
 
     fname=trim(gfile)//'.til'
-
     open (10,file=fname,status='old',action='read',form='formatted')
+
     fname='clsm/mosaic_veg_typs_fracs'
     open (20,file=fname,status='old',action='read',form='formatted')
 
@@ -1451,15 +1451,17 @@ SUBROUTINE modis_scale_para_high (ease_grid,MA,gfile)
     
     do n = 1,ip
       if (ease_grid) then     
-	 read(10,*,IOSTAT=ierr) typ,pfs,lon,lat,ig,jg,fr_gcm
+         read(10,*,IOSTAT=ierr) typ !,pfs,lont,latt,ig,jg,fr_gcm
       else
-      read(10,'(I10,3E20.12,9(2I10,E20.12,I10))',IOSTAT=ierr)     &    
-            typ,tarea,lont,latt,ig,jg,fr_gcm,indx_dum,pfs,j_dum,fr_cat,j_dum
+         !read(10,'(I10,3E20.12,9(2I10,E20.12,I10))',IOSTAT=ierr)     &    
+         !   typ,tarea,lont,latt,ig,jg,fr_gcm,indx_dum,pfs,j_dum,fr_cat,j_dum
+         read(10,*,IOSTAT=ierr) typ
       endif
        if (typ == 100) then
           ip2 = n 
-          read (20,'(i10,i8,2(2x,i3),2(2x,f6.4))')     &
-            indr1,indr1,vegcls(ip2),indr1,fr_gcm,fr_gcm
+          !read (20,'(i10,i8,2(2x,i3),2(2x,f6.4))')     &
+          !  indr1,indr1,vegcls(ip2),indr1,fr_gcm,fr_gcm
+          read (20,*,IOSTAT=ierr) indr1,indr1,vegcls(ip2)
        endif
        if(ierr /= 0)write (*,*)'Problem reading', n, ease_grid
     end do
@@ -1875,6 +1877,7 @@ SUBROUTINE modis_alb_on_tiles_high (nc_data,nr_data,rmap,MA,gfiler)
                    do j = jLL,jLL + nr_10 -1 
                       do i = iLL, iLL + nc_10 -1
                          pix_count = rmap%ij_index(i,j)
+                         if (pix_count ==0) cycle
                          if(net_data1(i-iLL +1 ,j - jLL +1) > 0) then
                             
                             if(rmap%map(pix_count)%nt > 0) then
@@ -2603,6 +2606,7 @@ SUBROUTINE hres_lai_no_gswp (nc_data,nr_data,rmap,gfiler,lai_name, merge)
                      do i = iLL, iLL + nc_10 -1 
                         if(net_data1(i-iLL +1 ,j - jLL +1) /= d_undef) then
                            pix_count = rmap%ij_index(i,j)
+                           if (pix_count ==0) cycle
                            if(rmap%map(pix_count)%nt > 0) then
                               do n = 1, rmap%map(pix_count)%nt
                                  if(vec_lai(rmap%map(pix_count)%tid(n)) == -9999.) vec_lai(rmap%map(pix_count)%tid(n)) = 0.                                 
@@ -2906,6 +2910,7 @@ SUBROUTINE hres_gswp2 (nc_data,nr_data,rmap, gfiler,lai_name,merge)
                  do i = iLL, iLL + nc_10 -1 
                     if(net_data1(i-iLL +1 ,j - jLL +1) /= d_undef) then
                        pix_count = rmap%ij_index(i,j)
+                       if (pix_count == 0) cycle
                        if(rmap%map(pix_count)%nt > 0) then
                           do n = 1, rmap%map(pix_count)%nt
                              if(vec_lai(rmap%map(pix_count)%tid(n)) == -9999.) vec_lai(rmap%map(pix_count)%tid(n)) = 0.                                 
@@ -4194,8 +4199,13 @@ SUBROUTINE soil_para_hwsd (nx,ny,gfiler)
                   j  = 0
                   i1 = n - k
                   i2 = n + k
-                  if((i1 >=     1).and.(soil_class_com (i1) >=1)) j = i1  ! tentatively use "lower" neighbor unless out of range
-                  if((i2 <=maxcat).and.(soil_class_com (i2) >=1)) j = i2  ! "upper" neighbor prevails unless out of range
+                  if(i1 >=     1) then
+                     if (soil_class_com (i1) >=1) j = i1  ! tentatively use "lower" neighbor unless out of range
+                  endif
+
+                  if(1 <= i2 .and. i2 <=maxcat) then
+                     if (soil_class_com (i2) >=1) j = i2  ! "upper" neighbor prevails unless out of range
+                  endif
 
                   if (j > 0) then
                      soil_class_com (n) = soil_class_com (j)
@@ -5915,9 +5925,9 @@ SUBROUTINE open_landparam_nc4_files(N_tile)
       write (time_stamp,'(i4.4,a1,i2.2,a1,i2.2,1x,a2,1x,i2.2,a1,i2.2,a1,i2.2)')      &
            date_time_values(1),'-',date_time_values(2),'-',date_time_values(3),'at', &
            date_time_values(5),':',date_time_values(6),':',date_time_values(7)
-!      call system('setenv    MYNAME `finger $USER | cut -d: -f3 | head -1`')
+!      call execute_command_line('setenv    MYNAME `finger $USER | cut -d: -f3 | head -1`')
 !      call sleep (5)
-      call getenv ("USER"        ,MYNAME        )
+      call get_environment_variable ("USER"        ,MYNAME        )
       status = NF_PUT_ATT_TEXT(NCCatOUTID  , NF_GLOBAL, 'CreatedBy', LEN_TRIM(MYNAME),  trim(MYNAME)      )
       status = NF_PUT_ATT_TEXT(NCCatOUTID  , NF_GLOBAL, 'Date'     , LEN_TRIM(time_stamp),trim(time_stamp))
       status = NF_PUT_ATT_TEXT(NCVegOUTID  , NF_GLOBAL, 'CreatedBy', LEN_TRIM(MYNAME),  trim(MYNAME)      )
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90
index fd381b50c..57092add8 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90
@@ -950,10 +950,10 @@ SUBROUTINE remove_tiny_tiles (                         &
     print *,trim(gname)
     string1 ='til/'//trim(gout)//'-Pfafstetter.til'//' '//&
     'clsm/'//trim(gout)//'-Pfaf.notiny.til'
-    call system ('cp '//trim(string1))
+    call execute_command_line ('cp '//trim(string1))
     string1 ='rst/'//trim(gout)//'-Pfafstetter.rst'//' '//&
     'clsm/'//trim(gout)//'-Pfaf.notiny.rst'
-    call system ('cp '//trim(string1))
+    call execute_command_line ('cp '//trim(string1))
     print *,'and, copied those those files to clsm/.'
 
     stop
@@ -2726,9 +2726,9 @@ SUBROUTINE cti_stat_file (ease_grid,gfile, MaskFile)
     
     do n = 1,ip
       if (ease_grid) then  
-         read(10,*,IOSTAT=ierr) typ,pfs,lon,lat,ig,jg,fr_gcm,i_dum
+         read(10,*,IOSTAT=ierr) typ,pfs !,lon,lat,ig,jg,fr_gcm,i_dum
       else
-	 read(10,'(I10,3E20.12,9(2I10,E20.12,I10))',IOSTAT=ierr)             &
+         read(10,'(I10,3E20.12,9(2I10,E20.12,I10))',IOSTAT=ierr)             &
               typ,tarea,lon,lat,ig,jg,fr_gcm,indx_dum,pfs,i_dum,fr_cat,j_dum
       endif
 
@@ -6215,9 +6215,9 @@ SUBROUTINE FUNCIDEP(                                         &
             endif
          enddo
          locmax=MAX(3,indimax10)
-
-      endif      ! if (nmax .ge. shift+1) 
-                  
+         ! add protection here in case nmax <3 . why 3 ?
+         if (locmax > nmax) locmax = nmax
+      endif      ! if (nmax .ge. shift+1)
  30   densmax=denstest(idep,locmax)
       aa(idep)=exp(1.)*densmax
 
@@ -6441,10 +6441,10 @@ SUBROUTINE FUNCZBAR(                                     &
                 densaux(n) .gt. densaux(n-11) .and.        &
                 densaux(n) .gt. densaux(n-12) .and.        &
                 densaux(n) .gt. densaux(n-13) .and.        &
-                densaux(n) .gt. densaux(n-14) .and.        &
-                densaux(n) .gt. densaux(n-15)) then
-               locmax=n
-               goto 30
+                densaux(n) .gt. densaux(n-14)) then ! .and.        &
+                !densaux(n) .gt. densaux(n-15)) then
+                locmax=n
+                goto 30
             endif
          enddo
 
@@ -6459,7 +6459,8 @@ SUBROUTINE FUNCZBAR(                                     &
             endif
          enddo
          locmax=MAX(3,indimax10)
-
+         ! in case nmax < 3. why hard coded 3?
+         if(locmax > nmax) locmax = nmax
       endif      ! if (nmax .ge. shift+1) 
          
  30   densmax=denstest(locmax)
@@ -7250,11 +7251,11 @@ SUBROUTINE REFORMAT_VEGFILES
   integer :: month
 
   tmp_string = 'mkdir -p '//'clsm/g5fmt'
-  call system(tmp_string)
+  call execute_command_line(tmp_string)
   tmp_string = '/bin/mv '//'clsm/lai.dat ' //'clsm/g5fmt/.'
-  call system(tmp_string) 
+  call execute_command_line(tmp_string) 
   tmp_string = '/bin/mv '//'clsm/green.dat ' //'clsm/g5fmt/.'
-  call system(tmp_string) 
+  call execute_command_line(tmp_string) 
 
   open (10,file='clsm/g5fmt/lai.dat'  , form = 'unformatted',   &
        convert='little_endian',status='old',action='read' )
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyTiles.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyTiles.F90
index af4840159..dfbad30e7 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyTiles.F90
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyTiles.F90
@@ -52,7 +52,7 @@ PROGRAM rmTinyTiles
     integer              :: II, JJ, Type
     logical              :: UseType = .false., DoZip=.false., Verb=.false.
     logical              :: Here = .false.
-    integer              :: I, J, status, iargc, nxt
+    integer              :: I, J, status, command_argument_count, nxt
     real*8               :: dx, dy, lon0
     real*8,  allocatable :: xs(:), ys(:)
     character*128        :: &
@@ -69,7 +69,7 @@ PROGRAM rmTinyTiles
 ! Process Arguments
 !------------------
 
-    I = iargc()
+    I = command_argument_count()
 
     if(I < 2 .or. I > 17) then
        print *, "Wrong Number of arguments: ", i
@@ -79,14 +79,14 @@ PROGRAM rmTinyTiles
 
     nxt = 1
 
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
 
     do while(arg(1:1)=='-')
        opt=arg(2:2)
        if(len(trim(arg))==2) then
           if(scan(opt,'zvh')==0) then
              nxt = nxt + 1
-             call getarg(nxt,arg)
+             call get_command_argument(nxt,arg)
           endif
        else
           arg = arg(3:)
@@ -116,18 +116,18 @@ PROGRAM rmTinyTiles
           call exit(1)
        end select
        nxt = nxt + 1
-       call getarg(nxt,arg)
+       call get_command_argument(nxt,arg)
     end do
 
     read(arg,'(i5)') ii
 
     nxt = nxt + 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
 
     read(arg,'(i5)') jj
 
     nxt = nxt + 1
-    call getarg(nxt,arg)
+    call get_command_argument(nxt,arg)
     Grid2 = ARG
  
     if(trim(Gridname) == '') then