Merge pull request #618 from GEOS-ESM/feature/borescan_snowalb_modis2

Add Option to Use MODIS-based Snow Albedo Climatology

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90

@@ -50,7 +50,7 @@ subroutine SetServices ( GC, RC )
     character(len=ESMF_MAXSTR)   :: CATCHCN_VERSION
     character(len=ESMF_MAXSTR)              :: SURFRC
     type(ESMF_Config)                       :: SCF
-    integer :: DO_GOSWIM, LSM_CHOICE, ATM_CO2
+    integer :: DO_GOSWIM, LSM_CHOICE, ATM_CO2, SNOW_ALBEDO_INFO
     character(len=ESMF_MAXSTR)              :: tmp
     integer                                 :: NUM_LDAS_ENSEMBLE, ens_id_width
 
@@ -82,6 +82,14 @@ subroutine SetServices ( GC, RC )
     SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS)
     call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS)
     call ESMF_ConfigGetAttribute (SCF, label='ATM_CO2:', value=ATM_CO2, DEFAULT=2, RC=STATUS) ; VERIFY_(STATUS)
+
+    ! SNOW ALBEDO -- so far, only parameterization based on look-up table is implemented for CatchCN
+    ! 0 : parameterization based on look-up table 
+    ! 1 : MODIS-derived snow albedo (backfilled with global land average snow albedo)
+    call ESMF_ConfigGetAttribute (SCF, label='SNOW_ALBEDO_INFO:', value=SNOW_ALBEDO_INFO, DEFAULT=0, RC=STATUS) ; VERIFY_(STATUS)
+
+    _ASSERT( SNOW_ALBEDO_INFO==0, "SNOW_ALBEDO_INFO must be 0 for CatchCN")
+
     call ESMF_ConfigGetAttribute (SCF, label='N_CONST_LAND4SNWALB:'  , value=DO_GOSWIM  , DEFAULT=0, RC=STATUS); VERIFY_(STATUS)
 
     if ( LSM_CHOICE == 2 ) then

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90

@@ -136,7 +136,8 @@ module GEOS_CatchGridCompMod
 end type CATCH_WRAP
 !#--
 
-integer :: USE_ASCATZ0, Z0_FORMULATION, AEROSOL_DEPOSITION, N_CONST_LAND4SNWALB,CHOOSEMOSFC
+integer :: USE_ASCATZ0, Z0_FORMULATION, AEROSOL_DEPOSITION, N_CONST_LAND4SNWALB
+integer :: CHOOSEMOSFC, SNOW_ALBEDO_INFO
 real    :: SURFLAY              ! Default (Ganymed-3 and earlier) SURFLAY=20.0 for Old Soil Params
                                 !         (Ganymed-4 and later  ) SURFLAY=50.0 for New Soil Params
 real    :: FWETC, FWETL
@@ -227,15 +228,21 @@ subroutine SetServices ( GC, RC )
        call MAPL_GetResource (SCF, FWETL, label='FWETL:', DEFAULT=0.025, __RC__ )
     endif
 
-    ! GOSWIM ANOW_ALBEDO 
+    ! SNOW ALBEDO 
+    ! 0 : parameterization based on look-up table 
+    ! 1 : MODIS-derived snow albedo (backfilled with global land average snow albedo)
+    call MAPL_GetResource (SCF, SNOW_ALBEDO_INFO,    label='SNOW_ALBEDO_INFO:',    DEFAULT=0, __RC__ )
+
+    ! GOSWIM SNOW_ALBEDO 
     ! 0 : GOSWIM snow albedo scheme is turned off
     ! 9 : i.e. N_CONSTIT in Stieglitz to turn on GOSWIM snow albedo scheme 
-    call MAPL_GetResource (SCF, N_CONST_LAND4SNWALB, label='N_CONST_LAND4SNWALB:', DEFAULT=0  , __RC__ )
+    call MAPL_GetResource (SCF, N_CONST_LAND4SNWALB, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ )
 
     ! 1: Use all GOCART aerosol values, 0: turn OFF everythying, 
     ! 2: turn off dust ONLY,3: turn off Black Carbon ONLY,4: turn off Organic Carbon ONLY
     ! __________________________________________
