Move overland roughness assignment out of time cycling disaggregation…

… routine (#553)

* Add RT disaggregation to read2dlsm()
* Added rt=.logical. to read2dlsm to call Aubrey's regrid_lowres_to_highres while LSM data are all on rank 0. Useful for OV_ROUGH2D and possibly others in the future.

Co-authored-by: Ryan Cabell <[email protected]>
Co-authored-by: Soren Rasmussen <[email protected]>

trunk/NDHMS/Routing/Noah_distr_routing.F

@@ -889,19 +889,17 @@ subroutine disaggregateDomain_drv(did)
                             RT_DOMAIN(did)%NEXP, &
                             rt_domain(did)%overland%properties%distance_to_neighbor, &
                             RT_DOMAIN(did)%INFXSWGT, &
-                            RT_DOMAIN(did)%OVROUGHRTFAC,RT_DOMAIN(did)%LKSATFAC, &
+                            RT_DOMAIN(did)%LKSATFAC, &
                             rt_domain(did)%overland%streams_and_lakes%ch_netrt, &
                             RT_DOMAIN(did)%SH2OWGT, &
                             RT_DOMAIN(did)%subsurface%grid_transform%smcrefrt, &
                             rt_domain(did)%overland%control%infiltration_excess, &
                             RT_DOMAIN(did)%subsurface%grid_transform%smcmaxrt, &
                             RT_DOMAIN(did)%subsurface%grid_transform%smcwltrt, &
                             rt_domain(did)%subsurface%grid_transform%smcrt, &
-                            rt_domain(did)%overland%properties%roughness, &
                             rt_domain(did)%overland%streams_and_lakes%lake_mask, &
                             rt_domain(did)%subsurface%properties%lksatrt, &
                             rt_domain(did)%subsurface%properties%nexprt, &
-                            RT_DOMAIN(did)%OV_ROUGH2d, &
                             RT_DOMAIN(did)%subsurface%properties%sldpth, &
                             RT_DOMAIN(did)%soiltypRT, RT_DOMAIN(did)%soiltyp, &
                             rt_domain(did)%ELRT, RT_DOMAIN(did)%iswater, &
@@ -926,9 +924,9 @@ end subroutine disaggregateDomain_drv
 !===================================================================================================
 subroutine disaggregateDomain(IX, JX, NSOIL, IXRT, JXRT, AGGFACTRT, &
                               SICE, SMC, SH2OX, INFXSRT, area_lsm, SMCMAX1, SMCREF1, &
-                              SMCWLT1, VEGTYP, LKSAT, NEXP, dist,INFXSWGT,OVROUGHRTFAC, &
-                              LKSATFAC, CH_NETRT,SH2OWGT,SMCREFRT, INFXSUBRT,SMCMAXRT, &
-                              SMCWLTRT,SMCRT, OVROUGHRT, LAKE_MSKRT, LKSATRT, NEXPRT, OV_ROUGH2d,  &
+                              SMCWLT1, VEGTYP, LKSAT, NEXP, dist, INFXSWGT, &
+                              LKSATFAC, CH_NETRT, SH2OWGT, SMCREFRT, INFXSUBRT, SMCMAXRT, &
+                              SMCWLTRT, SMCRT, LAKE_MSKRT, LKSATRT, NEXPRT,  &
                               SLDPTH, soiltypRT, soiltyp, elrt, iswater, impervfrac, imperv_adj)
 #ifdef MPP_LAND
    use module_mpp_land, only: left_id,down_id,right_id, &
@@ -957,7 +955,6 @@ subroutine disaggregateDomain(IX, JX, NSOIL, IXRT, JXRT, AGGFACTRT, &
    real, intent(in),  dimension(IX,JX)           :: SMCWLT1 ! coarse grid wilting point
    real, intent(in),  dimension(IX,JX)           :: LKSAT ! coarse grid lateral ksat (m/s)
    real, intent(in),  dimension(IX,JX)           :: NEXP ! coarse grid n exponent
-   real, intent(in), dimension(ix,jx)            :: OV_ROUGH2d ! overland roughness
    ! LSM states:
    real, intent(in),  dimension(IX,JX,NSOIL)     :: SMC ! total soil moisture (m3/m3)
    real, intent(in),  dimension(IX,JX,NSOIL)     :: SH2OX ! liquid soil moisture (m3/m3)
@@ -967,7 +964,6 @@ subroutine disaggregateDomain(IX, JX, NSOIL, IXRT, JXRT, AGGFACTRT, &
    real, intent(in), dimension(IXRT,JXRT,9)      :: dist
                                                     ! routing grid cell distances (m) and area (m2)
                                                     ! TODO: can we just pass in area since we don't need other distances?
-   real, intent(in), dimension(IXRT,JXRT)        :: OVROUGHRTFAC ! overland roughness adj factor
    real, intent(in), dimension(IXRT,JXRT)        :: LKSATFAC ! lateral ksat adj factor
    real, intent(in), dimension(IXRT,JXRT)        :: elrt ! elevation grid (m)
    integer, intent(in), dimension(IXRT,JXRT)     :: CH_NETRT ! channel network routing grid
@@ -987,7 +983,6 @@ subroutine disaggregateDomain(IX, JX, NSOIL, IXRT, JXRT, AGGFACTRT, &
    real, intent(out), dimension(IXRT,JXRT,NSOIL) :: SMCWLTRT ! wilting point on routing grid
    real, intent(out), dimension(IXRT,JXRT)       :: LKSATRT ! lateral ksat on the routing grid (m/s)
    real, intent(out), dimension(IXRT,JXRT)       :: NEXPRT ! n exponent on the routing grid
-   real, intent(out), dimension(IXRT,JXRT)       :: OVROUGHRT ! overland roughness on the routing grid
    ! States:
    real, intent(out), dimension(IX,JX,NSOIL)     :: SICE ! soil ice content on coarse grid (m3/m3)
    real, intent(out), dimension(IXRT,JXRT,NSOIL) :: SMCRT
@@ -1169,25 +1164,6 @@ subroutine disaggregateDomain(IX, JX, NSOIL, IXRT, JXRT, AGGFACTRT, &
 
                ! Now do simple grid remapping tasks
 
-               ! Overland roughness:
-               !DJG map ov roughness as function of land use provided in VEGPARM.TBL...
-               ! --- added extra check for VEGTYP for 'masked-out' locations...
-               ! --- out of basin locations (VEGTYP=0) have OVROUGH hardwired to 0.1
-               IF (VEGTYP(I,J).LE.0) then
-                  OVROUGHRT(IXXRT,JYYRT) = 0.1     !COWS mask test
-               ELSE
-                  OVROUGHRT(IXXRT,JYYRT) = OV_ROUGH2d(i,j)
-                  ! Modify based on impervious fraction
-                  ! See Liong et al. 1989 for linear weighting of "smoothness" (1/roughness)
-                  ! Assuming roughness of 0.02 for impervious and native cell roughness for pervious
-                  if (imperv_adj .ne. 0) then
-                     OVROUGHRT(IXXRT,JYYRT) = 1. / ((1./0.02)*impervfrac(IXXRT,JYYRT) + & ! impervious fraction
-                                                          (1./OVROUGHRT(IXXRT,JYYRT))*(1.-impervfrac(IXXRT,JYYRT))) ! pervious fraction
-                  endif
-                  ! Apply user-supplied adjustment factor
-                  OVROUGHRT(IXXRT,JYYRT) = OVROUGHRT(IXXRT,JYYRT)*OVROUGHRTFAC(IXXRT,JYYRT)
-               END IF
-
                ! Lateral ksat
                !DJG 6.12.08 Map lateral hydraulic conductivity and apply distributed scaling
                ! ---        factor that will be read in from hires terrain file
@@ -1283,7 +1259,6 @@ subroutine disaggregateDomain(IX, JX, NSOIL, IXRT, JXRT, AGGFACTRT, &
    call MPP_LAND_COM_REAL(INFXSUBRT,IXRT,JXRT,99)
    call MPP_LAND_COM_REAL(LKSATRT,IXRT,JXRT,99)
    call MPP_LAND_COM_REAL(NEXPRT,IXRT,JXRT,99)
-   call MPP_LAND_COM_REAL(OVROUGHRT,IXRT,JXRT,99)
    call MPP_LAND_COM_INTEGER(LAKE_MSKRT,IXRT,JXRT,99)
    do i = 1, NSOIL
       call MPP_LAND_COM_REAL(SMCMAXRT(:,:,i),IXRT,JXRT,99)

trunk/NDHMS/Routing/module_HYDRO_io.F

@@ -11171,16 +11171,16 @@ subroutine hdtbl_out(did)
       call hdtbl_out_nc(did,ncid, count,count_flag,"LKSAT",rt_domain(did)%LKSAT,"",ixd,jxd)
       call hdtbl_out_nc(did,ncid, count,count_flag,"NEXP",rt_domain(did)%NEXP,"",ixd,jxd)
    end do
-
 end subroutine hdtbl_out
+
 subroutine hdtbl_in_nc(did)
    implicit none
    integer :: did
    integer :: ierr
    call read2dlsm(did,trim(nlst(did)%hydrotbl_f),"SMCMAX1",rt_domain(did)%SMCMAX1,ierr)
    call read2dlsm(did,trim(nlst(did)%hydrotbl_f),"SMCREF1",rt_domain(did)%SMCREF1,ierr)
    call read2dlsm(did,trim(nlst(did)%hydrotbl_f),"SMCWLT1",rt_domain(did)%SMCWLT1,ierr)
-   call read2dlsm(did,trim(nlst(did)%hydrotbl_f),"OV_ROUGH2D",rt_domain(did)%OV_ROUGH2D,ierr)
+   call read2dlsm(did,trim(nlst(did)%hydrotbl_f),"OV_ROUGH2D",rt_domain(did)%overland%properties%roughness,ierr, rt=.true.)
    call read2dlsm(did,trim(nlst(did)%hydrotbl_f),"LKSAT",rt_domain(did)%LKSAT,ierr)
    call read2dlsm(did,trim(nlst(did)%hydrotbl_f),"NEXP",rt_domain(did)%NEXP,ierr)
    ! Letting this variable be optional and setting to global default value if not found
@@ -11190,37 +11190,99 @@ subroutine hdtbl_in_nc(did)
    endif
 end subroutine hdtbl_in_nc
 
-subroutine read2dlsm(did,file,varName,varOut,ierr)
+subroutine read2dlsm(did,file,varName,varOut,ierr,rt)
 use module_mpp_land,only: mpp_land_bcast_int1
   implicit none
   integer :: did, ncid , iret
   character(len=*) :: file,varName
   real,dimension(:,:) :: varOut
   character(len=256) :: units
   integer, intent(out) :: ierr
-#ifdef MPP_LAND
+  logical, optional, intent(in) :: rt
+  logical :: regrid
+
   real,allocatable,dimension(:,:) :: tmpArr
+
+#ifdef MPP_LAND
   if(my_id .eq. io_id) then
+#endif
      allocate(tmpArr(global_nx,global_ny))
      iret = nf90_open(trim(file), NF90_NOWRITE, ncid)
      call get_2d_netcdf(trim(varName), ncid, tmpArr, units, global_nx, global_ny, &
           .false., ierr)
      iret = nf90_close(ncid)
+#ifdef MPP_LAND
   else
      allocate(tmpArr(1,1))
   endif
-  call decompose_data_real (tmpArr,varOut)
+#endif
+
+  if (present(rt)) then
+    regrid = rt
+  else 
+    regrid = .false.
+  endif
+
+  if (regrid) then
+    call regrid_lowres_to_highres(did, tmpArr, varOut, rt_domain(did)%ixrt, rt_domain(did)%jxrt)
+  else
+    call decompose_data_real (tmpArr,varOut)
+  endif
+
+#ifdef MPP_LAND
   call mpp_land_bcast_int1(ierr)
+#endif
+
   deallocate(tmpArr)
+end subroutine read2dlsm
+
+subroutine regrid_lowres_to_highres(did, lowres_grid, highres_grid, ixrt, jxrt)
+
+  implicit none
+  integer :: did
+  integer :: ixrt, jxrt
+  real, dimension(global_nx, global_ny) :: lowres_grid
+  real, dimension(ixrt,jxrt) :: highres_grid
+  ! Local variables
+  integer :: i, j, irt, jrt, aggfacxrt, aggfacyrt
+
+#ifdef MPP_LAND
+  real,allocatable,dimension(:,:) :: tmpArr
+  if(my_id .eq. io_id) then
+     allocate(tmpArr(global_rt_nx, global_rt_ny))
+#endif
+
+      do j = 1,global_ny ! Start coarse grid j loop
+         do i = 1,global_nx ! Start coarse grid i loop
+
+            do aggfacyrt = nlst(did)%AGGFACTRT-1,0,-1 ! Start disagg fine grid j loop
+            do aggfacxrt = nlst(did)%AGGFACTRT-1,0,-1 ! Start disagg fine grid i loop
+
+               irt = i * nlst(did)%AGGFACTRT - aggfacxrt ! Define fine grid i
+               jrt = j * nlst(did)%AGGFACTRT - aggfacyrt ! Define fine grid j
+#ifdef MPP_LAND
+               ! if(left_id.ge.0) irt = irt + 1
+               ! if(down_id.ge.0) jrt = jrt + 1
+               tmpArr(irt,jrt) = lowres_grid(i,j)
 #else
-     iret = nf90_open(trim(file), NF90_NOWRITE, ncid)
-     call get_2d_netcdf(trim(varName), ncid, varOut, units, rt_domain(did)%ix, rt_domain(did)%jx, &
-          .false., ierr)
-     iret = nf90_close(ncid)
+               highres_grid(irt,jrt) = lowres_grid(i,j)
 #endif
 
-end subroutine read2dlsm
+            end do
+            end do
+
+         end do
+      end do
+
+#ifdef MPP_LAND
+  else
+     allocate(tmpArr(1,1))
+  endif
+  call decompose_RT_real(tmpArr, highres_grid, global_rt_nx, global_rt_ny, ixrt, jxrt)
+  deallocate(tmpArr)
+#endif
 
+end subroutine regrid_lowres_to_highres
 
 subroutine read_channel_only (olddateIn, hgrid, indir, dtbl)
 !use module_HYDRO_io, only: read_rst_crt_reach_nc

trunk/NDHMS/Routing/module_RT.F

@@ -660,6 +660,7 @@ subroutine LandRT_ini(did)
   use config_base, only: nlst, noah_lsm
   use module_RT_data, only: rt_domain
   use module_gw_gw2d_data, only: gw2d
+  use module_HYDRO_io, only: regrid_lowres_to_highres
 #ifdef HYDRO_D
   use module_HYDRO_io, only: output_lake_types
 #endif
@@ -1215,7 +1216,7 @@ subroutine LandRT_ini(did)
 
      !Apply calibration scaling factors to sfc roughness and retention depth here...
      rt_domain(did)%overland%properties%retention_depth = rt_domain(did)%overland%properties%retention_depth * rt_domain(did)%RETDEPRTFAC
-     ! Removing roughness parameter update since it has not been populated yet... currently happens in Noah_dist_routing (AD)
+     rt_domain(did)%overland%properties%roughness = rt_domain(did)%overland%properties%roughness * rt_domain(did)%OVROUGHRTFAC
 
    !ADCHANGE: Moved this channel cell setting from OV_RTNG so it is outside
    !of overland routine (frequently called) and time loop.
@@ -1244,6 +1245,7 @@ subroutine LandRT_ini(did)
 #ifdef MPP_LAND
    ! communicate the value to
      call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%retention_depth,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
+     call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%roughness,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
      call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%surface_slope_x,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
      call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%surface_slope_y,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
      do i = 1, 8