Merge pull request #1743 from agetiran/main

HyMAP-3 updates

Author: jvgeiger

ldt/configs/ldt.config.adoc

@@ -4493,6 +4493,8 @@ This section describes the parameters needed for the HYMAP, HYMAP2 and HYMAP3 ro
 
 `HYMAP node lat:` specifies the name of the grid node latitude data file. This parameter is optional and not used in HYMAP routing models. But it is a "good-to-have" for post-processing purposes.
 
+`HYMAP runoff bias correction map:` specifies the name of the runoff bias correction data file. This is an optional parameter that is required only when "`HYMAP3 runoff bias correction:`" is set to "`1`". If the runoff bias correction option is set to "`0`" or is not specified, the runoff bias correction map is not required.
+
 .Example _ldt.config_ entry
 ....
 HYMAP river width map:                   ../HYMAP_parms/lis_rivwth_0125.bin
@@ -4521,6 +4523,7 @@ HYMAP levee height map:                  ../HYMAP_parms/lis_levee_0125.bin
 Land grid area map:                      ../HYMAP_parms/lis_lgrdare_0125.bin
 HYMAP node lon:                          ../HYMAP_parms/lis_nodelon_0125.bin
 HYMAP node lat:                          ../HYMAP_parms/lis_nodelat_0125.bin
+HYMAP runoff bias correction map:        ../HYMAP_parms/lis_runoff_bc_0125.bin
 ....
 
 This section also outlines the domain specifications of the HYMAP parameter data. For the HYMAP parameters spatial transform option, only "`none`" is supported at this time, and the user is required to input the HYMAP parameters at the grid and resolution of interest.

ldt/params/HYMAP/HYMAP_parmsMod.F90

@@ -53,7 +53,7 @@ module HYMAP_parmsMod
      character*50  :: hymap_proj
      character*50  :: hymap_gridtransform
      integer       :: bfdwicode, rundwicode, rivflocode,urbdroutcode, &
-          leveehgtcode,landareacode,nodeloncode,nodelatcode
+          leveehgtcode,landareacode,nodeloncode,nodelatcode,runoffbiascode
      integer       :: baseflowdelaycode,runoffdelaycode,runoffdelaymcode
      character(LDT_CONST_PATH_LEN) :: riverwidthfile
      character(LDT_CONST_PATH_LEN) :: riverheightfile
@@ -81,6 +81,7 @@ module HYMAP_parmsMod
      character(LDT_CONST_PATH_LEN) :: landgridareafile
      character(LDT_CONST_PATH_LEN) :: nodelonfile
      character(LDT_CONST_PATH_LEN) :: nodelatfile
+     character(LDT_CONST_PATH_LEN) :: runoffbiasfile
      type(LDT_paramEntry) :: hymap_river_width
      type(LDT_paramEntry) :: hymap_river_height
      type(LDT_paramEntry) :: hymap_river_length
@@ -106,6 +107,7 @@ module HYMAP_parmsMod
      type(LDT_paramEntry) :: land_grid_area
      type(LDT_paramEntry) :: hymap_node_lon
      type(LDT_paramEntry) :: hymap_node_lat
+     type(LDT_paramEntry) :: hymap_runoff_bias
   end type hymap_type_dec
 
   type(hymap_type_dec), allocatable :: HYMAP_struc(:)
@@ -144,7 +146,7 @@ subroutine HYMAPparms_init
     integer       :: n
     integer       :: rc
     integer       :: bfdwicode, rundwicode, rivflocode,urbdroutcode, &
-         leveehgtcode,landareacode,nodeloncode,nodelatcode
+         leveehgtcode,landareacode,nodeloncode,nodelatcode,runoffbiascode
     integer       :: baseflowdelaycode,runoffdelaycode,runoffdelaymcode
     integer, allocatable :: nextx(:,:)
     integer, allocatable :: nexty(:,:)
@@ -233,6 +235,8 @@ subroutine HYMAPparms_init
             "HYMAP_node_lon")
        call set_param_attribs(HYMAP_struc(n)%hymap_node_lat,&
             "HYMAP_node_lat")
+       call set_param_attribs(HYMAP_struc(n)%hymap_runoff_bias,&
+            "HYMAP_runoff_bias_correction")
     end do
 
     call ESMF_ConfigFindLabel(LDT_config, &
@@ -348,6 +352,10 @@ subroutine HYMAPparms_init
             LDT_rc%lnc(n),LDT_rc%lnr(n), &
             HYMAP_struc(n)%hymap_node_lat%vlevels))
 