-    call MAPL_GetResource (SCF, AEROSOL_DEPOSITION, label='AEROSOL_DEPOSITION:', DEFAULT=0  , __RC__ )
+    call MAPL_GetResource (SCF, AEROSOL_DEPOSITION,  label='AEROSOL_DEPOSITION:',  DEFAULT=0, __RC__ )
+
     call ESMF_ConfigDestroy(SCF, __RC__)
 
 ! Set the Run entry points
@@ -1369,6 +1376,19 @@ subroutine SetServices ( GC, RC )
                                            RC=STATUS  ) 
   VERIFY_(STATUS)
 
+  if (SNOW_ALBEDO_INFO == 1) then
+    call MAPL_AddInternalSpec(GC                  ,&
+       LONG_NAME          = 'effective_snow_albedo'               ,&
+       UNITS              = '1'                         ,&
+       SHORT_NAME         = 'SNOWALB'                   ,&
+       FRIENDLYTO         = trim(COMP_NAME)             ,&
+       DIMS               = MAPL_DimsTileOnly           ,&
+       VLOCATION          = MAPL_VLocationNone          ,&
+       RESTART            = MAPL_RestartRequired        ,&
+                                           RC=STATUS  ) 
+     VERIFY_(STATUS)
+  endif
+
   call MAPL_AddInternalSpec(GC                  ,&
     LONG_NAME          = 'surface_heat_exchange_coefficient',&
     UNITS              = 'kg m-2 s-1'                ,&
@@ -3763,6 +3783,7 @@ subroutine Driver ( RC )
         real, dimension(:),   pointer :: psis
         real, dimension(:),   pointer :: bee
         real, dimension(:),   pointer :: poros
+        real, dimension(:),   pointer :: snowalb
         real, dimension(:),   pointer :: wpwet
         real, dimension(:),   pointer :: cond
         real, dimension(:),   pointer :: gnu
@@ -4130,6 +4151,7 @@ subroutine Driver ( RC )
         integer                       :: nv, nVars
         integer                       :: nDims,dimSizes(3)
         integer                       :: ldas_ens_id, ldas_first_ens_id
+
 !#---
 
         ! --------------------------------------------------------------------------
@@ -4350,8 +4372,8 @@ subroutine Driver ( RC )
         call MAPL_GetPointer(INTERNAL,CM         ,'CM'         ,RC=STATUS); VERIFY_(STATUS)
         call MAPL_GetPointer(INTERNAL,CQ         ,'CQ'         ,RC=STATUS); VERIFY_(STATUS)
         call MAPL_GetPointer(INTERNAL,FR         ,'FR'         ,RC=STATUS); VERIFY_(STATUS)
-        call MAPL_GetPointer(INTERNAL,DCQ        ,'DCQ'         ,RC=STATUS); VERIFY_(STATUS)
-        call MAPL_GetPointer(INTERNAL,DCH        ,'DCH'         ,RC=STATUS); VERIFY_(STATUS)
+        call MAPL_GetPointer(INTERNAL,DCQ        ,'DCQ'        ,RC=STATUS); VERIFY_(STATUS)
+        call MAPL_GetPointer(INTERNAL,DCH        ,'DCH'        ,RC=STATUS); VERIFY_(STATUS)
         if (N_CONST_LAND4SNWALB /= 0) then
            call MAPL_GetPointer(INTERNAL,RDU001     ,'RDU001'     , RC=STATUS); VERIFY_(STATUS)
            call MAPL_GetPointer(INTERNAL,RDU002     ,'RDU002'     , RC=STATUS); VERIFY_(STATUS)
@@ -4849,10 +4871,26 @@ subroutine Driver ( RC )
                  RHOFS,                                              &   
                  SNWALB_VISMAX, SNWALB_NIRMAX, SLOPE,                & 
                  WESNN, HTSNNN, SNDZN,                               &
-                 ALBVR, ALBNR, ALBVF, ALBNF, & ! instantaneous snow-free albedos on tiles
-                 SNOVR, SNONR, SNOVF, SNONF, &  ! instantaneous snow albedos on tiles
+                 ALBVR, ALBNR, ALBVF, ALBNF, &  ! instantaneous snow-free albedos on tiles
+                 SNOVR, SNONR, SNOVF, SNONF, &  ! instantaneous snow      albedos on tiles
                  RCONSTIT, UUU, TPSN1OUT1, DRPAR, DFPAR)    
 
