diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 index 2c290094f..d9596569e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 @@ -5167,13 +5167,10 @@ subroutine Driver ( RC ) ! Variables for FPAR ! -------------------------- real , allocatable, dimension (:,:) :: parzone + character(len=ESMF_MAXSTR) :: Co2_CycleFile IAm=trim(COMP_NAME)//"::RUN2::Driver" - ! Begin - - IAm=trim(COMP_NAME)//"Driver" - ! -------------------------------------------------------------------------- ! Get time step from configuration ! -------------------------------------------------------------------------- @@ -5656,7 +5653,11 @@ subroutine Driver ( RC ) call MPI_Info_create(info, STATUS); VERIFY_(status) call MPI_Info_set(info, "romio_cb_read", "automatic", STATUS); VERIFY_(status) - STATUS = NF_OPEN ('CO2_MonthlyMean_DiurnalCycle.nc4', NF_NOWRITE, CTfile); VERIFY_(status) + call MAPL_GetResource (MAPL, CO2_CycleFile, label = 'CO2_MonthlyMean_DiurnalCycle_FILE:', & + default = 'CO2_MonthlyMean_DiurnalCycle.nc4', RC=STATUS ) + VERIFY_(STATUS) + + STATUS = NF_OPEN (trim(CO2_CycleFile), NF_NOWRITE, CTfile); VERIFY_(status) allocate (CT_CO2V (1: NUNQ, 1:12, 1:8)) allocate (CTCO2_TMP (1:CT_grid_N_lon, 1:CT_grid_N_lat, 1:12, 1:8)) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 index 9e380d22a..5fa098a89 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 @@ -5115,7 +5115,7 @@ subroutine Driver ( RC ) real, allocatable, dimension(:,:,:) :: pft real, allocatable, dimension(:) :: lnfm - character(len=ESMF_MAXSTR) :: LNFMFile + character(len=ESMF_MAXSTR) :: LNFMFile, CO2_CycleFile integer :: ntile, nv, dpy, ierr, iok, ndt integer, save :: year_prev = -9999 @@ -5653,7 +5653,9 @@ subroutine Driver ( RC ) call MPI_Info_create(info, STATUS); VERIFY_(status) call MPI_Info_set(info, "romio_cb_read", "automatic", STATUS); VERIFY_(status) - STATUS = NF_OPEN ('CO2_MonthlyMean_DiurnalCycle.nc4', NF_NOWRITE, CTfile); VERIFY_(status) + call MAPL_GetResource (MAPL, CO2_CycleFile, label = 'CO2_MonthlyMean_DiurnalCycle_FILE:', default = 'CO2_MonthlyMean_DiurnalCycle.nc4', RC=STATUS ) + VERIFY_(STATUS) + STATUS = NF_OPEN (trim(CO2_CycleFile), NF_NOWRITE, CTfile); VERIFY_(status) allocate (CT_CO2V (1: NUNQ, 1:12, 1:8)) allocate (CTCO2_TMP (1:CT_grid_N_lon, 1:CT_grid_N_lat, 1:12, 1:8)) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt index fc6ef7bb5..add425bac 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt @@ -1,50 +1 @@ -esma_set_this(OVERRIDE raster) -set (srcs -date_time_util.F90 -leap_year.F90 -EASE_conv.F90 -mod_process_hres_data.F90 -rasterize.F90 -read_riveroutlet.F90 -CubedSphere_GridMod.F90 -rmTinyCatchParaMod.F90 -zip.c -util.c -) - -if(NOT FORTRAN_COMPILER_SUPPORTS_FINDLOC) - list(APPEND srcs findloc.F90) -endif () - -set_source_files_properties(mkMITAquaRaster.F90 PROPERTIES COMPILE_FLAGS "${BYTERECLEN}") - -esma_add_library(${this} SRCS ${srcs} DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared esmf NetCDF::NetCDF_Fortran OpenMP::OpenMP_Fortran OpenMP::OpenMP_C) - -if(NOT FORTRAN_COMPILER_SUPPORTS_FINDLOC) - target_compile_definitions(${this} PRIVATE USE_EXTERNAL_FINDLOC) -endif () - -# MAT NOTE This should use find_package(ZLIB) but Baselibs currently -# confuses find_package(). This is a hack until Baselibs is -# reorganized. -if (Baselibs_FOUND) - set (INC_ZLIB ${BASEDIR}/include/zlib) - target_include_directories(${this} PRIVATE ${INC_ZLIB}) -else () - find_package(ZLIB) - target_link_libraries(${this} PRIVATE ZLIB::zlib) -endif () - -ecbuild_add_executable (TARGET chk_clsm_params.x SOURCES chk_clsm_params.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET CombineRasters.x SOURCES CombineRasters.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkCatchParam.x SOURCES mkCatchParam.F90 LIBS MAPL ${this} OpenMP::OpenMP_Fortran) -ecbuild_add_executable (TARGET mkCubeFVRaster.x SOURCES mkCubeFVRaster.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkLandRaster.x SOURCES mkLandRaster.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkLatLonRaster.x SOURCES mkLatLonRaster.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkMITAquaRaster.x SOURCES mkMITAquaRaster.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkMOMAquaRaster.x SOURCES mkMOMAquaRaster.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET FillMomGrid.x SOURCES FillMomGrid.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mk_runofftbl.x SOURCES mk_runofftbl.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkEASETilesParam.x SOURCES mkEASETilesParam.F90 LIBS MAPL ${this}) - -install(PROGRAMS make_bcs clsm_plots.pro plot_curves.pro create_README.csh plot_curves.csh DESTINATION bin) +esma_add_subdirectories (makebcs preproc) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/Raster.h b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/Raster.h deleted file mode 100644 index 90a7e04d6..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/Raster.h +++ /dev/null @@ -1,6 +0,0 @@ - -#define VERIFY_(A) IF(A/=0)THEN;PRINT *,'ERROR AT LINE ', __LINE__;STOP;ENDIF -#define _ASSERT(A,msg) if(.not.A)then;print *,'Error:',__FILE__,__LINE__;stop;endif -#define REAL_ real(kind=8) -#define RASTER_PI 3.14159265358979323846264338_8 -#define RASTERUNDEF -999 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/leap_year.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/leap_year.F90 deleted file mode 100644 index 5382a570f..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/leap_year.F90 +++ /dev/null @@ -1,95 +0,0 @@ - -module leap_year - - implicit none - -contains - - integer function days_in_month(year, month) - - ! return the number of days in a given month - - implicit none - - integer :: year, month - - integer, dimension(12), parameter :: days_in_month_leap = & - (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) - - integer, dimension(12), parameter :: days_in_month_nonleap = & - (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) - - if (is_leap_year(year)) then - days_in_month = days_in_month_leap(month) - else - days_in_month = days_in_month_nonleap(month) - end if - - end function days_in_month - - ! ------------------------------------------------------------------ - - integer function days_in_year(year) - - ! return the number of days in a given year - - implicit none - - integer :: year - - if (is_leap_year(year)) then - days_in_year = 366 - else - days_in_year = 365 - end if - - end function days_in_year - - ! ------------------------------------------------------------------ - - logical function is_leap_year(year) - - implicit none - - integer :: year - - if (mod(year,4) /= 0) then - is_leap_year = .false. - else if (mod(year,400) == 0) then - is_leap_year = .true. - else if (mod(year,100) == 0) then - is_leap_year = .false. - else - is_leap_year = .true. - end if - - end function is_leap_year - - ! ------------------------------------------------------------------ - - integer function pentad_of_year(day_of_year, year) - - implicit none - - integer :: day_of_year, year - - ! determine pentad - - if ((is_leap_year(year)) .and. day_of_year>=59) then - - pentad_of_year = (day_of_year-2)/5+1 - - else - - pentad_of_year = (day_of_year-1)/5+1 - - end if - - end function pentad_of_year - -end module leap_year - - -! ========= EOF ========================================================= - - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt new file mode 100755 index 000000000..6d7b0328e --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt @@ -0,0 +1,49 @@ +esma_set_this() + +set (srcs +LDAS_DateTimeMod.F90 +EASE_conv.F90 +mod_process_hres_data.F90 +rasterize.F90 +read_riveroutlet.F90 +CubedSphere_GridMod.F90 +rmTinyCatchParaMod.F90 +zip.c +util.c +) + +if(NOT FORTRAN_COMPILER_SUPPORTS_FINDLOC) + list(APPEND srcs findloc.F90) +endif () + +set_source_files_properties(mkMITAquaRaster.F90 PROPERTIES COMPILE_FLAGS "${BYTERECLEN}") + +esma_add_library(${this} SRCS ${srcs} DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared esmf NetCDF::NetCDF_Fortran OpenMP::OpenMP_Fortran OpenMP::OpenMP_C) + +if(NOT FORTRAN_COMPILER_SUPPORTS_FINDLOC) + target_compile_definitions(${this} PRIVATE USE_EXTERNAL_FINDLOC) +endif () + +# MAT NOTE This should use find_package(ZLIB) but Baselibs currently +# confuses find_package(). This is a hack until Baselibs is +# reorganized. +if (Baselibs_FOUND) + set (INC_ZLIB ${BASEDIR}/include/zlib) + target_include_directories(${this} PRIVATE ${INC_ZLIB}) +else () + find_package(ZLIB) + target_link_libraries(${this} PRIVATE ZLIB::zlib) +endif () + +ecbuild_add_executable (TARGET CombineRasters.x SOURCES CombineRasters.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mkCatchParam.x SOURCES mkCatchParam.F90 LIBS MAPL ${this} OpenMP::OpenMP_Fortran) +ecbuild_add_executable (TARGET mkCubeFVRaster.x SOURCES mkCubeFVRaster.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mkLandRaster.x SOURCES mkLandRaster.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mkLatLonRaster.x SOURCES mkLatLonRaster.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mkMITAquaRaster.x SOURCES mkMITAquaRaster.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mkMOMAquaRaster.x SOURCES mkMOMAquaRaster.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET FillMomGrid.x SOURCES FillMomGrid.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mk_runofftbl.x SOURCES mk_runofftbl.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mkEASETilesParam.x SOURCES mkEASETilesParam.F90 LIBS MAPL ${this}) + +install(PROGRAMS make_bcs clsm_plots.pro create_README.csh DESTINATION bin) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CombineRasters.F90 similarity index 96% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CombineRasters.F90 index a5ba74ffc..3c3605c27 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CombineRasters.F90 @@ -1,12 +1,13 @@ -! $Id: - -#include "Raster.h" +#define I_AM_MAIN +#include "MAPL_ErrLog.h" program mkOverlaySimple use LogRectRasterizeMod use MAPL_SortMod use MAPL_HashMod + use MAPL_ExceptionHandling + use MAPL_Constants ! Overlay atmosphere, land, and ocean rasters, creating a .idx file. ! The ocean raster should be defined everywhere, or at least, everywhere @@ -27,7 +28,7 @@ program mkOverlaySimple integer, parameter :: TILUNIT1 = 22 integer, parameter :: TILUNIT2 = 23 - REAL_, parameter :: PI = RASTER_PI + real(kind=8), parameter :: PI = MAPL_PI_R8 integer :: command_argument_count integer :: nxt, argl, fill @@ -42,12 +43,12 @@ program mkOverlaySimple integer, allocatable :: RST2(: ) integer, allocatable :: iTable(:,:) - REAL_ , allocatable :: Table1(:,:) - REAL_ , allocatable :: Table2(:,:) - REAL_ , allocatable :: rTable(:,:) - REAL_ , allocatable :: cc(:), ss(:) - REAL_ :: dx, dy, area, xc, yc, d2r, vv(4) - REAL_ :: lats, lons, da + real(kind=8) , allocatable :: Table1(:,:) + real(kind=8) , allocatable :: Table2(:,:) + real(kind=8) , allocatable :: rTable(:,:) + real(kind=8) , allocatable :: cc(:), ss(:) + real(kind=8) :: dx, dy, area, xc, yc, d2r, vv(4) + real(kind=8) :: lats, lons, da logical :: DoZip logical :: Verb @@ -64,6 +65,7 @@ program mkOverlaySimple character*128 :: & Usage = "CombineRasters -v -h -z -t MT -g GF -f TYPE BOTTOMRASTER TOPRASTER" + character*128 :: Iam = "CombineRasters" integer :: Pix1, Pix2 ! Argument defaults diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CubedSphere_GridMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CubedSphere_GridMod.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CubedSphere_GridMod.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CubedSphere_GridMod.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/EASE_conv.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_conv.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/EASE_conv.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_conv.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/FillMomGrid.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/FillMomGrid.F90 similarity index 99% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/FillMomGrid.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/FillMomGrid.F90 index 3b596641f..36dbc1921 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/FillMomGrid.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/FillMomGrid.F90 @@ -1,4 +1,3 @@ -!#include "Raster.h" #define VERIFY_(A) IF(A/=0)THEN;PRINT *,'ERROR AT LINE ', __LINE__;STOP;ENDIF #define ASSERT_(A) if(.not.A)then;print *,'Error:',__FILE__,__LINE__;stop;endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/date_time_util.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/LDAS_DateTimeMod.F90 similarity index 61% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/date_time_util.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/LDAS_DateTimeMod.F90 index 9876e7236..8799ad792 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/date_time_util.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/LDAS_DateTimeMod.F90 @@ -1,7 +1,5 @@ -module date_time_util - - use leap_year +module LDAS_DateTimeMod implicit none @@ -10,6 +8,9 @@ module date_time_util private public :: date_time_type + public :: date_time_print + public :: is_leap_year + public :: days_in_month public :: get_dofyr_pentad public :: augment_date_time public :: datetime2_minus_datetime1 @@ -17,6 +18,8 @@ module date_time_util public :: datetime_le_refdatetime public :: datetime_lt_refdatetime public :: date_time2string + public :: datetime_to_J2000seconds + public :: J2000seconds_to_datetime ! --------------------------------------------------------------------- @@ -39,10 +42,106 @@ module date_time_util integer :: pentad ! pentad of year integer :: dofyr ! day of year end type date_time_type - + contains - ! ********************************************************** + + function date_time_print(dt) result (dtstr) + + implicit none + + type(date_time_type), intent(in) :: dt + character(len=19) :: dtstr ! output + + write(dtstr,298) dt%year, dt%month, dt%day, dt%hour, dt%min, dt%sec +298 format(i4.4,'-',i2.2,'-',i2.2,'T',i2.2,':',i2.2,':',i2.2) + + end function date_time_print + ! ********************************************************** + + integer function days_in_month(year, month) + + ! return the number of days in a given month + + implicit none + + integer :: year, month + + integer, dimension(12), parameter :: days_in_month_leap = & + (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) + + integer, dimension(12), parameter :: days_in_month_nonleap = & + (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) + + if (is_leap_year(year)) then + days_in_month = days_in_month_leap(month) + else + days_in_month = days_in_month_nonleap(month) + end if + + end function days_in_month + + ! ------------------------------------------------------------------ + + integer function days_in_year(year) + + ! return the number of days in a given year + + implicit none + + integer :: year + + if (is_leap_year(year)) then + days_in_year = 366 + else + days_in_year = 365 + end if + + end function days_in_year + + ! ------------------------------------------------------------------ + + logical function is_leap_year(year) + + implicit none + + integer :: year + + if (mod(year,4) /= 0) then + is_leap_year = .false. + else if (mod(year,400) == 0) then + is_leap_year = .true. + else if (mod(year,100) == 0) then + is_leap_year = .false. + else + is_leap_year = .true. + end if + + end function is_leap_year + + ! ------------------------------------------------------------------ + + integer function pentad_of_year(day_of_year, year) + + implicit none + + integer :: day_of_year, year + + ! determine pentad + + if ((is_leap_year(year)) .and. day_of_year>=59) then + + pentad_of_year = (day_of_year-2)/5+1 + + else + + pentad_of_year = (day_of_year-1)/5+1 + + end if + + end function pentad_of_year + + ! ------------------------------------------------------------------ subroutine get_dofyr_pentad( date_time ) @@ -254,12 +353,14 @@ integer function datetime2_minus_datetime1( d1, d2 ) type(date_time_type) :: de, dl - integer :: fac, secs, tmpint, y, secs_in_year_de, secs_in_year_dl + integer :: fac, secs, y, secs_in_year_de, secs_in_year_dl ! ------------------------------------------------------------------- ! ! make sure type integer is not out of range - + ! + ! [integer*4 only allows for differences of up to ~68 years] + if ( (abs(d2%year-d1%year)+1) > ((huge(secs)/86400)/366) ) then write (*,*) 'datetime2_minus_datetime1(): integer out of range.' @@ -485,10 +586,8 @@ end function datetime2_gt_eq_lt_datetime1 character(16) function date_time2string( date_time ) - ! Generates a string from date_time structure (ignore seconds) - - !!use driver_types - + ! Generate "YYYYMMDD_HHMMSSz" string from date_time structure + implicit none type(date_time_type) :: date_time @@ -498,22 +597,169 @@ character(16) function date_time2string( date_time ) character(2) :: char_day character(2) :: char_hour character(2) :: char_min + character(2) :: char_sec write(char_year, '(i4.4)') date_time%year write(char_month, '(i2.2)') date_time%month write(char_day, '(i2.2)') date_time%day write(char_hour, '(i2.2)') date_time%hour write(char_min, '(i2.2)') date_time%min + write(char_sec, '(i2.2)') date_time%sec - date_time2string = char_year // '.' // char_month // '.' & - // char_day // '.' // char_hour // '.' // char_min + date_time2string = char_year // char_month // char_day // '_' & + // char_hour // char_min // char_sec // 'z' end function date_time2string ! ******************************************************************** -end module date_time_util + real*8 function datetime_to_J2000seconds( date_time, epoch_id ) + + ! reichle, 14 Jan 2014 + + implicit none + + type(date_time_type) :: date_time + + character(4) :: epoch_id + + integer :: tmp_secs + + ! -------------------------------------------- + + ! get integer seconds + + tmp_secs = datetime2_minus_datetime1( J2000_epoch(epoch_id), date_time ) + + ! convert to double precision floating point value + + datetime_to_J2000seconds = real(tmp_secs,kind(0.0D0)) + + ! correct for 0.816 milliseconds (if using "UT12" definition of J2000 Epoch) + + if (epoch_id=='UT12') datetime_to_J2000seconds = datetime_to_J2000seconds - 0.816D0 + + end function datetime_to_J2000seconds + + ! ******************************************************************** + + type(date_time_type) function J2000seconds_to_datetime( J2000_seconds, epoch_id ) + + ! reichle, 14 Jan 2014 + + implicit none + + real*8 :: J2000_seconds + + type(date_time_type) :: date_time + + character(4) :: epoch_id + + integer :: tmp_secs + + ! -------------------------------------------- + + date_time = J2000_epoch( epoch_id ) + + tmp_secs = nint(J2000_seconds) + + if ( J2000_seconds > real(huge(tmp_secs),kind(0.0D0)) ) then + + write (*,*) 'J2000seconds_to_datetime(): J2000_seconds out of range.' + write (*,*) 'STOPPING.' + stop + + end if + + call augment_date_time( tmp_secs, date_time ) + + J2000seconds_to_datetime = date_time + + end function J2000seconds_to_datetime + + ! ******************************************************************** + + type(date_time_type) function J2000_epoch( epoch_id ) + + ! reichle, 30 Jun 2015 + + implicit none + + character(4) :: epoch_id + + character(len=*), parameter :: Iam = 'J2000_epoch' + + ! ------------------------------------- + ! + ! definition of J2000 epochs + ! + ! "J2000 seconds" are elapsed seconds since J2000 Epoch, which is either + ! + ! - "UT12": 11:58:55.816 on 1 Jan 2000 in Coordinated Universal Time (UTC), or + ! - "TT12": 12:00:00.000 on 1 Jan 2000 in Terrestrial Time (TT), or + ! - "UT00": 00:00:00.000 on 1 Jan 2000 in Coordinated Universal Time (UTC) + ! + ! NOTE: Per SMAP L1C_TB data products specs document, SMAP time stamps use "UT12" + ! but sample granules appear to be using "TT12". + ! NOTE: Per Clara Draper (30 Jun 2015), the nc4 ASCAT soil moisture retrieval + ! product uses "UT00". + + type(date_time_type), parameter :: J2000_UT12 = date_time_type( & + year = 2000, & + month = 1, & + day = 1, & + hour = 11, & + min = 58, & + sec = 55, & ! rounded down + pentad = 1, & + dofyr = 1 ) + + type(date_time_type), parameter :: J2000_TT12 = date_time_type( & + year = 2000, & + month = 1, & + day = 1, & + hour = 12, & + min = 0, & + sec = 0, & + pentad = 1, & + dofyr = 1 ) + + type(date_time_type), parameter :: J2000_UT00 = date_time_type( & + year = 2000, & + month = 1, & + day = 1, & + hour = 0, & + min = 0, & + sec = 0, & + pentad = 1, & + dofyr = 1 ) + + ! ---------------------------- + + select case (epoch_id) + + case ('UT12') + J2000_epoch = J2000_UT12 + + case ('TT12') + J2000_epoch = J2000_TT12 + + case ('UT00') + J2000_epoch = J2000_UT00 + + case default + + write (*,*) Iam // ' unknown J2000 epoch_id' + write (*,*) 'STOPPING.' + stop + + end select + + end function J2000_epoch + + ! ******************************************************************** +end module LDAS_DateTimeMod ! ****************************************************************** @@ -523,14 +769,19 @@ end module date_time_util program test - use date_time_util + use LDAS_DateTimeMod implicit none - type(date_time_type) :: start_time, end_time + type(date_time_type) :: start_time, end_time, date_time, date_time_tmp integer :: diff + real*8 :: J2000_seconds = -987303600.0D0 + + character(4) :: J2000_epoch_id + + start_time%year = 1992 ! 4-digit year start_time%month = 11 ! month in year start_time%day = 1 ! day in month @@ -540,21 +791,31 @@ program test start_time%pentad = -9999 ! pentad of year start_time%dofyr = -9999 ! day of year - - end_time%year = 1998 ! 4-digit year - end_time%month = 11 ! month in year - end_time%day = 3 ! day in month - end_time%hour = 0 ! hour of day - end_time%min = 0 ! minute of hour - end_time%sec = 0 ! seconds of minute + end_time%year = 2000 ! 4-digit year + end_time%month = 1 ! month in year + end_time%day = 1 ! day in month + end_time%hour = 11 ! hour of day + end_time%min = 58 ! minute of hour + end_time%sec = 55 ! seconds of minute end_time%pentad = -9999 ! pentad of year end_time%dofyr = -9999 ! day of year + date_time%year = 2014 ! 4-digit year + date_time%month = 1 ! month in year + date_time%day = 14 ! day in month + date_time%hour = 12 ! hour of day + date_time%min = 34 ! minute of hour + date_time%sec = 56 ! seconds of minute + date_time%pentad = -9999 ! pentad of year + date_time%dofyr = -9999 ! day of year + + write (*,*) huge(diff) call get_dofyr_pentad( start_time ) call get_dofyr_pentad( end_time ) + call get_dofyr_pentad( date_time ) write (*,*) start_time write (*,*) end_time @@ -605,7 +866,51 @@ program test write (*,*) start_time + ! ------------------------------ + + J2000_epoch_id = 'UT12' + + write (*,*) '-------------------------------' + + write (*,*) 'J2000_epoch_id = ', J2000_epoch_id + + write (*,*) '-------------------------------' + write (*,*) start_time + write (*,*) datetime_to_J2000seconds( start_time, J2000_epoch_id ) + write (*,*) '-------------------------------' + write (*,*) end_time + write (*,*) datetime_to_J2000seconds( end_time, J2000_epoch_id ) + write (*,*) '-------------------------------' + write (*,*) date_time + write (*,*) datetime_to_J2000seconds( date_time, J2000_epoch_id ) + write (*,*) '-------------------------------' + + ! ------------------------------ + + date_time_tmp = end_time + + write (*,*) end_time + write (*,*) date_time_tmp + write (*,*) date_time + + diff = datetime2_minus_datetime1( end_time, date_time ) + + write (*,*) diff + write (*,*) datetime_to_J2000seconds( date_time, J2000_epoch_id ) + + call augment_date_time( diff, date_time_tmp ) + + write (*,*) date_time_tmp + ! ---------------------------------- + + write (*,*) '-------------------------------' + write (*,*) J2000_seconds + + write (*,*) J2000seconds_to_datetime( J2000_seconds, J2000_epoch_id ) + write (*,*) '-------------------------------' + + end program test #endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/bcs_utils.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/bcs_utils.py similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/bcs_utils.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/bcs_utils.py diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/clsm_plots.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/clsm_plots.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/clsm_plots.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/clsm_plots.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/create_README.csh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/create_README.csh similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/create_README.csh rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/create_README.csh diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/findloc.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/findloc.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/findloc.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/findloc.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs.py similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs.py diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_cube_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_cube_bcs.py similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_cube_bcs.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_cube_bcs.py diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_ease_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_ease_bcs.py similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_ease_bcs.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_ease_bcs.py diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_latlon_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_latlon_bcs.py similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_latlon_bcs.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_latlon_bcs.py diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCubeFVRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCubeFVRaster.F90 similarity index 96% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCubeFVRaster.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCubeFVRaster.F90 index 02c7c3c90..43b9646c7 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCubeFVRaster.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCubeFVRaster.F90 @@ -1,6 +1,5 @@ -! $Id: - -#include "Raster.h" +#define I_AM_MAIN +#include "MAPL_ErrLog.h" program mkCubeFVRaster @@ -18,7 +17,7 @@ program mkCubeFVRaster ! use CubedSphere_GridMod use LogRectRasterizeMod - + use MAPL_ExceptionHandling !EOP implicit none @@ -36,6 +35,7 @@ program mkCubeFVRaster logical :: Here=.false. logical :: Verb=.false. character*128 :: Usage="mkCubeFVraster -x RX -y RY -z -h -v -g GN ncells" + character*128 :: Iam ="mkCubeFVraster" ! Process Arguments !------------------ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 similarity index 99% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkEASETilesParam.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 6f8c61443..6f6f4efa6 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -1,4 +1,5 @@ -#include "Raster.h" +#define I_AM_MAIN +#include "MAPL_ErrLog.h" PROGRAM mkEASETilesParam @@ -21,6 +22,7 @@ PROGRAM mkEASETilesParam use process_hres_data use MAPL_SortMod use MAPL_ConstantsMod + use MAPL_ExceptionHandling use LogRectRasterizeMod use netcdf @@ -75,7 +77,8 @@ PROGRAM mkEASETilesParam character*1 :: PF character(len=6) :: EASE_Version character(len=10) :: nc_string, nr_string - character(128) :: usage1, usage2 + character(len=128) :: usage1, usage2 + character(len=128) :: Iam = "mkEASETilesParam" call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLandRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLandRaster.F90 similarity index 97% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLandRaster.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLandRaster.F90 index 1cf56391b..5c2494f28 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLandRaster.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLandRaster.F90 @@ -1,14 +1,14 @@ - -#include "Raster.h" +#define I_AM_MAIN +#include "MAPL_ErrLog.h" Program MakeLandRaster - + use MAPL_ExceptionHandling use LogRectRasterizeMod use MAPL_HashMod use process_hres_data use MAPL_SortMod use rmTinyCatchParaMod, ONLY: SRTM_maxcat, MAKE_BCS_INPUT_DIR - + use MAPL_Constants, only: PI=>MAPL_PI_R8 ! Program to create a surface raster file at 2.5' that has ! the ocean divided with a regular lat-lon DE grid. Its inputs ! are Sarith's formatted 2.5' raster of the Pfafstetter catchments with @@ -28,11 +28,11 @@ Program MakeLandRaster integer :: type, maxtiles, nx, ny integer :: count0,count1,count_rate - REAL_ :: dx, dy, d2r ! Grid spacing of raster grid - REAL_ :: xmin, ymin, xmax, ymax, xs, ys, da + real(kind=8) :: dx, dy, d2r ! Grid spacing of raster grid + real(kind=8) :: xmin, ymin, xmax, ymax, xs, ys, da - REAL_, allocatable :: cc(:), ss(:) - REAL_ , allocatable :: rTable(:,:) + real(kind=8), allocatable :: cc(:), ss(:) + real(kind=8) , allocatable :: rTable(:,:) integer, pointer :: Raster(:,:) integer, allocatable, target :: Raster0(:,:) @@ -45,8 +45,7 @@ Program MakeLandRaster logical :: Verb logical :: regrid, reynolds_sst = .false. - REAL_ :: VV(4) - REAL_ :: PI=RASTER_PI + real(kind=8) :: VV(4) ! ESA/SRTM ocean/land/ice/lake mask parameters ! -------------------------------------------- @@ -72,6 +71,7 @@ Program MakeLandRaster character*128 :: MaskFile character*128 :: & Usage = "mkLandRaster -x nx -y ny -v -h -z -t maxtiles -l LandFile -g GridName" + character*128 :: Iam = "MakeLandRaster" include 'netcdf.inc' call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) @@ -498,7 +498,7 @@ subroutine RegridRaster(Rin,Rout) integer, intent(IN) :: Rin(:,:) integer, intent(OUT) :: Rout(:,:) - REAL_ :: xx, yy + real(kind=8) :: xx, yy integer :: i,j,ii,jj xx = size(Rin ,1)/float(size(Rout,1)) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLatLonRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLatLonRaster.F90 similarity index 97% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLatLonRaster.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLatLonRaster.F90 index 6812a6bfe..05c220ebd 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkLatLonRaster.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLatLonRaster.F90 @@ -1,6 +1,5 @@ -! $Id: - -#include "Raster.h" +#define I_AM_MAIN +#include "MAPL_ErrLog.h" program mkLatLonRaster @@ -46,6 +45,7 @@ program mkLatLonRaster ! in each box. use LogRectRasterizeMod + use MAPL_ExceptionHandling implicit none @@ -64,6 +64,7 @@ program mkLatLonRaster real*8, allocatable :: xs(:), ys(:), xv(:,:,:), yv(:,:,:) character*128 :: & Usage = "mkLatLonRaster -x nx -y ny -v -h -z -g Gridname -b lon0 -p pos -t type im jm" + character*128 :: Iam = "mkLatLonRaster" ! Process Arguments !------------------ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMITAquaRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkMITAquaRaster.F90 similarity index 98% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMITAquaRaster.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkMITAquaRaster.F90 index e2a64a68c..efeb89d49 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMITAquaRaster.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkMITAquaRaster.F90 @@ -1,10 +1,10 @@ - -#include "Raster.h" -!#define VERIFY_(A) if (A/=0) stop 'verify error' +#define I_AM_MAIN +#include "MAPL_ErrLog.h" program MAIN use LogRectRasterizeMod + use MAPL_ExceptionHandling implicit none @@ -101,6 +101,7 @@ program MAIN integer, dimension(MAXBLNKSZ) :: blankList real(kind=RKIND) :: areamin, xc, yc + character(len=128) :: Iam = "mkMITAquaRaster" NAMELIST /W2_EXCH2_PARM01/ sNx, SNy, blankList diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMOMAquaRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkMOMAquaRaster.F90 similarity index 91% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMOMAquaRaster.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkMOMAquaRaster.F90 index 4f7dff4a1..6532c91d3 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkMOMAquaRaster.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkMOMAquaRaster.F90 @@ -1,11 +1,9 @@ -! $Id: - -#include "Raster.h" - +#define I_AM_MAIN +#include "MAPL_ErrLog.h" program MOMraster use LogRectRasterizeMod - + use MAPL_ExceptionHandling implicit none ! this program builds a rasterized grid whose cells are 2.5 by 2.5 minutes @@ -15,9 +13,9 @@ program MOMraster ! via namelist hence can be changed at runtime integer :: im, jm ! dimensions of MOM grid - REAL_, pointer :: xvert(:,:,:) ! Lons of MOM's vertices - REAL_, pointer :: yvert(:,:,:) ! Lats of MOM's vertices - REAL_ :: xmin, xmax + real(kind=8), pointer :: xvert(:,:,:) ! Lons of MOM's vertices + real(kind=8), pointer :: yvert(:,:,:) ! Lats of MOM's vertices + real(kind=8) :: xmin, xmax integer :: i, j, nxt,k integer :: status, command_argument_count character*(128) :: GridFile @@ -26,6 +24,7 @@ program MOMraster character*(2) :: opt character*(128) :: & Usage = "mkMOMAquaRaster -x rx -y ry -z -v -g GridName -h GridSpecFile" + character*(128) :: Iam = "mkMOMAquaRaster" ! argument defaults @@ -35,7 +34,7 @@ program MOMraster integer :: Nc = 8640 integer :: NR = 4320 - REAL_ :: tol + real(kind=8) :: tol INCLUDE "netcdf.inc" ! Process Arguments @@ -145,13 +144,13 @@ end subroutine FieldSize subroutine ReadGridFile(FILE,XVERT,YVERT) character*(*), intent(IN ) :: FILE - REAL_, pointer :: XVERT(:,:,:) - REAL_, pointer :: YVERT(:,:,:) + real(kind=8), pointer :: XVERT(:,:,:) + real(kind=8), pointer :: YVERT(:,:,:) integer :: STATUS, NCID, VARID integer :: SIZ_XVERT_X, SIZ_XVERT_Y integer :: SIZ_YVERT_X, SIZ_YVERT_Y - REAL_, pointer :: VERTX(:,:),VERTY(:,:) + real(kind=8), pointer :: VERTX(:,:),VERTY(:,:) Status=NF_OPEN(FILE,NF_NOWRITE,NCID) _ASSERT(STATUS==NF_NOERR,'needs informative message') diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mk_runofftbl.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mk_runofftbl.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mod_process_hres_data.F90 similarity index 99% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mod_process_hres_data.F90 index e2d978d1a..8efc7f826 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mod_process_hres_data.F90 @@ -22,8 +22,7 @@ MODULE process_hres_data use rmTinyCatchParaMod use MAPL_SortMod -use date_time_util -use leap_year +use LDAS_DateTimeMod use MAPL_ConstantsMod use lsm_routines, ONLY: sibalb use MAPL_Base, ONLY: MAPL_UNDEF diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rasterize.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rasterize.F90 similarity index 94% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rasterize.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rasterize.F90 index 8a3370b26..c742a1f24 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rasterize.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rasterize.F90 @@ -1,11 +1,11 @@ - -#include "Raster.h" +#include "MAPL_ErrLog.h" module LogRectRasterizeMod use MAPL_SORTMOD use rmTinyCatchParaMod, ONLY: SRTM_maxcat - + use MAPL_ExceptionHandling + use MAPL_Constants, only: PI=>MAPL_PI_R8 implicit none private @@ -29,14 +29,13 @@ module LogRectRasterizeMod public WriteLine integer, parameter :: PUSHLEFT = 10000 - REAL_ , parameter :: Zero = 0.0 - REAL_ , parameter :: PI = RASTER_PI + real(kind=8) , parameter :: Zero = 0.0 integer, parameter :: NX = 8640 integer, parameter :: NY = 4320 - REAL_ :: garea_ + real(kind=8) :: garea_ integer :: ctg_ interface LRRasterize @@ -113,7 +112,7 @@ end subroutine ReadRaster subroutine SortTiling(Raster,rTable,iTable) integer, intent(INOUT) :: Raster(:,:), iTable(0:,:) - REAL_, intent(INOUT) :: rTable(:,:) + real(kind=8), intent(INOUT) :: rTable(:,:) integer, dimension(size(iTable,2)) :: old, new integer*8, dimension(size(iTable,2)) :: key, key0 @@ -163,15 +162,16 @@ subroutine SortTiling(Raster,rTable,iTable) return end subroutine SortTiling -subroutine WriteTilingIR(File, GridName, im, jm, ipx, nx, ny, iTable, rTable, Zip, Verb) +subroutine WriteTilingIR(File, GridName, im, jm, ipx, nx, ny, iTable, rTable, Zip, Verb, rc) character*(*), intent(IN) :: File character*(*), intent(IN) :: GridName(:) integer, intent(IN) :: nx,ny integer, intent(IN) :: iTable(0:,:) - REAL_, intent(IN) :: rTable(:,:) + real(kind=8), intent(IN) :: rTable(:,:) integer, intent(IN) :: IM(:), JM(:), ipx(:) logical, optional, intent(IN) :: Zip logical, optional, intent(IN) :: Verb + integer, optional, intent(out) :: rc ! Table variables ! @@ -192,10 +192,11 @@ subroutine WriteTilingIR(File, GridName, im, jm, ipx, nx, ny, iTable, rTable, Zi integer :: j, unit, ng, ip, l, i, k, ix character*1000 :: Line integer :: ii(size(GridName)), jj(size(GridName)), kk(size(GridName)) - REAL_ :: fr(size(GridName)) - REAL_ :: xc, yc, area - REAL_ :: garea, ctg(size(Gridname)) - REAL_ :: sphere, error + real(kind=8) :: fr(size(GridName)) + real(kind=8) :: xc, yc, area + real(kind=8) :: garea, ctg(size(Gridname)) + real(kind=8) :: sphere, error + integer :: status ip = size(iTable,2) ng = size(GridName) @@ -388,7 +389,7 @@ subroutine WriteLine(File, Unit, iTable, rTable, k, Zip, Verb) character*(*), intent(IN) :: File integer, intent(IN) :: Unit, k integer, intent(IN) :: iTable(0:) - REAL_, intent(IN) :: rTable(:) + real(kind=8), intent(IN) :: rTable(:) logical, optional, intent(IN) :: Zip logical, optional, intent(IN) :: Verb @@ -410,8 +411,8 @@ subroutine WriteLine(File, Unit, iTable, rTable, k, Zip, Verb) logical :: DoZip character*1000 :: Line integer :: ii, jj - REAL_ :: fr - REAL_ :: xc, yc, area + real(kind=8) :: fr + real(kind=8) :: xc, yc, area if(present(Zip)) then DoZip = Zip @@ -477,7 +478,7 @@ subroutine CloseTiling(FIle, Unit, ip, Zip, Verb) ! rTable(5) :: of 2nd grid box area logical :: DoZip - REAL_ :: sphere, error + real(kind=8) :: sphere, error character*1000 :: Line Line="" diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rasterize.H b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rasterize.H similarity index 87% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rasterize.H rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rasterize.H index e8c6bbde0..f2cf014ab 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rasterize.H +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rasterize.H @@ -1,6 +1,4 @@ - - #ifdef MESH #define LRRasterize2Mem LRRasterize2Mem0 #define LRRasterize2File LRRasterize2File0 @@ -8,12 +6,13 @@ #else #define POLYSIZE size(xv,3) #endif +#define RASTERUNDEF -999 !BOP ! !IROUTINE: LRRasterize -- Rasterizes a 2-D array of polygons. subroutine LRRasterize2File(GridName, xv,yv,nc,nr,xmn,xmx,ymn,ymx, & - SurfaceType, Verb, Here,jseg, tol ) + SurfaceType, Verb, Here,jseg, tol, rc ) ! This routine rasterizes a grid defined by a 2-D array of polygons. ! The raster value assigned each pixel is either sign((I*pushleft + J),sgn), @@ -55,22 +54,23 @@ character*(*), intent(INOUT) :: GridName ! Raster file name #ifdef MESH - REAL_, intent(INOUT) :: xv(:,:) ! X coordinates of vertices - REAL_, intent(INOUT) :: yv(:,:) ! Y coordinates of vertices + real(kind=8), intent(INOUT) :: xv(:,:) ! X coordinates of vertices + real(kind=8), intent(INOUT) :: yv(:,:) ! Y coordinates of vertices #else - REAL_, intent(INOUT) :: xv(:,:,:) ! X coordinates of vertices - REAL_, intent(INOUT) :: yv(:,:,:) ! Y coordinates of vertices + real(kind=8), intent(INOUT) :: xv(:,:,:) ! X coordinates of vertices + real(kind=8), intent(INOUT) :: yv(:,:,:) ! Y coordinates of vertices #endif integer, optional, intent(IN) :: nc,nr ! Raster field sizes - REAL_, optional, intent(IN) :: xmn ! LL x of LL raster cell (-180) - REAL_, optional, intent(IN) :: ymn ! LL y of LL raster cell ( -90) - REAL_, optional, intent(IN) :: xmx ! UR x of UR raster cell ( 180) - REAL_, optional, intent(IN) :: ymx ! UR y of UR raster cell ( 90) + real(kind=8), optional, intent(IN) :: xmn ! LL x of LL raster cell (-180) + real(kind=8), optional, intent(IN) :: ymn ! LL y of LL raster cell ( -90) + real(kind=8), optional, intent(IN) :: xmx ! UR x of UR raster cell ( 180) + real(kind=8), optional, intent(IN) :: ymx ! UR y of UR raster cell ( 90) logical, optional, intent(IN) :: verb ! Verbose logical, optional, intent(IN) :: here ! write here integer, optional, intent(IN) :: SurfaceType integer, optional, intent(IN) :: jseg - REAL_, optional :: tol + real(kind=8), optional :: tol + integer, optional, intent(out) :: rc character*(128) :: TileFile character*(128) :: TilDir='', RstDir='' @@ -172,7 +172,7 @@ subroutine LRRasterize2Mem(Raster,xv,yv,Tilefile, & xmn,xmx,ymn,ymx, & - SurfaceType,verb,jseg, tol ) + SurfaceType,verb,jseg, tol, rc ) ! This routine rasterizes a grid defined by a 2-D array of polygons. ! The raster value assigned is an index to the table. @@ -191,41 +191,42 @@ integer, intent(INOUT) :: Raster(:,:) ! Raster field to be filled #ifdef MESH - REAL_, intent(INOUT) :: xv(:,: ) ! X coordinates of vertices - REAL_, intent(INOUT) :: yv(:,: ) ! Y coordinates of vertices + real(kind=8), intent(INOUT) :: xv(:,: ) ! X coordinates of vertices + real(kind=8), intent(INOUT) :: yv(:,: ) ! Y coordinates of vertices #else - REAL_, intent(INOUT) :: xv(:,:,:) ! X coordinates of vertices - REAL_, intent(INOUT) :: yv(:,:,:) ! Y coordinates of vertices + real(kind=8), intent(INOUT) :: xv(:,:,:) ! X coordinates of vertices + real(kind=8), intent(INOUT) :: yv(:,:,:) ! Y coordinates of vertices #endif character*(*), intent(IN ) :: TileFile - REAL_, optional, intent(IN) :: xmn ! LL x of LL raster cell (-180) - REAL_, optional, intent(IN) :: ymn ! LL y of LL raster cell ( -90) - REAL_, optional, intent(IN) :: xmx ! UR x of UR raster cell ( 180) - REAL_, optional, intent(IN) :: ymx ! UR y of UR raster cell ( 90) + real(kind=8), optional, intent(IN) :: xmn ! LL x of LL raster cell (-180) + real(kind=8), optional, intent(IN) :: ymn ! LL y of LL raster cell ( -90) + real(kind=8), optional, intent(IN) :: xmx ! UR x of UR raster cell ( 180) + real(kind=8), optional, intent(IN) :: ymx ! UR y of UR raster cell ( 90) logical, optional, intent(IN) :: verb ! Verbose integer, optional, intent(IN) :: SurfaceType integer, optional, intent(IN) :: jseg - REAL_, optional :: tol + real(kind=8), optional :: tol + integer, optional, intent(out) :: rc ! X abd Y bounds of each polygon - REAL_ :: xmin, xmax - REAL_ :: ymin, ymax - REAL_ :: minx, miny - REAL_ :: maxx, maxy + real(kind=8) :: xmin, xmax + real(kind=8) :: ymin, ymax + real(kind=8) :: minx, miny + real(kind=8) :: maxx, maxy ! x and y coordinates of the Raster grid - REAL_, dimension(size(Raster,1)) :: xs, xcs, xss - REAL_, dimension(size(Raster,2)) :: ys, ycs, yss, da + real(kind=8), dimension(size(Raster,1)) :: xs, xcs, xss + real(kind=8), dimension(size(Raster,2)) :: ys, ycs, yss, da integer :: IM, JM, NV ! Shape of input grid - REAL_ :: dx, dy, dxi, dyi ! Grid spacing of raster grid + real(kind=8) :: dx, dy, dxi, dyi ! Grid spacing of raster grid integer :: xsize, ysize ! Dimensions of Raster grid integer :: i, j, jn, n, ib, jb, fill, uType, js, k - REAL_ :: range, d2r, r2d, ddx, grid_ymin, grid_ymax, xc, yc, Area, xx + real(kind=8) :: range, d2r, r2d, ddx, grid_ymin, grid_ymax, xc, yc, Area, xx logical :: DoZip, uVerb integer :: idx, ct integer :: count0,count1,count_rate @@ -234,18 +235,18 @@ character*(128) :: GridName, TilFile integer, allocatable :: iTable(:,:) - REAL_, allocatable :: rTable(:,:) + real(kind=8), allocatable :: rTable(:,:) integer :: useg, unit, fq integer, dimension(POLYSIZE) & :: nxt - REAL_, dimension(POLYSIZE) & + real(kind=8), dimension(POLYSIZE) & :: xvc, yvc, xvs, yvs, xrd, yrd, x3, y3, z3, & dx3, dy3, dz3, x31, x32, y31, y32, z31, z32, & dx4, dy4, x4, y4, xu, yu - REAL_ :: tol_ + real(kind=8) :: tol_ ! Process optionals @@ -343,8 +344,8 @@ dxi = 1.0/dx dyi = 1.0/dy - d2r = (2._8*RASTER_PI)/range - r2d = range/(2._8*RASTER_PI) + d2r = (2._8*PI)/range + r2d = range/(2._8*PI) ! Report @@ -557,14 +558,14 @@ subroutine FillPoly(sh) - REAL_, intent(IN) :: sh + real(kind=8), intent(IN) :: sh logical :: IsIn integer :: i1, i2, jj1 integer :: ii, jj, n1, n2, jx integer, save :: j1, j2 integer :: HALO=10 - REAL_ :: x0, y0 + real(kind=8) :: x0, y0 if (abs(miny+90._8) < 1.e-10_8) then diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/read_riveroutlet.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/read_riveroutlet.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/read_riveroutlet.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/read_riveroutlet.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 similarity index 99% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 index 61a8f55e0..9171524e6 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 @@ -7,8 +7,7 @@ module rmTinyCatchParaMod - use date_time_util - use leap_year + use LDAS_DateTimeMod use MAPL_ConstantsMod use MAPL_Base, ONLY: MAPL_UNDEF use lsm_routines, ONLY: sibalb diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/util.c b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/util.c similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/util.c rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/util.c diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/zip.c b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/zip.c similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/zip.c rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/zip.c diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/README b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/README new file mode 100644 index 000000000..c1af76617 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/README @@ -0,0 +1,2 @@ +Directory [..]/GEOSsurface_GridComp/Utils/Raster/misc contains miscellaneous debugging and plotting scripts. + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/chk_clsm_params.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/chk_clsm_params.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/chk_clsm_params.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/chk_clsm_params.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/compare_bcs.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/compare_bcs.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/compare_bcs.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/compare_bcs.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/get_frac.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/get_frac.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/get_frac.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/get_frac.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/modis_scale_factor.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/modis_scale_factor.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/modis_scale_factor.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/modis_scale_factor.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mosaic_classes_on_tiles.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/mosaic_classes_on_tiles.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mosaic_classes_on_tiles.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/mosaic_classes_on_tiles.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/plot_curves.csh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/plot_curves.csh similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/plot_curves.csh rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/plot_curves.csh diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/plot_curves.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/plot_curves.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/plot_curves.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/plot_curves.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/plot_geos5_grid.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/plot_geos5_grid.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/plot_geos5_grid.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/plot_geos5_grid.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/soil_types_on_tiles.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/soil_types_on_tiles.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/soil_types_on_tiles.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/misc/soil_types_on_tiles.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/ConvertAlb.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/ConvertAlb.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/ConvertAlb.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/ConvertAlb.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/asia_tiles.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/asia_tiles.pro similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/asia_tiles.pro rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/asia_tiles.pro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/compile_mk_runoff b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/compile_mk_runoff similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/compile_mk_runoff rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/compile_mk_runoff diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_tiles b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/make_tiles similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_tiles rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/make_tiles diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyTiles.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/rmTinyTiles.F90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyTiles.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/obsolete/rmTinyTiles.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt new file mode 100644 index 000000000..15405d696 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt @@ -0,0 +1 @@ +esma_add_subdirectories (soil) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/README b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/README new file mode 100644 index 000000000..efe175265 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/README @@ -0,0 +1 @@ +Directory [..]/GEOSsurface_GridComp/Utils/Raster/preproc/ contains programs and scripts needed to generate *input* files for make_bcs. diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/get_lat_lon4tils.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/get_lat_lon4tils.pro new file mode 100755 index 000000000..bdb3558b3 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/get_lat_lon4tils.pro @@ -0,0 +1,81 @@ +pro get_lat_lon4tils + +; code to derive lat/lons for each MOD10A1 grid point. The code +; loops over MODIS files (hereinafter refered to as tiles), its +; grid boxes and calculates lat/lon values using the same approach +; as given in C function provided by MODIS team (publically available). +; The lat/lon values for each gridbox are written out in 10-degx10-deg files +; (to match MODIS files; each 2400x2400 elements), as well as in two +; global files holding lat and lon each. The code will print +; out even the tiles (and gridboxes) with no valid data (miss val:-999.) + +; to run in terminal window: +; type 'idl', then 'r. get_lat_lon4tils', then 'get_lat_lon4tils' +; requirements: a subdirectory 'MODIS_lat_lon' in the local space +; created: March 2022 Biljana Orescanin SSAI@NASA + +mis_val=-999. + +lat_stiched=make_array(2400l*36, 2400l*18, value=mis_val,/float) +lon_stiched=make_array(2400l*36, 2400l*18, value=mis_val,/float) ; Array[864, 432] + +; Set Map projection +sinusMap = MAP_PROJ_INIT('Sinusoidal', DATUM=8, CENTER_LAT=0., CENTER_LON=0,SPHERE_RADIUS=6371007.181) + +; loop over MOIDS tiles +for ivtil=0,18-1 do begin + for ihtil=0,36-1 do begin + + ; get 2-digit strings for vert and horiz tile/file numbering + v_str=strmid('0'+strtrim(ivtil,2),1,2,/reverse) + h_str=strmid('0'+strtrim(ihtil,2),1,2,/reverse) + + ; declare 2D arays to store lat/lons for the current file + lat2d=make_array(2400l, 2400l, value=mis_val,/float) ; fltarr(2400,2400) + lon2d=make_array(2400l, 2400l, value=mis_val,/float) ; + + ; loop over grid boxes within the current file + for ivgrid=0,2400-1 do begin + for ihgrid=0,2400-1 do begin + + ; get lats and lons from UV coordinates (requires the map projection loaded above) + xOrigin=-20015109.354d +((ihtil)*2400l+ihgrid)*463.31271653d + yOrigin= 10007554.677d -((ivtil)*2400l+ivgrid)*463.31271653d + lonlat = MAP_PROJ_INVERSE([xOrigin], [yOrigin], MAP_STRUCTURE=sinusMap) + lat =lonlat[1] + lon =lonlat[0] + + if lat eq lat then begin + lat2d(ihgrid,ivgrid)=lat + lon2d(ihgrid,ivgrid)=lon + i_ind=(ihtil)*2400l+ihgrid + j_ind=(ivtil)*2400l+ivgrid + lat_stiched(i_ind,j_ind)=lat + lon_stiched(i_ind,j_ind)=lon + endif + + endfor ; ihgrid + endfor ; ivgrid + + ;write 2D lat/lons into binary file + openw, 10, 'MODIS_lat_lon/MODIS_hdres_lat_lon_v'+v_str+'_h'+h_str+'.bin.gz', /compress + writeu,10, lat2D,lon2D + close, 10 + + ; place the current tile in global gird + lat_stiched((ihtil)*2400l:(ihtil)*2400l+2399,(ivtil)*2400l:(ivtil)*2400l+2399l)=lat2d + lon_stiched((ihtil)*2400l:(ihtil)*2400l+2399,(ivtil)*2400l:(ivtil)*2400l+2399l)=lon2d + + endfor ; ihtil +endfor ; ivril + +; write stiched lat/lons into a text file +openw , 20, 'MODIS_lat_lon/MODIS_lat_stiched_hdres_global.txt.gz', /compress +printf, 20, lat_stiched +close , 20 + +openw , 30, 'MODIS_lat_lon/MODIS_lon_stiched_hdres_global.txt.gz', /compress +printf, 30, lon_stiched +close , 30 + +end diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/get_snow_alb_mod10a1_30arcsec.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/get_snow_alb_mod10a1_30arcsec.pro new file mode 100755 index 000000000..487fe8c4d --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/get_snow_alb_mod10a1_30arcsec.pro @@ -0,0 +1,569 @@ +pro get_snow_alb_mod10a1_30arcsec, year=year, h_s=h_s, h_e=h_e + +; Code to stitch MODIS files, such as: MOD10A1.A2022002.h35v08.006.2022004050029.hdf +; into a static global 30arcsec grid and produce inputs for GEOS make_bcs package, such as: +; snow_alb_MOD10A1.061_30arcsec_H36V05.nc +; The original MODIS input files are on a sinusoidal (~500-m) grid. The final snow albedo +; product is on a 30arcsec regular lat/lon grid (date-line-edge and pole-edge) + +; It first reads the MODIS-tile files to generate snow cover PDFs over a period of time. +; These PDFs are made based on grid cells that exceed a snow cover value (e.g. 60%). This is done +; so a decent amount of snow over a grid box ensures a better estimate of albedo. +; The PDFs are then used to locate a user-defined cutoff for a %ile of snow coverage +; (e.g. top 10 %tile) and filter snow albedo values. + +; Input argument 'year' in the call specifies for what year(s) to run the code (to allow multiple +; serial runs). Only the PDF building can be processed using multiple serial runs, the remainder +; of the code has to execute on a single CPU. If no year(s) is provided on the input, the default +; year is used. Input argument can be an integer or an array of integers indicating year(s) to +; process (e.g. year=[2019,2020,2021]) + +; Input arguments 'h_s' and 'h_e' are starting and ending horizontal tiles. These are optional +; and if provided define the range of horizontal MODIS tiles to be processed (it saves time if +; certain tiles need re-processing; otherwise should be ignored) + +; dependencies: +; grid.pro (performes gridding) +; read_hdf_sd.pro (reads hdf files) +; read_mod10a1_hdf.pro (reads MODIS files) +; get_lat_lon4tils.pro (perform conversion of original MODIS sinusoidal projection into lat/lon space) +; +; [sub]directories to be created prior executing this code: +; data/ - output directory (user needs to create) +; MODIS_lat_lon/ - holds previousely created lat/lon info (user needs to create this with get_lat_lon4tils.pro listed in dependencies) +; MOD10A1_data/ - holds MODIS input data (user has to bring it in; follow the dir structure) + +; There are five steps (runs) to process. Steps have to be processed in order. The current step must +; complete before the next is initiated. Use the 'goto' commands to control which step is exectued. +; When the current step is completed, uncomment the next-step goto command on lines 86-92 + +; To execture the code in terminal window type (without quotations): +; 'idl', then '.r get_snow_alb_mod10a1_30arcsec', then 'get_snow_alb_mod10a1_30arcsec' + +; Created April 2022 Biljana Orescanin SSIA@NASA + +; path to MODIS data (file name example: MOD10A1.A2022001.h24v06.061.2022003035827.hdf) + +path_in='./MOD10A1_data/' ; here is only one file to provide an example + +print, "You must make sure the 'path_in' is of the expected structure" + +; set parameters to loop over +if keyword_set(year) then begin + year=year[sort(year)] ; set an increasing order + years=strtrim(year,2) + print, ' Processing year(s): ', years +endif else begin + print,'No year to run specified on the input. Processing for default year - 2020!' + years =['2020'] +endelse + +snow_cvr_min = 0. ; minimum snow cover to consider [%] +cvr_cutoff = 0 ; snow cover cutoff in %s (to consider only tiles gt this value) +top_alb_limit= 99 ; the top %ile for snow albedo (I choose this value) to be considered (e.g. top 10% of the PDF) +top_limit = 10 ; the top %ile for snow cover (I choose this value) to be considered (e.g. top 10% of the PDF) +mis_val =-999. ; has to be in range [0,255] b/c snow_alb_til is byte array +snow_min_str=strtrim(snow_cvr_min,2) +year_start =years(0) +year_end =years(-1) +print, 'Start/End Years:',year_start,'/',year_end + +if keyword_set(h_s) and keyword_set(h_e) then begin + h_start=h_s + h_end =h_e +endif else begin + h_start=0 + h_end =35 +endelse +print, 'h_start:',h_start,' h_end:',h_end + +v_start=0 +v_end =17 + +dim_1d=2400l*2400l +dim_2d=[2400l,2400l] + +; ********* STEPS *********** +; [un]comment out the following goto's depending on what part of the code is to be executed +;goto, skip_extracting_snow_cvr ; *** STEP 2 +;goto, skip_calculating_snow_cvr ; *** STEP 3 +;goto, skip_reading_snow_alb ; *** STEP 4 +;goto, skip2ncoutput ; *** STEP 5 +; ********* STEPS *********** + +; loop over years +for iyear=0,n_elements(years)-1 do begin + + ; loop over horiontal tiles + for ih=h_start,h_end do begin + ih_str=strmid('0'+strtrim(ih,2),1,2,/reverse) + + ; loop over vertical tiles + for iv=v_start,v_end do begin + iv_str=strmid('0'+strtrim(iv,2),1,2,/reverse) + + ; get all files for the current h/v tile + filename=file_search(path_in+years(iyear)+ $ + '/n5eil01u.ecs.nsidc.org/DP4/MOST/MOD10A1.061/'+ $ + years(iyear)+'.'+'*.*/MOD10A1.A'+years(iyear)+ $ + '*.h'+ih_str+'v'+iv_str+'.*.hdf',count=nfiles) + + ; if no files available for this tile, skip to the next one (note: ~ 1/3 of tiles will be empty) + if nfiles eq 0 then begin + print, 'No data for tile: h'+ih_str+'v'+iv_str, ' skipping' + continue + endif + + ; declare arrays to hold all days info + snow_cvr_pdf=make_array(dim_1d,101,/long,value=0l) ; to hold pdf of snow cover values + + ; loop over files (i.e. days) for this tile to stitch + for ifile=0,nfiles-1 do begin + + ; read file + data=read_mod10a1_hdf(filename(ifile)) + + ; keep only what you need + ; See bottom of the code for details + snow_cvr_til=data.NDSI_Snow_Cover ; values [0,100] + snow_cvr_qc =data.NDSI_Snow_Cover_Basic_QA ; values [0,255] 0=best, 1=good, 2=ok, 3=poor-not used + + ; reform all arrays into 1D arrays so I can apply indexes from where command + snow_cvr_til=reform(snow_cvr_til,dim_2d[0]*dim_2d[1],/overwrite) + snow_cvr_qc =reform(snow_cvr_qc ,dim_2d[0]*dim_2d[1],/overwrite) + + ; If data present then accumulate non-missing values and their counts + ind_val=where(snow_cvr_til gt snow_cvr_min and snow_cvr_til le 100. and $ + snow_cvr_qc eq 0,c_val) + + ; If no valida data of albedo and snow cover in this tile, skip. + if c_val eq 0 then begin + continue + endif + + ; Get PDFs of snow fraction and albedo + for ival=0,c_val-1 do begin + snow_cvr_pdf(ind_val(ival),fix(snow_cvr_til[ind_val(ival)]))++ + endfor + + endfor ; ifile + + openw, lun, 'data/data_out/snow_cov_pdfs_titch08_OD10A1.A.'+years(iyear)+ $ + '.h'+ih_str+'v'+iv_str+'.bin.gz', /get_lun,/compress + writeu,lun, snow_cvr_pdf + free_lun,lun + + endfor ; iv + endfor ; ih +endfor ; iyear + +print, 'Done creating PDFs of snow cover and albedo' +stop + +skip_extracting_snow_cvr: +print, 'Starting snow fraction top percentile calculations' + + ; -- get snow fractions for the top n-percentile of PDFs at each pixel + ; This will take a while since it has to be done pixel by pixel + + ; Few notes: + ; - Needed is a mean albedo corresponding to the top 10% (by count) snow fraction values + ; but only among those vaues that correspond to 60% or more of snow over a grid box. + ; - So, used is a snow cover cutoff value of 60 to get the top 90the percentile of snow cover values + +; loop over horiontal tiles +for ih=h_start,h_end do begin + ih_str=strmid('0'+strtrim(ih,2),1,2,/reverse) + + ; loop over vertical tiles + for iv=v_start,v_end do begin + iv_str=strmid('0'+strtrim(iv,2),1,2,/reverse) + + snow_cvr_pdf_tmp=make_array(dim_1d,101,/long,value=0l) ; to hold current pdf of snow cover values + snow_cvr_pdf_all=make_array(dim_1d,101,/long,value=0l) ; to hold cummulatinve pdfs of snow cover values + + ; First read PDFs of all the years for each tile. Then, find the top %ile and write out the cutoffs + + ; get all files for the current h/v tile (all available years) + filename_cvr=file_search('data/data_out/snow_cov_pdfs_titch08_OD10A1.A.*.h'+ih_str+'v'+iv_str+'.bin.gz',count=nfiles2) + + ; if no files available for this tile, skip to the next one (note: ~ 1/3 of tiles will be empty) + if nfiles2 eq 0 then begin + print, 'No data for tile: h'+ih_str+'v'+iv_str, ' skipping' + continue + endif + + for ifile=0,nfiles2-1 do begin + openr, lun, filename_cvr(ifile), /get_lun,/compress + readu, lun, snow_cvr_pdf_tmp + free_lun,lun + + snow_cvr_pdf_all=snow_cvr_pdf_all+snow_cvr_pdf_tmp + endfor + + ; find the snow cover cutoff + sf_lim_tail_mean=make_array(dim_1d,/float,value=mis_val) ; array to store the top 90th snow fraction %ile cutoff for each pixel + + for ipix=0,dim_1d-1 do begin + + tot_pix=total(snow_cvr_pdf_all(ipix,*)) ; tot # of valid snow fraction obs at this pixel + + if tot_pix lt 10 then continue ; if lt 10 valid elements leave as missing (no enough data) + + tot_sf_pix=0l + for ibin=100,30,-1 do begin ; accumulate from the top bin down till you get enough data. + ; Yet, don't go below 30th bin (i.e. snow fraction lt 30%) + ; b/c there is not enough snow to make it a "reliable albedo estimate" + tot_sf_pix=tot_sf_pix+snow_cvr_pdf_all(ipix,ibin)*ibin + if total(snow_cvr_pdf_all(ipix,ibin:100))/tot_pix*100. gt top_limit then begin + sf_lim_tail_mean(ipix)=tot_sf_pix/total(snow_cvr_pdf_all(ipix,ibin:100)) + break ; ibin loop + endif + endfor ; ibin + + endfor ; ipix + + ; write out the snow cover cutoff %ile for this MODIS tile + openw, lun, 'data/data_out/snow_cov_cutoff_titch08_MOD10A1.A.h'+ih_str+'v'+iv_str+'_'+year_start+'_'+year_end+'.bin.gz', /get_lun,/compress + writeu,lun, sf_lim_tail_mean ;,sf_max ; both are [dim_1d] float arrays + free_lun,lun + + endfor ; iv +endfor ;ih + +print, 'Done snow fraction top percentile calculations' +stop + +skip_calculating_snow_cvr: +print, 'Starting reading in the albedo data to get mean value for those grids that have above the limit snow cover' + +; loop over years +for iyear=0,n_elements(years)-1 do begin + + ; loop over horiontal tiles + for ih=h_start,h_end do begin + ih_str=strmid('0'+strtrim(ih,2),1,2,/reverse) + + ; loop over vertical tiles + for iv=v_start,v_end do begin + iv_str=strmid('0'+strtrim(iv,2),1,2,/reverse) + + ; get all files for the current h/v tile + filename=file_search(path_in+years(iyear)+'/n5eil01u.ecs.nsidc.org/DP4/MOST/MOD10A1.061/'+years(iyear)+'.'+ $ + '*.*/MOD10A1.A'+years(iyear)+'*.h'+ih_str+'v'+iv_str+'.*.hdf',count=nfiles3) + + ; if no files available for this tile, skip to the next one (note: ~ 1/3 of tiles will be empty) + if nfiles3 eq 0 then begin + print, 'No data for tile: h'+ih_str+'v'+iv_str, ' skipping' + continue + endif + + ; read in the snow cover cutoff %ile for this MODIS tile + sf_lim_tail_mean=make_array(dim_1d,/float,value=mis_val) ; to store the top snow fraction %ile cutoff for each pixel + openr, lun, 'data/data_out/snow_cov_cutoff_titch08_MOD10A1.A.h'+ih_str+'v'+iv_str+'_'+year_start+'_'+year_end+'.bin.gz', /get_lun,/compress + readu,lun, sf_lim_tail_mean + free_lun,lun + + ; declare arrays to hold all days info + accu_alb=make_array(dim_1d,/float,value=0l) ; to hold accumulated snow albedo values + accu_cnt=make_array(dim_1d,/float,value=0l) ; to hold couhnts of accumulated values + + ; loop over files (i.e. days) for this tile to stitch + for ifile=0,nfiles3-1 do begin + + ; read file + data=read_mod10a1_hdf(filename(ifile)) + + ; keep only what you need ; see bottom of the code for details + snow_alb_til=data.SNOW_ALBEDO_DAILY_TILE ; values [0,100] + snow_cvr_til=data.NDSI_Snow_Cover ; values [0,100] + snow_cvr_qc =data.NDSI_Snow_Cover_Basic_QA ; values [0,255] 0=best, 1=good, 2=ok, 3=poor-not used + + ; reform all arrays into 1D arrays so I can apply indexes from where command + snow_alb_til=reform(snow_alb_til,dim_2d[0]*dim_2d[1],/overwrite) + snow_cvr_til=reform(snow_cvr_til,dim_2d[0]*dim_2d[1],/overwrite) + snow_cvr_qc =reform(snow_cvr_qc ,dim_2d[0]*dim_2d[1],/overwrite) + + ; if no valid data of abledo and snow cover in this tile, skip. + ; If data present then accumulate non-missing values and their counts + ind_val=where(snow_alb_til gt 0 and snow_alb_til le 100 and $ + snow_cvr_til gt 0 and snow_cvr_til le 100 and $ + sf_lim_tail_mean gt 0 and snow_cvr_qc eq 0,c_val) + + if c_val eq 0 then begin + continue + endif + + ; get accumulated albedo + for ipix=0,c_val-1 do begin + if snow_cvr_til(ind_val[ipix]) ge sf_lim_tail_mean(ind_val[ipix]) then begin + accu_alb(ind_val[ipix])=accu_alb(ind_val[ipix])+snow_alb_til(ind_val[ipix]) + accu_cnt(ind_val[ipix])=accu_cnt(ind_val[ipix])+1. + endif + endfor ; ipix + + endfor ; ifile + + ; write out the counts, cumulative and mean values of albedo and PDFs of albedo + ; for this tile (these will be stitched once all completed) + openw, lun, 'data/data_out/snow_alb_pdfs_08_MOD10A1.A.'+years(iyear)+ $ + '.h'+ih_str+'v'+iv_str+'_'+strtrim(100-top_limit,2)+'%ile_cover2_gt_'+ $ + strtrim(cvr_cutoff,2)+'.bin.gz', /get_lun,/compress + writeu,lun, accu_alb,accu_cnt + free_lun,lun + + endfor ; iv + endfor ;ih +endfor ; iyear +print, 'Done with calculating mean snow albedo' +stop + +skip_reading_snow_alb: + +; -- All snow albedo PDFs are in, only in yearly files. Read them all +; and make the cumulative stats for mean albedo at a chosen %ile +; This is to be done on a single CPU as a single sbatch job +print, 'Reading in all Snow Albedo accum and counts to form the overall mean albedo values' + +; need to read in all the years for each tile, find top %ile and write out the mean albedo +; loop over horizontal tiles +for ih=h_start,h_end do begin + ih_str=strmid('0'+strtrim(ih,2),1,2,/reverse) + + ; loop over vertical tiles + for iv=v_start,v_end do begin + iv_str=strmid('0'+strtrim(iv,2),1,2,/reverse) + + ; get all files for the current h/v tile (all available years) + filename_alb=file_search('data/data_out/snow_alb_pdfs_08_MOD10A1.A.*.h'+ $ + ih_str+'v'+iv_str+'_'+strtrim(100-top_limit,2) + $ + '%ile_cover2_gt_'+strtrim(cvr_cutoff,2)+'.bin.gz',count=nfiles4) + + ; if no files available for this tile, skip to the next one (note: ~ 1/3 of tiles will be empty) + if nfiles4 eq 0 then begin + print, 'No data for tile: h'+ih_str+'v'+iv_str, ' skipping' + continue + endif + + print, 'tile: h'+ih_str+'v'+iv_str + + ; declare arrays to hold all days info + accu_alb_tmp=make_array(dim_1d,/float,value=0l) ; to hold accumulated snow albedo values + accu_cnt_tmp=make_array(dim_1d,/float,value=0l) ; to hold counts of accumulated values + snow_alb_accu_all=0. + snow_alb_cnt_all =0. + + for ifile=0,nfiles4-1 do begin + openr, lun, filename_alb(ifile), /get_lun,/compress + readu, lun, accu_alb_tmp,accu_cnt_tmp + free_lun,lun + snow_alb_accu_all=snow_alb_accu_all+accu_alb_tmp + snow_alb_cnt_all =snow_alb_cnt_all +accu_cnt_tmp + endfor + + ; get mean albedo + mean_alb=snow_alb_accu_all/snow_alb_cnt_all + mean_alb[where(snow_alb_cnt_all lt 10,/null)]=mis_val + + ; write out the snow albedo max and cutoff %ile for this MODIS tile + openw, lun, 'data/data_out/snow_alb_08_'+strtrim(top_alb_limit,2)+'_cutoff_MOD10A1.A.h'+ $ + ih_str+'v'+iv_str+'_'+year_start+'_'+year_end+'.bin.gz', /get_lun,/compress + writeu,lun, mean_alb + free_lun,lun + + ; reset arrays for accumulating snow albedo + snow_alb_accu_all(*)=0l + snow_alb_cnt_all (*)=0l + + endfor ; iv +endfor ;ih + +print, 'Done with snow albedo cutoff %ile calculations' +stop + +skip2ncoutput: +; -- Now all the albedo vaules are in. Output them on +; GEOS-friendly 10x10 deg tiles using nc format + +; set new missing value +mis_val= 1.e15 + +print, 'Starting stitching' + +; array to store the top albedo %ile cutoff for each pixel at a given MODIS tile +alb_lim_tail_mean=make_array(dim_1d,/float,value=mis_val) + +; create 36x18 tiles at 30sec arc resolution (each tile is 10x10deg with 1200x1200 grid boxes) +; loop over horiontal tiles +for ih=h_start,h_end do begin + ih_str=strmid('0'+strtrim(ih+1,2),1,2,/reverse) + + ; loop over vertical tiles + for iv=v_start,v_end do begin + iv_str=strmid('0'+strtrim(iv+1,2),1,2,/reverse) + + ; declare array to store the output for this tile; fill with missing + alb_30sec_grid=make_array(1200l, 1200l, value=mis_val,/float) + + print, 'Creating tile: h'+ih_str+'v'+iv_str + + ; get min/max and all lat/lon values for this tile + minlat = iv *10.-90. + maxlat =(iv+1)*10.-90. + minlon = ih *10.-180.0 + maxlon =(ih+1)*10.-180.0 + lat_positions=minlat+indgen(1200)*10./1200. + lon_positions=minlon+indgen(1200)*10./1200. + + ; have to read +/- 1 MODIS tile in vertical direction + ; and as many as needed to cover all valid lat/lons in horizontal direction + + alb_lim_tail_mean_all=[] ; To store all 500m resolution albedo values + lat_all =[] ; corresponding to this tile to be used in + lon_all =[] ; gridding + + ; loop over vertical tiles + for iiv=17-iv-1,17-iv+1 do begin + iiv_str=strmid('0'+strtrim(iiv,2),1,2,/reverse) + + ; loop over horizontal tiles + for iih=0,35 do begin + iih_str=strmid('0'+strtrim(iih,2),1,2,/reverse) + + ; read in the cumulative and mean values for this tile (these will be stitched once all completed) + filename_alb=file_search('data/data_out/snow_alb_08_'+strtrim(top_alb_limit,2)+ $ + '_cutoff_MOD10A1.A.h'+iih_str+'v'+iiv_str+'_'+year_start+ $ + '_'+year_end+'.bin.gz',count=n_files5) + + if n_files5 ne 1 then continue + + openr, lun, filename_alb, /get_lun,/compress + readu,lun, alb_lim_tail_mean + free_lun,lun + + ; Reform the arrays back to 2 dimensions + alb_lim_tail_mean=reform(alb_lim_tail_mean,dim_2d,/overwrite) + + ; read lat and lon for this tile + lat2d = fltarr(2400l,2400l) + lon2d = fltarr(2400l,2400l) + file_lat_lon = 'MODIS_lat_lon/MODIS_hdres_lat_lon_v'+iiv_str+'_h'+iih_str+'.bin.gz' + openr , lun, file_lat_lon, /get_lun,/compress + readu , lun, lat2d,lon2d + free_lun, lun + + ; if no valid values for snow albedo on this tile, write all missing vals + ind_fit=where(lon2d ge minlon and lon2d le maxlon and $ + lat2d ge minlat and lat2d le maxlat and $ + alb_lim_tail_mean ge 0 and alb_lim_tail_mean le 100 , c_fit) + + ; if no grids with corresponding lat/lons to this tile, skip it + if c_fit eq 0 then continue + + ; if valid grids corresponding to this tile exist then remember them + alb_lim_tail_mean_all=[alb_lim_tail_mean_all,alb_lim_tail_mean[ind_fit]] + lon_all =[lon_all ,lon2d [ind_fit]] + lat_all =[lat_all ,lat2d [ind_fit]] + + endfor ; iih + endfor ; iiv + + ; if there are any valid values, grid them to populate the tile + if n_elements(alb_lim_tail_mean_all) gt 0 then begin + + print, 'Start gridding' + + ; regrid from 2400x2400 sinusoidal grid to 1200x1200 (30 arc sec) regular lat/lon grid (dateline-on-edge, pole-on-edge) + alb_30sec_grid=grid(alb_lim_tail_mean_all,lat_all,lon_all,nlat=1200,nlon=1200, $ + region=[minlat,maxlat,minlon,maxlon],mis_val=mis_val) + + ; scale to range [0.0,1.0] + ind_val=where(alb_30sec_grid ne mis_val,c_val) + if c_val gt 0 then alb_30sec_grid[ind_val]=alb_30sec_grid[ind_val]/100. + + endif + + ; write snow albedo in NetCDF format for the current MODIS tile + ; *** Create a NCDF file + + ; Set up the file & handler + nc_file='data/data_out/snow_alb_all_08_Top'+strtrim(top_alb_limit,2)+ $ + 'th_percentile_MOD10A1.A_30arcsec_'+year_start+'_'+year_end+ $ + '_H'+ih_str+'V'+iv_str+'.nc' + + ; create the file + fid=ncdf_create(nc_file,/CLOBBER,/NETCDF4_FORMAT) ; erese existing file and make a new one + + ; write global attributes + NCDF_ATTPUT, fid, 'N_lon_global' , '43200' , /GLOBAL,/char ; Total number of grids in i-direction + NCDF_ATTPUT, fid, 'N_lat_global' , '21600' , /GLOBAL,/char ; Total number of grids in j-direction + NCDF_ATTPUT, fid, 'i_ind_offset_LL', strtrim(ih*1200l+1,2), /GLOBAL,/char ; GEOS-friendly H grid box + NCDF_ATTPUT, fid, 'j_ind_offset_LL', strtrim(iv*1200l+1,2), /GLOBAL,/char ; GEOS-friendly V grid box + NCDF_ATTPUT, fid, 'CellSize=' , '30arcsec' , /GLOBAL,/char ; grid size + NCDF_ATTPUT, fid, 'CreatedBy' , 'borescan' , /GLOBAL,/char ; user info + spawn,'date',date1 ; get the time form the system + NCDF_ATTPUT, fid, 'Date' , strtrim(date1,2) , /GLOBAL,/char ; time and date + NCDF_ATTPUT, fid, 'Data_Resolution',' 1200 x 1200 ' , /GLOBAL,/char ; resolution of the run + NCDF_ATTPUT, fid, 'Region:' ,'MODIS tile H'+ih_str+'V'+iv_str, /GLOBAL,/char ; covearge + + ;Set Dimensions (there will be two dimentsions: lat and lon, i.e. x and y) + d=indgen(2,/LONG) ; number of dimensions that will be created + d[0]=ncdf_dimdef(fid,'N_lon',1200l) ; number of longitudes + d[1]=ncdf_dimdef(fid,'N_lat',1200l) ; number of latitudes + + ; Define variables to be stored + var_id1=ncdf_vardef(fid,'lon' ,d[0],/FLOAT) + var_id2=ncdf_vardef(fid,'lat' ,d[1],/FLOAT) + var_id3=ncdf_vardef(fid,'Snow_Albedo',d ,/FLOAT) + + ; Change modes (to enter data section/group) + ncdf_control,fid,/ENDEF + + ; Write the data + ncdf_varput,fid,'lat', lat_positions + ncdf_varput,fid,'lon', lon_positions + + ; missing value for MAPL is 1.e+15f + ncdf_varput,fid,'Snow_Albedo', alb_30sec_grid + NCDF_ATTPUT,fid,'Snow_Albedo', "scale_factor", 1.0 + NCDF_ATTPUT,fid,'Snow_Albedo', "offset", 0.0 + NCDF_ATTPUT,fid,'Snow_Albedo', "missing_value", 1.e15 + NCDF_ATTPUT,fid,'Snow_Albedo', "long_name", string("Snow Albedo"),/char + NCDF_ATTPUT,fid,'Snow_Albedo', "units", "-",/char + + ; Close the file & release the handler + ncdf_close,fid + + print, 'Done creating NCDF file' + + endfor ; iv +endfor ;ih + +; values in the SNOW_ALBEDO_DAILY_TILE are as following: + ; 1-100: snow albedo + ; 101: no decision + ; 111: night + ; 125: land + ; 137: inland water + ; 139: ocean + ; 150: cloud + ; 151: cloud detected as snow 250: missing + ; 251: self-shadowing + ; 252: land mask mismatch 253: BRDF failure + ; 254: non-production mask + +; NDSI snow cover general quality value ; values in the SNOW_ALBEDO_DAILY_TILE are as following: + ; 0=best, ; 1–100: snow albedo + ; 1=good, ; 101: no decision + ; 2=ok, ; 111: night + ; 3=poor-not used, ; 125: land + ; 4=other-not used, ; 137: inland water + ; 211=night, ; 139: ocean + ; 239=ocean, ; 150: cloud + ; 255=unusable L1B data or no data ; 151: cloud detected as snow 250: missing + ; 251: self-shadowing + ; 252: land mask mismatch 253: BRDF failure + ; 254: non-production mask + +stop +end diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/grid.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/grid.pro new file mode 100755 index 000000000..4b34d807d --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/grid.pro @@ -0,0 +1,80 @@ +function grid,data,lat,lon,region=rgn,nlat=nlat,nlon=nlon, $ + npts=npts,mis_val=mval + +; function to [re]grid data given lat and lon position of the elements. +; Output is a 2D array on a lat/lon grid, holding means of the elements that fell +; into each gridbox. Default regridded data resolution is 1deg (both lat/lon) + +; Optional arguments include: +; - number of points along lat-direction (nlat) and lon-direction (nlon) - controles the resolution of the output +; - a user-defined missing value (mval) - defines missing value +; - a user-defined region to [re]grid over. [min_lat,mat_lat, min_lon, max_lon] +; - flag '/npts'. Will add the number of elements per grid box to the output + +if (n_elements(rgn) EQ 0) then rgn=[-90.0,90.0,-180.0,180.0] +if (n_elements(nlat) EQ 0) then nlat=180 +if (n_elements(nlon) EQ 0) then nlon=360 +if (n_elements(mval) EQ 0) then mval=-999.0 +if (n_elements(eps) EQ 0) then eps=0.1 +if max(lon) gt 180. then lon[where(lon gt 180.,/null)]=lon[where(lon gt 180.,/null)]-360. +minlat = rgn[0] +maxlat = rgn[1] +minlon = rgn[2] +maxlon = rgn[3] +xinc=(maxlon-minlon)/float(NLON) +yinc=(maxlat-minlat)/float(NLAT) +ngrd=lonarr(NLON,NLAT) +rgrd=fltarr(NLON,NLAT) +sgrd=fltarr(NLON,NLAT) + +; Eliminate missing data + +lat2=lat +lon2=lon +data2=data +m=where(abs(data2-mval) GT eps AND abs(lat2-mval) GT eps AND abs(lon2-mval) GT eps,cnt) +if (cnt GT 0) then begin + lat2=lat[m] + lon2=lon[m] + data2=data[m] +endif + +; Determine if data is inside specified grid region + +xind=long((lon2 - minlon)/double(xinc)) +yind=long((lat2 - minlat)/double(yinc)) +m=where(xind GE 0 AND xind LT NLON AND yind GE 0 AND yind LT NLAT,cnt) + +lat2=lat2[m] +lon2=lon2[m] +data2=data2[m] +xind=xind[m] +yind=yind[m] + +ind=xind+NLON*yind +h=histogram(ind,MIN=0,MAX=long(NLON)*long(NLAT),locations=x) +m=where(h GT 0,mcnt) +for i=0L,mcnt-1 do begin + ind2=where(ind EQ x[m[i]],ncnt) + ngrd[m[i]] = ngrd[m[i]] + ncnt + rgrd[m[i]] = rgrd[m[i]] + total(data2[ind2]) +endfor + +; Compute grid averages + +ind=where(ngrd EQ 0,cnt) +if (cnt GT 0) then rgrd[ind] = mval +ind=where(ngrd GT 0,cnt) +if (cnt GT 0) then rgrd[ind] = rgrd[ind] / float(ngrd[ind]) + +n=1 +if (keyword_set(npts)) then n=n+1 +a=fltarr(NLON,NLAT,n) +a[*,*,0] = rgrd +n=1 +if (keyword_set(npts)) then begin + a[*,*,n] = float(ngrd) + n = n + 1 +endif +return,a +end diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/read_hdf_sd.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/read_hdf_sd.pro new file mode 100755 index 000000000..401a4cd8b --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/read_hdf_sd.pro @@ -0,0 +1,90 @@ +function read_hdf_sd,filename,vars=varname,list=list + +FileID=HDF_SD_START(filename,/READ) +HDF_SD_FILEINFO,FileID,nvar,nattr + +nvarname = n_elements(varname) +if (nvarname EQ 0) then varname='' +a=create_struct('FileName',filename) + +var=strarr(nvar) +typ=strarr(nvar) +ndim=intarr(nvar) +dim=strarr(nvar,10) +varlen = 0 +typlen = 0 +dimnum = 0 +ndmax = 0 +index = 0 +ivar = 0 +while (index LT nvar) do begin + sds_id=HDF_SD_SELECT(FileID,index) + HDF_SD_GETINFO, sds_id, NAME=name, DIMS=dims, NDIMS=ndims, TYPE=type + ind = where(name EQ varname, cnt) + if (cnt GT 0 OR nvarname EQ 0) then begin + if (ndims GT 1 OR dims[0] GT 0) then begin + var[ivar] = name + typ[ivar] = type + ndim[ivar] = ndims + nd=0 + flag=0 + for i=0,ndims-1 do begin + dim[ivar,i] = dims[i] + nd=nd+alog10(dim[ivar,i])+3 + if (dims[i] EQ 0) then flag=1 + endfor + if (flag EQ 0) then begin + if (nd GT ndmax) then ndmax=nd + if (strlen(var[ivar]) GT varlen) then varlen=strlen(var[ivar]) + if (strlen(typ[ivar]) GT typlen) then typlen=strlen(typ[ivar]) + if (ndims GT dimnum) then dimnum=ndims + HDF_SD_GETDATA, sds_id, data + name = var[ivar] + c=['-','.','/','(',')','+'] + for n=0,n_elements(c)-1 do begin + pos=strpos(name,c[n]) + while (pos GE 0) do begin + strput,name,' ',pos + pos=strpos(name,c[n]) + endwhile + endfor + name = strcompress(name,/remove_all) + fchar=strmid(name,0,1) + if (fchar GE '0' AND fchar LE '9') then name='Var_'+name + + svar=string(var[ivar],ivar,format='(a,"_",i3.3)') + a = create_struct(a, name, data) + ivar = ivar + 1 + endif + endif + endif + index = index + 1 +endwhile +nvar = ivar + +if (nvar GT 0) then begin + s1=string(fix(alog10(nvar))+1,format='("i",i1)') + s2='2x,a' + s3=string(typlen+1,format='("a",i2.2)') + ndmax=0 + if (keyword_set(list)) then openw,2,'var.list' + for index=0,nvar-1 do begin + fmt='("SD[",'+s1+',"]: ",'+s3+','+s2 + s4=strarr(ndim[index]) + for i=0,ndim[index]-1 do begin + s4[i] = string(alog10(dim[index,i])+1,format='("i",i1)') + fmt=fmt+',"[",'+s4[i]+',"]"' + endfor + fmt=fmt+')' + if (keyword_set(list)) then $ + printf,2,index,typ[index],var[index],dim[index,0:ndim[index]-1],format=fmt + endfor +endif else begin + a = -1 +endelse + +if (keyword_set(list)) then close,2 +HDF_SD_END, FileID +return,a +end + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/read_mod10a1_hdf.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/read_mod10a1_hdf.pro new file mode 100755 index 000000000..5e30667c8 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/albedo/snow/read_mod10a1_hdf.pro @@ -0,0 +1,30 @@ +function read_mod10a1_hdf, filename + +; code to read MODIS snow albedo files, like: +; MOD10A1.A2020365.h34v08.061.2021006142942.hdf +; or +; MOD10A1.A2022002.h12v11.006.2022004045329.hdf + +; to be used in the stitching process + +; Created March 2022 Biljana Orescanin NASA@SSAI + +vars2read= [ $ + 'NDSI_Snow_Cover' , $ + 'NDSI_Snow_Cover_Basic_QA' , $ + 'NDSI_Snow_Cover_Algorithm_Flags_QA', $ + 'NDSI' , $ + 'Snow_Albedo_Daily_Tile' , $ + 'orbit_pnt' , $ + 'granule_pnt' $ + ] + +; read the hdf file; no scale and offsets, so it's ok to use read_hdf_sd +; note: read_hdf_sd does not allow for choosing which variables to read, +; so I read them all +d=read_hdf_sd(filename) + +return, d + +stop +end diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/CMakeLists.txt new file mode 100644 index 000000000..72c5bf489 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/CMakeLists.txt @@ -0,0 +1,7 @@ +# build without installation + +add_library(m_loss_during_routines m_loss_during_routines.f90) +add_executable(loss_during_day.x loss_during_day.f90) +target_link_libraries(loss_during_day.x m_loss_during_routines) +add_executable(loss_surf_5cm_gensoil.x loss_surf_5cm_gensoil.f90) +target_link_libraries(loss_surf_5cm_gensoil.x m_loss_during_routines) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/PEATCLSM_fitting_CLSM_params.R b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/PEATCLSM_fitting_CLSM_params.R similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/PEATCLSM_fitting_CLSM_params.R rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/PEATCLSM_fitting_CLSM_params.R diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_during_day.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_during_day.f90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_during_day.f90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_during_day.f90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_perday.sh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_perday.sh similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_perday.sh rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_perday.sh diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_perhour.sh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_perhour.sh similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_perhour.sh rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_perhour.sh diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_surf_5cm_gensoil.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_surf_5cm_gensoil.f90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/loss_surf_5cm_gensoil.f90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/loss_surf_5cm_gensoil.f90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/m_loss_during_routines.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/m_loss_during_routines.f90 similarity index 100% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/m_loss_during_routines.f90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/soil/m_loss_during_routines.f90