make_bcs fix for bad HWSDv1.21 data in Argentina "peatland"

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc

@@ -120,7 +120,7 @@
 #           - backfilled with global land average snow albedo where unavailable
 #           - must use compatible bcs version that includes MODIS-based snow albedo (e.g., v06, v08, v09, ...)
 #           - NOTE: bcs v06, v08, and v09 used approximate averaging of MODIS-based snow albedo to tile space;
-#             later bcs versions employ more accurate, raster-based averaging       
+#                   bcs v11 and v12 employ more accurate, raster-based averaging.
 #
 # GEOSagcm=>SNOW_ALBEDO_INFO: 0
 # GEOSldas=>SNOW_ALBEDO_INFO: 0

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt

@@ -51,4 +51,9 @@ ecbuild_add_executable (TARGET mkEASETilesParam.x SOURCES mkEASETilesParam.F90 L
 
 install(PROGRAMS clsm_plots.pro create_README.csh  DESTINATION bin)
 file(GLOB MAKE_BCS_PYTHON CONFIGURE_DEPENDS "./make_bcs*.py")
+list(FILTER MAKE_BCS_PYTHON EXCLUDE REGEX "make_bcs_shared.py")
 install(PROGRAMS  ${MAKE_BCS_PYTHON} DESTINATION bin)
+
+set(file ./make_bcs_shared.py)
+configure_file(${file} ${file} @ONLY)
+install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/${file} DESTINATION bin)

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/create_README.csh

@@ -620,7 +620,7 @@ APPENDIX I - mkCatchParam input options and log ................................
 	 file name: soil_param.dat
 	 do n = 1, ${NTILES}
 _EOI_
-if( $mysoil == HWSD ) then
+if( $mysoil == HWSD || $mysoil == HWSD_b ) then
 cat << _EOS1_ > clsm/soil
 		read ([UNIT],'(i10,i8,i4,i4,3f8.4,f12.8,f7.4,f10.4,3f7.3,4f7.3,2f10.4, f8.4)')       &
                      tile_index, pfaf_code, soil_class_top, soil_class_com, BEE,                     &

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py

@@ -194,6 +194,7 @@ def ask_questions(default_grid="Cubed-Sphere"):
    "v09 : NL3 + PEATMAP + MODIS snow alb", \
    "v10 : NL3 + PEATMAP + MODIS snow alb v2", \
    "v11 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2", \
+   "v12 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix", \
    "ICA : Icarus        (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus/)", \
    "GM4 : Ganymed-4_0   (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Ganymed-4_0/)", \
    "F25 : Fortuna-2_5   (archived*: n/a)"], 

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py

@@ -5,18 +5,23 @@
 import os
 import glob
 
+BUILT_ON_SLES15 = "@BUILT_ON_SLES15@"
 def get_script_head() :
 
-   return  """#!/bin/csh -x
+  head =  """#!/bin/csh -x
 
 #SBATCH --output={EXPDIR}/{TMP_DIR}/logs/{GRIDNAME}/{GRIDNAME2}.log
 #SBATCH --error={EXPDIR}/{TMP_DIR}/logs/{GRIDNAME}/{GRIDNAME2}.err
 #SBATCH --account={account}
 #SBATCH --time=12:00:00
 #SBATCH --nodes=1
 #SBATCH --job-name={GRIDNAME2}.j
-#SBATCH --constraint=sky|cas
+"""
+  constraint = "#SBATCH --constraint=sky|cas"
+  if BUILT_ON_SLES15 :
+     constraint = "#SBATCH --constraint=mil"
 
+  head = head + constraint + """
 echo "-----------------------------" 
 echo "make_bcs starts date/time" 
 echo `date` 