+        if (SNOW_ALBEDO_INFO == 1) then
+
+           ! use MODIS-derived snow albedo from bcs (via Catch restart)
+           ! 
+           ! as a restart parameter from the bcs, snow albedo must not have no-data-values 
+           ! (checks for unphysical values should be in the make_bcs package)
+
+           call MAPL_GetPointer(INTERNAL,SNOWALB,'SNOWALB',RC=STATUS); VERIFY_(STATUS)
+
+           SNOVR = SNOWALB
+           SNONR = SNOWALB
+           SNOVF = SNOWALB
+           SNONF = SNOWALB
+
+        endif
+
         ! --------------------------------------------------------------------------
         ! albedo/swnet partitioning
         ! --------------------------------------------------------------------------
@@ -5531,6 +5569,20 @@ subroutine Driver ( RC )
                  SNOVR, SNONR, SNOVF, SNONF, & ! instantaneous snow albedos on tiles
                  RCONSTIT, UUU, TPSN1OUT1,DRPAR, DFPAR)   
 
+        if (SNOW_ALBEDO_INFO == 1) then
+
+           ! use MODIS-derived snow albedo from bcs (via Catch restart)
+           ! 
+           ! as a restart parameter from the bcs, snow albedo must not have no-data-values 
+           ! (checks for unphysical values should be in the make_bcs package)
+
+           SNOVR = SNOWALB
+           SNONR = SNOWALB
+           SNOVF = SNOWALB
+           SNONF = SNOWALB
+
+        endif
+
         ALBVR   = ALBVR    *(1.-ASNOW) + SNOVR    *ASNOW
         ALBVF   = ALBVF    *(1.-ASNOW) + SNOVF    *ASNOW
         ALBNR   = ALBNR    *(1.-ASNOW) + SNONR    *ASNOW

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc

@@ -114,6 +114,16 @@
 #
 # GEOSagcm=>MODIS_DVG: 0                     # NOT yet used in GEOSldas
 
+# ---- SNOW ALBEDO  
+#
+#    0  : Snow albedo parameterization based on look-up table (default)
+#    1  : Snow albedo derived from MODIS Collection MOD10A1.061 (Feb/2000 - Mar/2022)  
+#           - backfilled with global land average snow albedo where unavailable
+#           - must use compatible bcs version (e.g., v06, v08, v09)
+#           - NOT YET AVAILABLE FOR CatchCN
+#
+# GEOSagcm=>SNOW_ALBEDO_INFO: 0
+# GEOSldas=>SNOW_ALBEDO_INFO: 0
 
 #--------------------------------------------------------#
 #    GOSWIM aerosol deposition on surface snow           #
@@ -234,3 +244,5 @@
 #
 # GEOSagcm=>PRESCRIBE_DVG: 0
 # GEOSldas=>PRESCRIBE_DVG: 0
+
+

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs

@@ -103,13 +103,19 @@ if ( $HELPMODE != YES ) then
 else
    echo "Options for land BCs version:"
 endif
-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 "   ${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}      = NL3 + JPL veg height{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}      = NL3 + JPL veg height + PEATMAP${CR}"
+echo "   ${C2}v06   = NL3 + JPL veg height + PEATMAP + MODIS snow alb${CR}"
+echo "   ${C2}v07   = NL3 + PEATMAP${CR}"
+echo "   ${C2}v08   = NL3 + MODIS snow alb${CR}"
+echo "   ${C2}v09   = NL3 + PEATMAP + MODIS snow alb${CR}"
 echo "  "
 if ( $HELPMODE != YES ) then        
    echo "                         ${C1}       *BCs produced by this code will differ from BCs in archived directories\!\!\! ${CR}"