+       allocate(HYMAP_struc(n)%hymap_runoff_bias%value( &
+            LDT_rc%lnc(n),LDT_rc%lnr(n), &
+            HYMAP_struc(n)%hymap_runoff_bias%vlevels))
+
     enddo
 
     call ESMF_ConfigFindLabel(LDT_config,"HYMAP river width map:",rc=rc)
@@ -567,6 +575,16 @@ subroutine HYMAPparms_init
        HYMAP_struc(n)%nodelatcode = nodelatcode
     enddo
 
+    call ESMF_ConfigFindLabel(LDT_config,"HYMAP runoff bias correction map:",rc=runoffbiascode)
+    do n=1,LDT_rc%nnest
+       call ESMF_ConfigGetAttribute(LDT_config, &
+            HYMAP_struc(n)%runoffbiasfile,rc=rc)
+       if (runoffbiascode.eq.0) then
+          call LDT_verify(rc,'HYMAP runoff bias correction map: not specified')
+       endif
+       HYMAP_struc(n)%runoffbiascode = runoffbiascode
+    enddo
+
     write(LDT_logunit,*) &
          " - - - - - - - - HYMAP Router Parameters - - - - - - - - - -"
 
@@ -889,6 +907,18 @@ subroutine HYMAPparms_init
                'Done reading '//trim(HYMAP_struc(n)%nodelatfile)
        endif
 
+       if (HYMAP_struc(n)%runoffbiascode.eq.0) then
+          call LDT_gridOptChecks( n, "HYMAP runoff bias correction", &
+               HYMAP_struc(n)%hymap_gridtransform, hymap_proj, &
+               hymapparms_gridDesc(n,9) )
+          write(LDT_logunit,*) &
+               'reading '//trim(HYMAP_struc(n)%runoffbiasfile)
+          call read_HYMAP_runoff_bias( &
+               n,HYMAP_struc(n)%hymap_runoff_bias%value(:,:,1))
+          write(LDT_logunit,*) &
+               'Done reading '//trim(HYMAP_struc(n)%runoffbiasfile)
+       endif
+
     enddo
 
   end subroutine HYMAPparms_init
@@ -995,6 +1025,11 @@ subroutine HYMAPparms_writeHeader(n,ftn,dimID,monthID)
        call LDT_writeNETCDFdataHeader(n,ftn,dimID, &
             HYMAP_struc(n)%land_grid_area)
     endif
+               
+    if (HYMAP_struc(n)%runoffbiascode.eq.0) then
+       call LDT_writeNETCDFdataHeader(n,ftn,dimID, &
+            HYMAP_struc(n)%hymap_runoff_bias)
+    endif      
 
   end subroutine HYMAPparms_writeHeader
 
@@ -1065,6 +1100,9 @@ subroutine HYMAPparms_writeData(n,ftn)
     if (HYMAP_struc(n)%nodelatcode.eq.0) then
       call LDT_writeNETCDFdata(n,ftn,HYMAP_struc(n)%hymap_node_lat)
     endif
+    if (HYMAP_struc(n)%runoffbiascode.eq.0) then
+      call LDT_writeNETCDFdata(n,ftn,HYMAP_struc(n)%hymap_runoff_bias)
+    endif
 
   end subroutine HYMAPparms_writeData
 

ldt/params/HYMAP/read_HYMAP_runoff_bias.F90

