Merge pull request #447 from karsenau/fix/410_ldtbuffer

Fix update for domain-added buffer in LDT

ldt/configs/ldt.config.adoc

@@ -389,6 +389,24 @@ LDT output directory:       OUTPUT
 Undefined value:             -9999.0
 ....
 
+`Add buffer to parameter grid domain:` adds a set buffer around a parameter file target domain. 
+Acceptable values are:
+
+[cols="<,<",]
+|===
+|Value |Description
+
+|"0" |No buffer added
+|"1" |Buffer included
+|===
+
+The default value is 1.
+
+.Example _ldt.config_ entry
+....
+Add buffer to parameter grid domain:     0
+....
+
 `Number of ensembles per tile:` specifies the number of ensembles of tiles to be used. The value should be greater than or equal to 1.
 
 .Example _ldt.config_ entry

ldt/core/LDT_PRIV_rcMod.F90

@@ -34,6 +34,7 @@ module LDT_PRIV_rcMod
      character*100, allocatable :: paramAttribsFile(:)
      character*50           :: lis_map_proj
      real                   :: lis_map_resfactor
+     integer                :: add_buffer
 
 ! -- Land surface input parameters:
      integer                :: max_model_types 

ldt/core/LDT_gridmappingMod.F90

@@ -77,6 +77,9 @@ subroutine LDT_RunDomainPts( n, param_proj, param_grid, &
    real,allocatable :: lats(:)
    real             :: diff_lon
    integer          :: nlats
+   integer          :: buffer_flag
+   real             :: rlat1, rlat2, rlon1, rlon2
+   real             :: plat1, plat2, plon1, plon2
 ! _____________________________________________________
 
    glpnr = 0
@@ -131,17 +134,49 @@ subroutine LDT_RunDomainPts( n, param_proj, param_grid, &
    lisdom_yres_ur = abs(rlat(LDT_rc%lnc(n),  LDT_rc%lnr(n))  &
                   - rlat(LDT_rc%lnc(n),  LDT_rc%lnr(n)-1) )
 
+   buffer_flag = 1
+   ! Checks for when to implement a "buffer" for specific domains and
+   !  projections ...
+   rlat1=rlat(1,1);  rlat2=rlat(LDT_rc%lnc(n),LDT_rc%lnr(n))
+   rlon1=rlon(1,1);  rlon2=rlon(LDT_rc%lnc(n),LDT_rc%lnr(n))
+   plat1=param_grid(4);  plat2=param_grid(7) 
+   plon1=param_grid(5);  plon2=param_grid(8)
+
+  !-- Checks for turning off buffer domain ...
+   ! Parameter geographic extents match LIS run domain extents:
+   if( rlat1==plat1 .and. rlat2==plat2 .and. &
+       rlon1==plon1 .and. rlon2==plon2 ) then
+     buffer_flag = 0   ! Buffer not needed
+   endif
+   ! Turn off buffer when latlon x/y resolutions are the same:
+   if( param_proj == "latlon" .and. &
+       param_grid(9) == lisdom_xres_ll.and. &
+       param_grid(10) == lisdom_yres_ll )then
+     buffer_flag = 0   ! Buffer not needed
+   endif 
+   if( LDT_rc%add_buffer == 0 ) then   ! New LDT config option
+     buffer_flag = 0
+   endif
+
    ! Incorporate buffer for target grid (e.g., helps with curvilinear projections)
-   lisdom_min_lat = max((rlat(1,1)-(5*lisdom_yres_ll)),-90.0)
-   lisdom_max_lat = min((rlat(LDT_rc%lnc(n),LDT_rc%lnr(n))+(5*lisdom_yres_ll)),90.0)
-
-   ! Account for crossing IDL:
-   if( rlon(1,1) <= rlon(LDT_rc%lnc(n),LDT_rc%lnr(n)) ) then
-     lisdom_min_lon = max((rlon(1,1)-(5*lisdom_xres_ll)),-180.0)
-     lisdom_max_lon = min((rlon(LDT_rc%lnc(n),LDT_rc%lnr(n))+(5*lisdom_xres_ll)),180.0)
-   else  
-     lisdom_min_lon = max((rlon(1,1)-(5*lisdom_xres_ll)),0.0)
-     lisdom_max_lon = min((rlon(LDT_rc%lnc(n),LDT_rc%lnr(n))+(5*lisdom_xres_ll)),0.0)
+   if( buffer_flag == 1 ) then 
+     write(LDT_logunit,*) "[INFO] Incorporating buffer around subsetted grid domain ... "
+     lisdom_min_lat = max((rlat(1,1)-(5*lisdom_yres_ll)),-90.0)
+     lisdom_max_lat = min((rlat(LDT_rc%lnc(n),LDT_rc%lnr(n))+(5*lisdom_yres_ll)),90.0)
+     ! Account for crossing IDL:
+     if( rlon(1,1) <= rlon(LDT_rc%lnc(n),LDT_rc%lnr(n)) ) then
+       lisdom_min_lon = max((rlon(1,1)-(5*lisdom_xres_ll)),-180.0)
+       lisdom_max_lon = min((rlon(LDT_rc%lnc(n),LDT_rc%lnr(n))+(5*lisdom_xres_ll)),180.0)
+     else  
+       lisdom_min_lon = max((rlon(1,1)-(5*lisdom_xres_ll)),0.0)
+       lisdom_max_lon = min((rlon(LDT_rc%lnc(n),LDT_rc%lnr(n))+(5*lisdom_xres_ll)),0.0)
+     endif
+   else  ! No buffer applied
+!     write(LDT_logunit,*) "[INFO] No buffer incorporated ... "
+     lisdom_min_lat = rlat(1,1)
+     lisdom_max_lat = rlat(LDT_rc%lnc(n),LDT_rc%lnr(n))
+     lisdom_min_lon = rlon(1,1)
+     lisdom_max_lon = rlon(LDT_rc%lnc(n),LDT_rc%lnr(n))
    endif
 
 ! -------------------------------------------------------------
@@ -157,7 +192,6 @@ subroutine LDT_RunDomainPts( n, param_proj, param_grid, &
 
   !- Calculate total points for global (or "full") parameter file domains:
      glpnr = nint((param_grid(7)-param_grid(4))/param_grid(10)) + 1
-!     glpnc = nint((param_grid(8)-param_grid(5))/param_grid(9) ) + 1
      glpnc = nint(diff_lon(param_grid(8),param_grid(5))/param_grid(9) ) + 1
 
   !- LIS RUN DOMAIN GRID INFORMATION: 
@@ -331,7 +365,6 @@ subroutine LDT_RunDomainPts( n, param_proj, param_grid, &
    !! Assemble the parameter subdomain extents and other info:
 
      ! Estimate final number of subsetted parameter points:
-!     subparam_nc = nint((subparm_urlon_ext - subparm_lllon_ext)/param_grid(9) ) + 1
      subparam_nc = nint(diff_lon(subparm_urlon_ext,subparm_lllon_ext)/param_grid(9) ) + 1
      subparam_nr = nint((subparm_urlat_ext - subparm_lllat_ext)/param_grid(10)) + 1
 
@@ -408,21 +441,21 @@ subroutine LDT_RunDomainPts( n, param_proj, param_grid, &
 
    case ( "gaussian" )
 
-      !glpnc = nint((param_grid(8)-param_grid(5))/param_grid(9)) + 1
-      glpnc = nint(diff_lon(param_grid(8),param_grid(5))/param_grid(9)) + 1
-      glpnr = param_grid(3)
+     !glpnc = nint((param_grid(8)-param_grid(5))/param_grid(9)) + 1
+     glpnc = nint(diff_lon(param_grid(8),param_grid(5))/param_grid(9)) + 1
+     glpnr = param_grid(3)
 
-    ! Assume that parameter grid resolution SAME as LIS run grid resolution
-      subparm_lllat_ext = lisdom_min_lat
-      subparm_lllon_ext = lisdom_min_lon
-      subparm_urlat_ext = lisdom_max_lat
-      subparm_urlon_ext = lisdom_max_lon
+     ! Assume that parameter grid resolution SAME as LIS run grid resolution
+     subparm_lllat_ext = lisdom_min_lat
+     subparm_lllon_ext = lisdom_min_lon
+     subparm_urlat_ext = lisdom_max_lat
+     subparm_urlon_ext = lisdom_max_lon
 
-  !- Estimate final number of subsetted parameter points:
-      subparam_nc = LDT_rc%lnc(n)
-      subparam_nr = LDT_rc%lnr(n) 
+     ! Estimate final number of subsetted parameter points:
+     subparam_nc = LDT_rc%lnc(n)
+     subparam_nr = LDT_rc%lnr(n) 
 
-  !- Set subsetted parameter grid array inputs:
+     ! Set subsetted parameter grid array inputs:
      subparam_gridDesc(1)  = 4  
      subparam_gridDesc(2)  = float(subparam_nc)
      subparam_gridDesc(3)  = float(subparam_nr)

ldt/core/LDT_readConfig.F90

@@ -382,6 +382,11 @@ subroutine LDT_readConfig(configfile)
        rc=rc)
   call LDT_verify(rc,'Undefined value: not defined')
 
+! Option to add buffer around parameter grid domain:
+  call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%add_buffer,&
+       label="Add buffer to parameter grid domain:",&
+       default=1,rc=rc)
+  call LDT_verify(rc,'Add buffer to parameter grid domain: not defined')
 
 ! Set number of processors along x and y directions:
   call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%npesx,&

ldt/domains/latlon/readinput_latlon.F90

@@ -130,11 +130,10 @@ subroutine readinput_latlon
           LDT_rc%gridDesc(n,4),LDT_rc%gridDesc(n,7),&
           LDT_rc%gridDesc(n,5),LDT_rc%gridDesc(n,8)
 
-     LDT_rc%gridDesc(n,1) = 0
-     LDT_rc%gridDesc(n,9) = run_dd(n,5)
-
+     LDT_rc%gridDesc(n,1) = 0                ! latlon grid
+     LDT_rc%gridDesc(n,9) = run_dd(n,5)      ! dx
      if(LDT_rc%gridDesc(n,1).eq.0) then 
-        LDT_rc%gridDesc(n,10) = run_dd(n,6)
+        LDT_rc%gridDesc(n,10) = run_dd(n,6)  ! dy
         LDT_rc%gridDesc(n,6)  = 128
         LDT_rc%gridDesc(n,11) = 64
         LDT_rc%gridDesc(n,20) = 64
@@ -152,11 +151,11 @@ subroutine readinput_latlon
                              LDT_rc%gridDesc(n,8),LDT_rc%gridDesc(n,5)
      endif
 
-     ! Difference in number of longitudes (dx):
+     ! Difference in number of longitudes (nx):
      LDT_rc%gridDesc(n,2) = nint(diff_lon(LDT_rc%gridDesc(n,8),LDT_rc%gridDesc(n,5))&
                           / LDT_rc%gridDesc(n,9)) + 1   
 
-     ! Difference in number of latitudes (dy):
+     ! Difference in number of latitudes (ny):
      LDT_rc%gridDesc(n,3) = nint((LDT_rc%gridDesc(n,7)-LDT_rc%gridDesc(n,4))&
                           / LDT_rc%gridDesc(n,10)) + 1    
 

ldt/params/HYMAP/HYMAP_parmsMod.F90

@@ -15,7 +15,7 @@ module HYMAP_parmsMod
 !
 !  29 Mar 2013: Sujay Kumar; Initial implementation
 !   2 Dec 2015: Augusto Getirana: Included drainage area and basin maps
-!   1 Nov 2017: Augusto Getirana: Included flow type maps and, baseflow and and surface runoff dwi maps
+!   1 Nov 2017: Augusto Getirana: Included flow type maps, baseflow and surface runoff dwi maps
 !
   use ESMF
   use LDT_coreMod

ldt/params/HYMAP/read_HYMAP_baseflow_delay.F90

@@ -48,7 +48,8 @@ subroutine read_HYMAP_baseflow_delay(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%baseflowdelayfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Baseflowdelay map ',trim(HYMAP_struc(n)%baseflowdelayfile),' not found'
+     write(LDT_logunit,*) '[ERR] Baseflowdelay map, ',&
+           trim(HYMAP_struc(n)%baseflowdelayfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_baseflow_dwi_ratio.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_baseflow_dwi_ratio(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%baseflowdwiratiofile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'baseflow DWI ratio map ',trim(HYMAP_struc(n)%baseflowdwiratiofile),' not found'
+     write(LDT_logunit,*) '[ERR] baseflow DWI ratio map, ',&
+          trim(HYMAP_struc(n)%baseflowdwiratiofile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_basin.F90

@@ -47,7 +47,8 @@ subroutine read_HYMAP_basin(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%basinfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Basin map ',trim(HYMAP_struc(n)%basinfile),' not found'
+     write(LDT_logunit,*) '[ERR] Basin map, ',&
+           trim(HYMAP_struc(n)%basinfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_basin_mask.F90

@@ -47,7 +47,8 @@ subroutine read_HYMAP_basin_mask(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%basinmaskfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Basin mask map ',trim(HYMAP_struc(n)%basinmaskfile),' not found'
+     write(LDT_logunit,*) '[ERR] Basin mask map, ',&
+           trim(HYMAP_struc(n)%basinmaskfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_drain_area.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_drain_area(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%drainareafile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Drainarea map ',trim(HYMAP_struc(n)%drainareafile),' not found'
+     write(LDT_logunit,*) '[ERR] Drainage area map, ',&
+           trim(HYMAP_struc(n)%drainareafile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_fld_height.F90

@@ -55,7 +55,8 @@ subroutine read_HYMAP_fld_height(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%fldheightfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Fldheight map ',trim(HYMAP_struc(n)%fldheightfile),' not found'
+     write(LDT_logunit,*) '[ERR] Fldheight map, ',trim(HYMAP_struc(n)%fldheightfile),&
+                          ', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_fld_z.F90

@@ -48,7 +48,8 @@ subroutine read_HYMAP_fld_z(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%fldzfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Fldz map ',trim(HYMAP_struc(n)%fldzfile),' not found'
+     write(LDT_logunit,*) '[ERR] Fldz map, ',&
+           trim(HYMAP_struc(n)%fldzfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_flow_dir_x.F90

@@ -48,7 +48,8 @@ subroutine read_HYMAP_flow_dir_x(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%flowdirxfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Flowdirx map ',trim(HYMAP_struc(n)%flowdirxfile),' not found'
+     write(LDT_logunit,*) '[ERR] Flowdirx map, ',&
+           trim(HYMAP_struc(n)%flowdirxfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_flow_dir_y.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_flow_dir_y(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%flowdiryfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Flowdiry map ',trim(HYMAP_struc(n)%flowdiryfile),' not found'
+     write(LDT_logunit,*) '[ERR] Flowdiry map, ',&
+          trim(HYMAP_struc(n)%flowdiryfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_flow_type.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_flow_type(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%flowtypefile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Flow type map ',trim(HYMAP_struc(n)%flowtypefile),' not found'
+     write(LDT_logunit,*) '[ERR] Flow type map, ',&
+           trim(HYMAP_struc(n)%flowtypefile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_grid_area.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_grid_area(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%gridareafile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Gridarea map ',trim(HYMAP_struc(n)%gridareafile),' not found'
+     write(LDT_logunit,*) '[ERR] Gridarea map, ',&
+           trim(HYMAP_struc(n)%gridareafile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_grid_dist.F90

@@ -48,7 +48,8 @@ subroutine read_HYMAP_grid_dist(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%griddistfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Griddist map ',trim(HYMAP_struc(n)%griddistfile),' not found'
+     write(LDT_logunit,*) '[ERR] Griddist map, ',&
+           trim(HYMAP_struc(n)%griddistfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_grid_elev.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_grid_elev(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%gridelevfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Gridelev map ',trim(HYMAP_struc(n)%gridelevfile),' not found'
+     write(LDT_logunit,*) '[ERR] Gridelev map, ',&
+           trim(HYMAP_struc(n)%gridelevfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_refq.F90

@@ -47,7 +47,8 @@ subroutine read_HYMAP_refq(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%refqfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Refq map ',trim(HYMAP_struc(n)%refqfile),' not found'
+     write(LDT_logunit,*) '[ERR] Refq map, ',&
+           trim(HYMAP_struc(n)%refqfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_river_height.F90

@@ -48,7 +48,8 @@ subroutine read_HYMAP_river_height(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%riverheightfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Riverheight map ',trim(HYMAP_struc(n)%riverheightfile),' not found'
+     write(LDT_logunit,*) '[ERR] River height map, ',&
+           trim(HYMAP_struc(n)%riverheightfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_river_length.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_river_length(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%riverlengthfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Riverlength map ',trim(HYMAP_struc(n)%riverlengthfile),' not found'
+     write(LDT_logunit,*) '[ERR] River length map, ',&
+           trim(HYMAP_struc(n)%riverlengthfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_river_width.F90

@@ -47,7 +47,8 @@ subroutine read_HYMAP_river_width(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%riverwidthfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Riverwidth map ',trim(HYMAP_struc(n)%riverwidthfile),' not found'
+     write(LDT_logunit,*) '[ERR] River width map, ',&
+           trim(HYMAP_struc(n)%riverwidthfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_river_z.F90

@@ -49,7 +49,8 @@ subroutine read_HYMAP_river_z(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%riverzfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Riverz map ',trim(HYMAP_struc(n)%riverzfile),' not found'
+     write(LDT_logunit,*) '[ERR] Riverz map, ',&
+           trim(HYMAP_struc(n)%riverzfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif

ldt/params/HYMAP/read_HYMAP_runoff_delay.F90

@@ -48,7 +48,8 @@ subroutine read_HYMAP_runoff_delay(n, array)
 
   inquire(file=trim(HYMAP_struc(n)%runoffdelayfile), exist=file_exists)
   if(.not.file_exists) then 
-     write(LDT_logunit,*) 'Runoffdelay map ',trim(HYMAP_struc(n)%runoffdelayfile),' not found'
+     write(LDT_logunit,*) '[ERR] Runoff delay map, ',&
+           trim(HYMAP_struc(n)%runoffdelayfile),', not found.'
      write(LDT_logunit,*) 'Program stopping ...'
      call LDT_endrun
   endif