@@ -128,16 +134,18 @@ if ( $HELPMODE != YES ) then
 
    set dummy = `echo $<`
 
-   set dummy = `echo $dummy | tr "[:lower:]" "[:upper:]"`
 
    if( $dummy == 'F25'  | \
        $dummy == 'GM4'  | \
        $dummy == 'ICA'  | \
        $dummy == 'NL3'  | \
        $dummy == 'NL4'  | \
        $dummy == 'NL5'  | \
-       $dummy == 'DEV') then
-      set lbcsv        = $dummy
+       $dummy == 'v06'  | \
+       $dummy == 'v07'  | \
+       $dummy == 'v08'  | \
+       $dummy == 'v09') then
+      set lbcsv = $dummy
    else if ( $dummy == '' ) then
       echo $lbcsv
    else

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

@@ -4,14 +4,14 @@ PROGRAM mkCatchParam
 !
 ! !ARGUMENTS:
 !
-!  Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBSV -e EASE"       
+!  Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE"       
 !     -x: Size of longitude dimension of input raster. DEFAULT: 8640
 !     -y: Size of latitude dimension of input raster.  DEFAULT: 4320
 !     -b: position of the dateline in the first box. DEFAULT: DC 
 !     -g: Gridname  (name of the .til or .rst file without file extension)  
-!     -v: LBCSV : Choose bcs version (ICA, NL3, NL4, NL5, or development)             
-!     -e: EASE : This is optional if catchment.def file is available already or                    
-!         the til file format is pre-Fortuna-2.                                                    
+!     -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09)
+!     -e: EASE  : This is optional if catchment.def file is available already or                    
+!                  the til file format is pre-Fortuna-2.                                                    
 !     
 !
 ! This program is good to generate  
@@ -21,7 +21,7 @@ PROGRAM mkCatchParam
 ! Sarith Mahanama - March 23, 2012 
 ! Email: sarith.p.mahanama@nasa.gov
 
-  use rmTinyCatchParaMod
+  use rmTinyCatchParaMod 
   use process_hres_data
   !   use module_irrig_params, ONLY : create_irrig_params
 
@@ -33,7 +33,7 @@ PROGRAM mkCatchParam
   ! NC and NR are typically overwritten through command-line arguments "-x nx -y ny".
 
   integer              :: NC = i_raster, NR = j_raster    
-  character*4          :: LBSV = 'DEF'
+  character*5          :: LBCSV = 'UNDEF'
   character*128        :: GridName = ''
   character*128        :: ARG, MaskFile
   character*256        :: CMD
@@ -60,6 +60,7 @@ PROGRAM mkCatchParam
   character*200        :: tmpstring, tmpstring1, tmpstring2
   character*200        :: fname_tmp, fname_tmp2, fname_tmp3, fname_tmp4
   integer              :: N_tile
+  logical              :: process_snow_albedo = .false. 
 
 ! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------
 !
@@ -102,14 +103,14 @@ PROGRAM mkCatchParam
 !   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                                     "
-    USAGE(3) ="     -y: Size of latitude dimension of input raster.  DEFAULT: 4320                                     "
-    USAGE(4) ="     -g: Gridname  (name of the .til or .rst file without file extension)                               "
-    USAGE(5) ="     -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC                         "
-    USAGE(6) ="     -e: EASE : This is optional if catchment.def file is available already or                          "          
-    USAGE(7) ="                the til file format is pre-Fortuna-2.                                                   "
-    USAGE(8) ="     -v  LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, or DEV)                              "
+    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                   "
+    USAGE(3) ="     -y: Size of latitude dimension of input raster.  DEFAULT: 4320                   "
+    USAGE(4) ="     -g: Gridname  (name of the .til or .rst file without file extension)             "
+    USAGE(5) ="     -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC       "
+    USAGE(6) ="     -e: EASE  : This is optional if catchment.def file is available already or       "          
+    USAGE(7) ="                  the til file format is pre-Fortuna-2.                               "
+    USAGE(8) ="     -v  LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09)  "
 
 ! Process Arguments                            
 !------------------ 