@@ -0,0 +1,78 @@
+!-----------------------BEGIN NOTICE -- DO NOT EDIT-----------------------
+! NASA Goddard Space Flight Center
+! Land Information System Framework (LISF)
+! Version 7.5
+!
+! Copyright (c) 2024 United States Government as represented by the
+! Administrator of the National Aeronautics and Space Administration.
+! All Rights Reserved.
+!-------------------------END NOTICE -- DO NOT EDIT-----------------------
+!BOP
+!
+! !ROUTINE: read_HYMAP_runoff_bias
+! \label{read_HYMAP_runoff_bias}
+!
+! !REVISION HISTORY:
+!  13 Nov 2025: Augusto Getirana; Initial Specification
+!
+! !INTERFACE:
+subroutine read_HYMAP_runoff_bias(n, array)
+
+! !USES:
+  use ESMF
+  use LDT_coreMod,       only : LDT_rc, LDT_domain
+  use LDT_logMod,        only : LDT_logunit, LDT_getNextUnitNumber, &
+       LDT_releaseUnitNumber, LDT_endrun
+  use LDT_fileIOMod,     only : readLISdata 
+  use HYMAP_parmsMod
+
+  implicit none
+! !ARGUMENTS: 
+  integer,      intent(in) :: n
+  real,      intent(inout) :: array(LDT_rc%lnc(n),LDT_rc%lnr(n),1)
+
+! !DESCRIPTION:
+!  This subroutine retrieves runoff bias correction maps
+!  
+!  The arguments are:
+!  \begin{description}
+!  \item[n]
+!   index of the nest
+!  \item[mo]
+!   time index (month or quarter)
+!  \item[array]
+!   output field with the retrieved greenness fraction
+!  \end{description}
+!
+!EOP      
+  integer :: ftn
+  integer :: c,r
+  logical :: file_exists
+
+  ftn = LDT_getNextUnitNumber()
+
+  inquire(file=trim(HYMAP_struc(n)%runoffbiasfile), exist=file_exists)
+  if(.not.file_exists) then 
+     write(LDT_logunit,*) '[ERR] HYMAP runoff bias correction map, ',&
+           trim(HYMAP_struc(n)%runoffbiasfile),', not found.'
+     write(LDT_logunit,*) 'Program stopping ...'
+     call LDT_endrun
+  endif
+
+  open(ftn, file=trim(HYMAP_struc(n)%runoffbiasfile), access='direct',&
+       status='old', form="unformatted", convert="big_endian", recl=4)
+  
+  call readLISdata(n, ftn, HYMAP_struc(n)%hymap_proj, &
+       HYMAP_struc(n)%hymap_gridtransform, &
+       HYMAP_struc(n)%hymapparms_gridDesc(:), 1, array)  ! 1 indicates 2D layer
+  
+  do r=1,LDT_rc%lnr(n)
+     do c=1,LDT_rc%lnc(n)
+        if(array(c,r,1).lt.0) then
+           array(c,r,1) = LDT_rc%udef
+        endif
+     enddo
+  enddo
+  call LDT_releaseUnitNumber(ftn)
+
+end subroutine read_HYMAP_runoff_bias

lis/configs/lis.config.adoc

@@ -11163,6 +11163,17 @@ active restart file.
 `HYMAP3 2-way coupling flooded fraction threshold:` specifies the flooded fraction
 threshold used in the 2-way coupling mode.
 
+`HYMAP3 runoff bias correction:` specifies whether to account for a runoff bias correction. This option is elective and, if not defined, HYMAP3 does not account for runoff bias correction. This option requires a runoff bias correction field provided in the LIS domain and parameter data file (HYMAP_runoff_bias_correction).
+
+Acceptable values are:
+
+|====
+|Value | Description
+
+|0     | Do not use
+|1     | Use 
+|====
+
 .Example _lis.config_ entry
 ....
 HYMAP3 routing model time step: "30mn"
@@ -11185,6 +11196,7 @@ HYMAP3 routing model restart interval: "1mo"
 HYMAP3 routing model restart file: "./OL/ROUTING/197901/LIS_RST_HYMAP3_router_197901030000.d01.nc
 HYMAP3 enable 2-way coupling: 0
 HYMAP3 2-way coupling flooded fraction threshold: 0.1
+HYMAP3 runoff bias correction: 1
 ....
 
 

lis/routing/HYMAP2_router/runoffdata/LISoutput/LISrunoffdataMod.F90

@@ -49,7 +49,7 @@ module LISrunoffdataMod
      integer             :: nc,nr
      character(len=LIS_CONST_PATH_LEN)       :: previous_filename
      real                :: datares
-     real, allocatable   :: qs(:,:),qsb(:,:),evap(:,:)
+     real, allocatable   :: qs(:,:),qsb(:,:)
 
   end type LISrunoffdatadec
 
@@ -214,10 +214,8 @@ subroutine LISrunoffdata_init
       LISrunoffdata_struc(n)%previous_filename='none'
       allocate(LISrunoffdata_struc(n)%qs(LIS_rc%lnc(n),LIS_rc%lnr(n)))
       allocate(LISrunoffdata_struc(n)%qsb(LIS_rc%lnc(n),LIS_rc%lnr(n)))
-      allocate(LISrunoffdata_struc(n)%evap(LIS_rc%lnc(n),LIS_rc%lnr(n)))
       LISrunoffdata_struc(n)%qs=LIS_rc%udef
       LISrunoffdata_struc(n)%qsb=LIS_rc%udef