@@ -37,6 +42,7 @@ def get_script_head() :
   mkdir -p geometry land/shared til rst data/MOM5 data/MOM6 clsm/plots
 endif
 """
+  return head
 
 def get_change_til_file(grid_type):
   script = ""
@@ -116,6 +122,11 @@ def get_script_mv(grid_type):
 /bin/mv clsm/lnfm.dat    land/{GRIDNAME}/lnfm_clim_{RC}.data
 /bin/mv clsm/ndvi.dat    land/{GRIDNAME}/ndvi_clim_{RC}.data
 
+# vegdyn_{RC}.dat file is nc4; for clarification, create link with proper file name extension
+cd land/{GRIDNAME}
+ln -s vegdyn_{RC}.dat vegdyn_{RC}.nc4
+cd ../../
+
 /bin/mv clsm/ar.new \\
         clsm/bf.dat \\
         clsm/ts.dat \\

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90

@@ -545,7 +545,7 @@ PROGRAM mkCatchParam
        ! this block is for n_threads>1
        !==============================
 
-       if(SOILBCS=='NGDC') then
+       if(trim(SOILBCS)=='NGDC') then
           write (log_file,'(a)')'Creating (intermediate) NGDC soil types file...'
           call create_soil_types_files (nc,nr,ease_grid,fnameTil,fnameRst)    
           write (log_file,'(a)')'         Done.'
@@ -562,11 +562,11 @@ PROGRAM mkCatchParam
        inquire(file=trim(fname_tmp), exist=file_exists)
        if (.not.file_exists) then
           write (log_file,'(a)')'         Creating file...'
-          if(SOILBCS=='NGDC')  then 
+          if(trim(SOILBCS)=='NGDC')  then 
              if(     F25Tag) call soil_para_high (nc,nr,regrid,fnameRst,F25Tag=F25Tag)
              if(.not.F25Tag) call soil_para_high (nc,nr,regrid,fnameRst)
           endif
-          if(SOILBCS=='HWSD') call soil_para_hwsd (nc,nr,fnameRst)
+          if(SOILBCS(1:4)=='HWSD') call soil_para_hwsd (nc,nr,fnameRst)
           write (log_file,'(a)')'         Done.'           
        else
           write (log_file,'(a,a)')'         Using existing file.'
@@ -586,8 +586,8 @@ PROGRAM mkCatchParam
        inquire(file=trim(fname_tmp4), exist=file_exists4)
        if ((.not.file_exists).or.(.not.file_exists2).or.(.not.file_exists3).or.(.not.file_exists4)) then
           write (log_file,'(a)')'         Creating files...'
-          if(SOILBCS=='NGDC') call create_model_para (MaskFile)
-          if(SOILBCS=='HWSD') call create_model_para_woesten (MaskFile) 
+          if(trim(SOILBCS)=='NGDC') call create_model_para(        MaskFile)
+          if(SOILBCS(1:4) =='HWSD') call create_model_para_woesten(MaskFile) 
           write (log_file,'(a)')'         Done.'           
        else
           write (log_file,'(a,a)')'         Using existing files.'

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mod_process_hres_data.F90

@@ -3608,7 +3608,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,fnameRst)
       character*100 :: fout
       character*200 :: fname
       character*10 :: string
-      character*2 :: VV,HH
+      character*2 :: VV,HH, tmpversion
 
       logical, allocatable, dimension(:,:) :: land_pixels
       integer, allocatable, dimension (:,:) :: &
@@ -3869,7 +3869,16 @@ SUBROUTINE soil_para_hwsd (nx,ny,fnameRst)
       ! get info common to all H[xx]V[yy] rectangles (could in theory differ from that
       !   of soildepth data read above but is the same as of 29 Apr 2022).
 
-      fname =trim(MAKE_BCS_INPUT_DIR)//'/land/soil/SOIL-DATA/soil_properties/v2/SoilProperties_H11V13.nc'
+      if      (trim(SOILBCS)=='HWSD_b')  then
+         tmpversion = 'v3'
+      else if (trim(SOILBCS)=='HWSD')    then
+         tmpversion = 'v2'
+      else
+         print *, 'Unknown SOILBCS: ', SOILBCS  
+         stop         
+      end if
+
+      fname =trim(MAKE_BCS_INPUT_DIR)//'/land/soil/SOIL-DATA/soil_properties/' // tmpversion // '/SoilProperties_H11V13.nc'
       status = NF_OPEN(trim(fname),NF_NOWRITE, ncid); VERIFY_(STATUS)
       !status = NF_GET_att_INT(ncid,NF_GLOBAL,'i_ind_offset_LL',iLL); VERIFY_(STATUS)  ! cannot be needed here
       !status = NF_GET_att_INT(ncid,NF_GLOBAL,'j_ind_offset_LL',jLL); VERIFY_(STATUS)  ! cannot be needed here
@@ -3925,7 +3934,7 @@ SUBROUTINE soil_para_hwsd (nx,ny,fnameRst)
       	 do ix = 1,36
 	    write (vv,'(i2.2)')jx
 	    write (hh,'(i2.2)')ix 
-	    fname = trim(MAKE_BCS_INPUT_DIR)//'/land/soil/SOIL-DATA/soil_properties/v2/SoilProperties_H'//hh//'V'//vv//'.nc'
+	    fname = trim(MAKE_BCS_INPUT_DIR)//'/land/soil/SOIL-DATA/soil_properties/' // tmpversion // '/SoilProperties_H'//hh//'V'//vv//'.nc'
             status = NF_OPEN(trim(fname),NF_NOWRITE, ncid)
 	    if(status == 0) then
 		status = NF_GET_att_INT  (ncid, NF_GLOBAL,'i_ind_offset_LL',iLL); VERIFY_(STATUS)

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90

@@ -87,8 +87,10 @@ SUBROUTINE init_bcs_config (LBCSV)
     !   MODC061   : Static snow albedo derived from MODIS Collection 6.1 data where available, fill value of 0.56 elsewhere. 
     !   MODC061v2 : Same as MODC061 but using tile ID instead of tile bounding box for mapping from raster to tile.
     !
-    ! SOILBCS: Soil parameter data.             DEFAULT : HWSD                                                       
-    !   HWSD      : Merged HWSD-STATSGO2 soil properties on 43200x21600 with Woesten et al. (1999) parameters   
+    ! SOILBCS: Soil parameter data.             DEFAULT : HWSD    
+    !   NGDC      : Soil parameters from Reynolds et al. 2000, doi:10.1029/2000WR900130 (MERRA-2, Fortuna, Ganymed, Icarus)
+    !   HWSD      : Merged HWSDv1.21-STATSGO2 soil properties on 43200x21600 with Woesten et al. (1999) parameters   
+    !   HWSD_b    : As in HWSD but with surgical fix of Argentina peatland issue (38S,60W)
 
     implicit none
 
@@ -194,6 +196,16 @@ SUBROUTINE init_bcs_config (LBCSV)
        GNU     = 1.0
        use_PEATMAP = .true.
        jpl_height  = .true.
+
+     case ("v12")   
+       LAIBCS  = 'MODGEO'
+       SOILBCS = 'HWSD_b'
+       MODALB  = 'MODIS2'
+       SNOWALB = 'MODC061v2'
+       GNU     = 1.0
+       use_PEATMAP = .true.
+       jpl_height  = .true.
+
     case default
 
        print *,'init_bcs_config(): unknown land boundary conditions version (LBCSV)'