@@ -153,9 +154,9 @@ PROGRAM mkCatchParam
        case ('g')
           GridName = trim(arg)
        case ('v')
-          LBSV = trim(arg)
+          LBCSV = trim(arg)
           if (trim(arg).eq."F25") F25Tag = .true.
-          call init_bcs_config (trim(LBSV))
+          call init_bcs_config (trim(LBCSV))       ! get bcs details from version string
        case ('b')
           DL = trim(arg)
        case ('e')
@@ -191,11 +192,14 @@ PROGRAM mkCatchParam
     if(use_PEATMAP)  PEATSOURCE   = 'PEATMAP'
     if(jpl_height)   VEGZSOURCE   = 'JPL'
 
+    if (trim(SNOWALB)=='MODC061')  process_snow_albedo=.true.
+
     if(n_threads == 1) then
 
        write (log_file,'(a)')trim(LAIBCS)
        write (log_file,'(a)')trim(MODALB)
        write (log_file,'(a)')trim(SOILBCS)   
+       write (log_file,'(a)')trim(SNOWALB)
        write (log_file,'(a)')trim(MaskFile)
        write (log_file,'(a)')trim(PEATSOURCE)
        write (log_file,'(a)')trim(VEGZSOURCE)
@@ -249,7 +253,7 @@ PROGRAM mkCatchParam
         close (10, status = 'keep')
 
        inquire(file='clsm/catch_params.nc4', exist=file_exists)
-       if (.not.file_exists) CALL open_landparam_nc4_files(N_tile)  
+       if (.not.file_exists) CALL open_landparam_nc4_files(N_tile,process_snow_albedo)  
 
        ! Creating cti_stats.dat 
        ! ----------------------
@@ -564,7 +568,7 @@ PROGRAM mkCatchParam
              if(     F25Tag) call soil_para_high (nc,nr,regrid,gridnamer,F25Tag=F25Tag)
              if(.not.F25Tag) call soil_para_high (nc,nr,regrid,gridnamer)
           endif
-          if(SOILBCS=='HWSD')  call soil_para_hwsd (nc,nr,gridnamer)
+          if(SOILBCS=='HWSD') call soil_para_hwsd (nc,nr,gridnamer)
           write (log_file,'(a)')'         Done.'           
        else
           write (log_file,'(a,a)')'         Using existing file.'
@@ -653,6 +657,15 @@ PROGRAM mkCatchParam
        endif
        write (log_file,'(a)')' '
 
+       if(process_snow_albedo)then
+          tmpstring = 'Step 14: Static snow albedo from MODIS' 
+          write (log_file,'(a)') trim(tmpstring)
+          write (log_file,'(a)')'         Creating file...'
+          call MODIS_snow_alb ( )
+          write (log_file,'(a)')'         Done.'           
+          write (log_file,'(a)')' '
+       endif 
+
        !      inquire(file='clsm/irrig.dat', exist=file_exists)
        !      if (.not.file_exists) call create_irrig_params (nc,nr,gridnamer)
        !      write (log_file,'(a)')'Done computing irrigation model parameters ...............13'

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90

@@ -16,7 +16,7 @@ PROGRAM mkSMAPTilesPara_v2
       use LogRectRasterizeMod
 
       implicit none
-      character*4          :: LBSV = 'DEF'
+      character*5          :: LBCSV = 'UNDEF'
       integer i,j,ig,jg,i0,iop,n,d1,d2,j1,j2,i1,i2,ix, jx,icount,pcount
       integer :: NC = i_raster, NR = j_raster, NT = 16330000, ND = 10000, ND_raster = 10000
 
@@ -95,7 +95,7 @@ PROGRAM mkSMAPTilesPara_v2
 
          elseif ( trim(arg) == '-v' ) then
             i = i+1
-            call get_command_argument(i,LBSV)
+            call get_command_argument(i,LBCSV)
 
          else ! stop for any other arguments
 
@@ -753,7 +753,7 @@ PROGRAM mkSMAPTilesPara_v2
       ! now run mkCatchParam
       ! --------------------
 
-      tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBSV)
+      tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBCSV)
       write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr
       tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1)
       print *,trim(tmpstring)