-      LISrunoffdata_struc(n)%evap=LIS_rc%udef
     enddo
 
   end subroutine LISrunoffdata_init

lis/routing/HYMAP2_router/runoffdata/LISoutput/readLISrunoffdata.F90

@@ -33,7 +33,6 @@ subroutine readLISrunoffdata(n,surface_runoff, baseflow)
   integer,          intent(in) :: n
   real                         :: surface_runoff(LIS_rc%lnc(n),LIS_rc%lnr(n))
   real                         :: baseflow(LIS_rc%lnc(n),LIS_rc%lnr(n))
-  real                         :: total_evapotranspiration(LIS_rc%lnc(n),LIS_rc%lnr(n))
 
   !caveats 
   ! 1) assumes the LIS outputs are the same output interval as that of
@@ -46,20 +45,19 @@ subroutine readLISrunoffdata(n,surface_runoff, baseflow)
   ! 
   ! 5) LIS outputs are in NETCDF format. 
   !
-  !Added total evapotranspiration (Evap)
   integer                       :: c,r,t
   integer, allocatable      :: nqs(:,:)
 
   !ag - 17Mar2016
   !real,   allocatable       :: qs(:,:),qs_t(:)
   real                  :: qs2d(LISrunoffdata_struc(n)%nc,LISrunoffdata_struc(n)%nr)
   real                  :: qsb2d(LISrunoffdata_struc(n)%nc,LISrunoffdata_struc(n)%nr)
-  real,   allocatable   :: qs_t(:),qsb_t(:),evap_t(:)
+  real,   allocatable   :: qs_t(:),qsb_t(:)
   logical*1             :: lb(LISrunoffdata_struc(n)%nc*LISrunoffdata_struc(n)%nr)
   real                  :: var_input(LISrunoffdata_struc(n)%nc*LISrunoffdata_struc(n)%nr)
   logical*1             :: lo(LIS_rc%lnc(n)*LIS_rc%lnr(n))
   real                  :: var_out(LIS_rc%lnc(n)*LIS_rc%lnr(n))
-  integer               :: ios, nid,qsid,qsbid,evapid
+  integer               :: ios, nid,qsid,qsbid
   character(len=LIS_CONST_PATH_LEN) :: filename
   logical               :: file_exists
   logical               :: check_Flag
@@ -94,9 +92,6 @@ subroutine readLISrunoffdata(n,surface_runoff, baseflow)
       ios = nf90_inq_varid(nid,'Qsb_tavg',qsbid)
       call LIS_verify(ios,'failed to read Qsb_tavg field in readLISrunoffdata')
 
-      ios = nf90_inq_varid(nid,'Evap_tavg',evapid)
-      call LIS_verify(ios,'failed to read Evap_tavg field in readLISrunoffdata')
-
       if(LISrunoffdata_struc(n)%domainCheck) then 
          if(LIS_rc%wopt.eq."2d gridspace") then 
             ios = nf90_get_var(nid,qsid,LISrunoffdata_struc(n)%qs, &
@@ -109,10 +104,6 @@ subroutine readLISrunoffdata(n,surface_runoff, baseflow)
                  count=(/LIS_rc%lnc(n),LIS_rc%lnr(n)/))
             call LIS_verify(ios, 'failed to read Qsb_tavg field in readLISrunoffdata')
 
-            ios = nf90_get_var(nid,evapid,LISrunoffdata_struc(n)%evap,&
-                 start=(/LIS_ews_halo_ind(n,LIS_localPet+1),LIS_nss_halo_ind(n,LIS_localPet+1)/),&
-                 count=(/LIS_rc%lnc(n),LIS_rc%lnr(n)/))
-            call LIS_verify(ios, 'failed to read Evap_tavg field in readLISrunoffdata')
          else
             write(LIS_logunit,*) "Stand-alone HYMAP is only supported for '2d gridspace' outputs currently"
             call LIS_endrun()
@@ -247,11 +238,9 @@ subroutine readLISrunoffdata(n,surface_runoff, baseflow)
   where(LISrunoffdata_struc(n)%qs/=LIS_rc%udef)
     surface_runoff = LISrunoffdata_struc(n)%qs
     baseflow = LISrunoffdata_struc(n)%qsb
-    total_evapotranspiration = LISrunoffdata_struc(n)%evap
   else where
     surface_runoff = 0.0
     baseflow = 0.0
-    total_evapotranspiration = 0.0
   end where
 
 end subroutine readLISrunoffdata