diff --git a/CMakeLists.txt b/CMakeLists.txt index 9618bd2ea..2f93095f7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -299,14 +299,23 @@ else() src/offline/cbl_model_driver_offline.F90 src/offline/landuse_inout.F90 src/offline/spincasacnp.F90 + src/util/aggregator.F90 src/util/cable_climate_type_mod.F90 src/util/masks_cbl.F90 src/util/cable_array_utils.F90 + src/util/cable_enum.F90 + src/util/cable_grid_reductions.F90 + src/util/cable_timing_utils.F90 src/util/netcdf/cable_netcdf_decomp_util.F90 src/util/netcdf/cable_netcdf.F90 src/util/netcdf/cable_netcdf_internal.F90 src/util/netcdf/cable_netcdf_stub_types.F90 src/util/netcdf/nf90/cable_netcdf_nf90.F90 + src/util/io/output/cable_output_core.F90 + src/util/io/output/cable_output_definitions.F90 + src/util/io/output/cable_output_types.F90 + src/util/io/output/cable_output_utils.F90 + src/util/io/output/cable_output.F90 ) target_link_libraries(cable_common PRIVATE PkgConfig::NETCDF) diff --git a/src/offline/cable_abort.F90 b/src/offline/cable_abort.F90 index ef9962e51..c1bf27343 100644 --- a/src/offline/cable_abort.F90 +++ b/src/offline/cable_abort.F90 @@ -21,7 +21,6 @@ MODULE cable_abort_module USE iso_fortran_env, ONLY: error_unit - USE cable_IO_vars_module, ONLY: check, logn USE cable_mpi_mod, ONLY: mpi_grp_t IMPLICIT NONE @@ -106,72 +105,4 @@ SUBROUTINE nc_abort(ok, message) END SUBROUTINE nc_abort - !============================================================================== - ! - ! Name: range_abort - ! - ! Purpose: Prints an error message and localisation information then stops the - ! code - ! - ! CALLed from: write_output_variable_r1 - ! write_output_variable_r2 - ! - ! MODULEs used: cable_def_types_mod - ! cable_IO_vars_module - ! - !============================================================================== - - SUBROUTINE range_abort(vname, ktau, met, value, var_range, i, xx, yy) - - USE cable_def_types_mod, ONLY: met_type - USE cable_IO_vars_module, ONLY: latitude, longitude, & - landpt, lat_all, lon_all - - ! Input arguments - CHARACTER(LEN=*), INTENT(IN) :: vname - - INTEGER, INTENT(IN) :: & - ktau, & ! time step - i ! tile number along mp - - REAL, INTENT(IN) :: & - xx, & ! coordinates of erroneous grid square - yy ! coordinates of erroneous grid square - - TYPE(met_type), INTENT(IN) :: met ! met data - - REAL(4), INTENT(IN) :: value ! value deemed to be out of range - - REAL, INTENT(IN) :: var_range(2) ! appropriate var range - - INTEGER :: iunit - - IF (check%exit) THEN - iunit = 6 - ELSE - iunit = logn ! warning - END IF - - WRITE (iunit, *) "in SUBR range_abort: Out of range" - WRITE (iunit, *) "for var ", vname ! error from subroutine - - ! patch(i)%latitude, patch(i)%longitude - WRITE (iunit, *) 'Site lat, lon:', xx, yy - WRITE (iunit, *) 'Output timestep', ktau, & - ', or ', met%hod(i), ' hod, ', & - INT(met%doy(i)), 'doy, ', & - INT(met%year(i)) - - WRITE (iunit, *) 'Specified acceptable range (cable_checks.f90):', & - var_range(1), 'to', var_range(2) - - WRITE (iunit, *) 'Value:', value - - IF (check%exit) THEN - STOP - END IF - - END SUBROUTINE range_abort - - !============================================================================== END MODULE cable_abort_module diff --git a/src/offline/cable_checks.F90 b/src/offline/cable_checks.F90 index 5ad32517c..b57537dd0 100644 --- a/src/offline/cable_checks.F90 +++ b/src/offline/cable_checks.F90 @@ -30,8 +30,9 @@ MODULE cable_checks_module ! particular sections of the code - largely for diagnostics/fault finding. ! rh_sh - converts relative to sensible humidity if met file units require it ! - USE cable_IO_vars_module, ONLY: patch - USE cable_abort_module, ONLY: range_abort + USE iso_fortran_env, ONLY: error_unit + USE cable_IO_vars_module, ONLY: patch, check, logn + USE cable_abort_module, ONLY: cable_abort USE cable_def_types_mod USE cable_common_module, ONLY: cable_user @@ -212,6 +213,52 @@ MODULE cable_checks_module CONTAINS + SUBROUTINE range_abort(vname, ktau, met, value, var_range, i, xx, yy) + !! Prints an error message and localisation information then stops the code + + CHARACTER(LEN=*), INTENT(IN) :: vname + + INTEGER, INTENT(IN) :: & + ktau, & ! time step + i ! tile number along mp + + REAL, INTENT(IN) :: & + xx, & ! coordinates of erroneous grid square + yy ! coordinates of erroneous grid square + + TYPE(met_type), INTENT(IN) :: met ! met data + + REAL(4), INTENT(IN) :: value ! value deemed to be out of range + + REAL, INTENT(IN) :: var_range(2) ! appropriate var range + + INTEGER :: iunit + + IF (check%exit) THEN + iunit = error_unit + ELSE + iunit = logn ! warning + END IF + + WRITE (iunit, *) "in SUBR range_abort: Out of range" + WRITE (iunit, *) "for var ", vname ! error from subroutine + + ! patch(i)%latitude, patch(i)%longitude + WRITE (iunit, *) 'Site lat, lon:', xx, yy + WRITE (iunit, *) 'Output timestep', ktau, & + ', or ', met%hod(i), ' hod, ', & + INT(met%doy(i)), 'doy, ', & + INT(met%year(i)) + + WRITE (iunit, *) 'Specified acceptable range (cable_checks.f90):', & + var_range(1), 'to', var_range(2) + + WRITE (iunit, *) 'Value:', value + + IF (check%exit) CALL cable_abort("Aborting...") + + END SUBROUTINE range_abort + SUBROUTINE check_range_d1(vname, parameter_r1, parameter_range, ktau, met) CHARACTER(LEN=*) :: vname diff --git a/src/offline/cable_define_types.F90 b/src/offline/cable_define_types.F90 index df411ca31..c720c918a 100644 --- a/src/offline/cable_define_types.F90 +++ b/src/offline/cable_define_types.F90 @@ -24,7 +24,8 @@ !#define UM_BUILD yes MODULE cable_def_types_mod -USE cable_climate_type_mod, ONLY: climate_type + USE cable_climate_type_mod, ONLY: climate_type + USE aggregator_mod, ONLY: aggregator_real32_1d_t, new_aggregator ! Contains all variables which are not subroutine-internal @@ -531,6 +532,8 @@ MODULE cable_def_types_mod ! vh_js ! !litter thermal conductivity (Wm-2K-1) and vapour diffusivity (m2s-1) REAL(r_2), DIMENSION(:), POINTER :: kthLitt, DvLitt + type(aggregator_real32_1d_t), allocatable :: tscrn_max_daily + type(aggregator_real32_1d_t), allocatable :: tscrn_min_daily END TYPE canopy_type @@ -1186,6 +1189,9 @@ SUBROUTINE alloc_canopy_type(var, mp) ALLOCATE (var % kthLitt(mp)) ALLOCATE (var % DvLitt(mp)) + var%tscrn_max_daily = new_aggregator(source_data=var%tscrn); CALL var%tscrn_max_daily%init(method="max") + var%tscrn_min_daily = new_aggregator(source_data=var%tscrn); CALL var%tscrn_min_daily%init(method="min") + END SUBROUTINE alloc_canopy_type ! ------------------------------------------------------------------------------ @@ -1811,6 +1817,9 @@ SUBROUTINE dealloc_canopy_type(var) DEALLOCATE (var % kthLitt) DEALLOCATE (var % DvLitt) + DEALLOCATE(var%tscrn_max_daily) + DEALLOCATE(var%tscrn_min_daily) + END SUBROUTINE dealloc_canopy_type ! ------------------------------------------------------------------------------ diff --git a/src/offline/cable_io_decomp.F90 b/src/offline/cable_io_decomp.F90 index 84055c1c7..4da99efe6 100644 --- a/src/offline/cable_io_decomp.F90 +++ b/src/offline/cable_io_decomp.F90 @@ -269,40 +269,40 @@ subroutine cable_io_decomp_init(io_decomp) io_decomp%patch_to_x_y_patch_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch, var_shape_x_y_patch, CABLE_NETCDF_INT) io_decomp%patch_to_x_y_patch_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch, var_shape_x_y_patch, CABLE_NETCDF_FLOAT) io_decomp%patch_to_x_y_patch_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch, var_shape_x_y_patch, CABLE_NETCDF_DOUBLE) - io_decomp%patch_soil_to_x_y_patch_soil_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_x_y_soil, CABLE_NETCDF_INT) - io_decomp%patch_soil_to_x_y_patch_soil_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_x_y_soil, CABLE_NETCDF_FLOAT) - io_decomp%patch_soil_to_x_y_patch_soil_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_x_y_soil, CABLE_NETCDF_DOUBLE) - io_decomp%patch_snow_to_x_y_patch_snow_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_x_y_snow, CABLE_NETCDF_INT) - io_decomp%patch_snow_to_x_y_patch_snow_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_x_y_snow, CABLE_NETCDF_FLOAT) - io_decomp%patch_snow_to_x_y_patch_snow_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_x_y_snow, CABLE_NETCDF_DOUBLE) - io_decomp%patch_rad_to_x_y_patch_rad_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_x_y_rad, CABLE_NETCDF_INT) - io_decomp%patch_rad_to_x_y_patch_rad_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_x_y_rad, CABLE_NETCDF_FLOAT) - io_decomp%patch_rad_to_x_y_patch_rad_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_x_y_rad, CABLE_NETCDF_DOUBLE) - io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_x_y_plantcarbon, CABLE_NETCDF_INT) - io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_x_y_plantcarbon, CABLE_NETCDF_FLOAT) - io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_x_y_plantcarbon, CABLE_NETCDF_DOUBLE) - io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_x_y_soilcarbon, CABLE_NETCDF_INT) - io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_x_y_soilcarbon, CABLE_NETCDF_FLOAT) - io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_x_y_soilcarbon, CABLE_NETCDF_DOUBLE) + io_decomp%patch_soil_to_x_y_patch_soil_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_x_y_patch_soil, CABLE_NETCDF_INT) + io_decomp%patch_soil_to_x_y_patch_soil_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_x_y_patch_soil, CABLE_NETCDF_FLOAT) + io_decomp%patch_soil_to_x_y_patch_soil_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_x_y_patch_soil, CABLE_NETCDF_DOUBLE) + io_decomp%patch_snow_to_x_y_patch_snow_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_x_y_patch_snow, CABLE_NETCDF_INT) + io_decomp%patch_snow_to_x_y_patch_snow_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_x_y_patch_snow, CABLE_NETCDF_FLOAT) + io_decomp%patch_snow_to_x_y_patch_snow_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_x_y_patch_snow, CABLE_NETCDF_DOUBLE) + io_decomp%patch_rad_to_x_y_patch_rad_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_x_y_patch_rad, CABLE_NETCDF_INT) + io_decomp%patch_rad_to_x_y_patch_rad_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_x_y_patch_rad, CABLE_NETCDF_FLOAT) + io_decomp%patch_rad_to_x_y_patch_rad_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_x_y_patch_rad, CABLE_NETCDF_DOUBLE) + io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_x_y_patch_plantcarbon, CABLE_NETCDF_INT) + io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_x_y_patch_plantcarbon, CABLE_NETCDF_FLOAT) + io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_x_y_patch_plantcarbon, CABLE_NETCDF_DOUBLE) + io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_int32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_x_y_patch_soilcarbon, CABLE_NETCDF_INT) + io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_real32 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_x_y_patch_soilcarbon, CABLE_NETCDF_FLOAT) + io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_real64 = io_decomp_patch_to_x_y_patch(land_x, land_y, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_x_y_patch_soilcarbon, CABLE_NETCDF_DOUBLE) io_decomp%patch_to_land_patch_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch, var_shape_land_patch, CABLE_NETCDF_INT) io_decomp%patch_to_land_patch_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch, var_shape_land_patch, CABLE_NETCDF_FLOAT) io_decomp%patch_to_land_patch_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch, var_shape_land_patch, CABLE_NETCDF_DOUBLE) - io_decomp%patch_soil_to_land_patch_soil_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_land_soil, CABLE_NETCDF_INT) - io_decomp%patch_soil_to_land_patch_soil_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_land_soil, CABLE_NETCDF_FLOAT) - io_decomp%patch_soil_to_land_patch_soil_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_land_soil, CABLE_NETCDF_DOUBLE) - io_decomp%patch_snow_to_land_patch_snow_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_land_snow, CABLE_NETCDF_INT) - io_decomp%patch_snow_to_land_patch_snow_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_land_snow, CABLE_NETCDF_FLOAT) - io_decomp%patch_snow_to_land_patch_snow_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_land_snow, CABLE_NETCDF_DOUBLE) - io_decomp%patch_rad_to_land_patch_rad_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_land_rad, CABLE_NETCDF_INT) - io_decomp%patch_rad_to_land_patch_rad_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_land_rad, CABLE_NETCDF_FLOAT) - io_decomp%patch_rad_to_land_patch_rad_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_land_rad, CABLE_NETCDF_DOUBLE) - io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_land_plantcarbon, CABLE_NETCDF_INT) - io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_land_plantcarbon, CABLE_NETCDF_FLOAT) - io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_land_plantcarbon, CABLE_NETCDF_DOUBLE) - io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_land_soilcarbon, CABLE_NETCDF_INT) - io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_land_soilcarbon, CABLE_NETCDF_FLOAT) - io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_land_soilcarbon, CABLE_NETCDF_DOUBLE) + io_decomp%patch_soil_to_land_patch_soil_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_land_patch_soil, CABLE_NETCDF_INT) + io_decomp%patch_soil_to_land_patch_soil_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_land_patch_soil, CABLE_NETCDF_FLOAT) + io_decomp%patch_soil_to_land_patch_soil_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soil, var_shape_land_patch_soil, CABLE_NETCDF_DOUBLE) + io_decomp%patch_snow_to_land_patch_snow_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_land_patch_snow, CABLE_NETCDF_INT) + io_decomp%patch_snow_to_land_patch_snow_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_land_patch_snow, CABLE_NETCDF_FLOAT) + io_decomp%patch_snow_to_land_patch_snow_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_snow, var_shape_land_patch_snow, CABLE_NETCDF_DOUBLE) + io_decomp%patch_rad_to_land_patch_rad_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_land_patch_rad, CABLE_NETCDF_INT) + io_decomp%patch_rad_to_land_patch_rad_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_land_patch_rad, CABLE_NETCDF_FLOAT) + io_decomp%patch_rad_to_land_patch_rad_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_rad, var_shape_land_patch_rad, CABLE_NETCDF_DOUBLE) + io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_land_patch_plantcarbon, CABLE_NETCDF_INT) + io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_land_patch_plantcarbon, CABLE_NETCDF_FLOAT) + io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_plantcarbon, var_shape_land_patch_plantcarbon, CABLE_NETCDF_DOUBLE) + io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_int32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_land_patch_soilcarbon, CABLE_NETCDF_INT) + io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_real32 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_land_patch_soilcarbon, CABLE_NETCDF_FLOAT) + io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_real64 = io_decomp_patch_to_land_patch(land_decomp_start, landpt(:)%cstart, landpt(:)%nap, mem_shape_patch_soilcarbon, var_shape_land_patch_soilcarbon, CABLE_NETCDF_DOUBLE) io_decomp%patch_to_patch_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch, var_shape_patch, CABLE_NETCDF_INT) io_decomp%patch_to_patch_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch, var_shape_patch, CABLE_NETCDF_FLOAT) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index cccdcedea..2551cee83 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -91,11 +91,22 @@ MODULE cable_mpimaster compare_consistency_check_values USE cable_mpicommon USE cable_IO_vars_module, ONLY : NO_CHECK + use cable_io_vars_module, only: patch USE cable_io_decomp_mod, ONLY: io_decomp_t, cable_io_decomp_init USE casa_cable USE casa_inout_module USE cable_checks_module, ONLY: constant_check_range USE cable_mpi_mod, ONLY: mpi_grp_t + use cable_output_prototype_v2_mod, only: cable_output_mod_init + use cable_output_prototype_v2_mod, only: cable_output_mod_end + use cable_output_prototype_v2_mod, only: cable_output_register_output_variables + use cable_output_prototype_v2_mod, only: cable_output_profiles_init + use cable_output_prototype_v2_mod, only: cable_output_update + use cable_output_prototype_v2_mod, only: cable_output_write + use cable_output_prototype_v2_mod, only: cable_output_write_parameters + use cable_output_prototype_v2_mod, only: cable_output_write_restart + use cable_output_prototype_v2_mod, only: cable_output_core_outputs + use cable_netcdf_mod, only: cable_netcdf_mod_init, cable_netcdf_mod_end IMPLICIT NONE @@ -334,6 +345,10 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) type(io_decomp_t) :: io_decomp + integer :: start_year + + call cable_netcdf_mod_init(mpi_grp_master) + ! END header ! outer loop - spinup loop no. ktau_tot : @@ -633,6 +648,13 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) ENDIF call cable_io_decomp_init(io_decomp) + + if (.not. casaonly) then + call cable_output_mod_init(io_decomp) + call cable_output_register_output_variables(cable_output_core_outputs(canopy, soil)) + call cable_output_profiles_init() + end if + ! MPI: mostly original serial code follows... ENDIF ! CALL1 @@ -769,6 +791,16 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) !$ CALL POPLUC_set_patchfrac(POPLUC,LUC_EXPT) !$ ENDIF + ! TODO(Sean): this is a hack for determining if the current time step + ! is the last of the month. Better way to do this? + IF(ktau == 1) THEN + !MC - use met%year(1) instead of CABLE_USER%YearStart for non-GSWP forcing and leap years + IF ( TRIM(cable_user%MetType) .EQ. '' ) THEN + start_year = met%year(1) + ELSE + start_year = CABLE_USER%YearStart + ENDIF + END IF IF ( .NOT. CASAONLY ) THEN @@ -787,6 +819,12 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) ENDIF ENDIF + if (ktau > kstart .and. mod(ktau - kstart, ktauday) == 0) then + ! Reset daily aggregators if previous time step was the end of day + call canopy%tscrn_max_daily%reset() + call canopy%tscrn_min_daily%reset() + end if + ! MPI: receive this time step's results from the workers CALL master_receive (ocomm, oktau, recv_ts) ! CALL MPI_Waitall (wnp, recv_req, recv_stats, ierr) @@ -820,6 +858,9 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) ENDIF ENDIF + call canopy%tscrn_max_daily%accumulate() + call canopy%tscrn_min_daily%accumulate() + ELSE IF ( MOD((ktau-kstart+1+koffset),ktauday)==0 ) THEN CALL master_send_input (icomm, casa_dump_ts, iktau ) @@ -862,11 +903,16 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) casamet,ssnow, & rad, bal, air, soil, veg, CSBOLTZ, & CEMLEAF, CEMSOIL ) + if (ktau_tot == kstart) call cable_output_write_parameters(kstart, patch, landpt, met) + call cable_output_update(ktau_tot, dels, leaps, start_year, met) + call cable_output_write(ktau_tot, dels, leaps, start_year, met, patch, landpt) CASE DEFAULT CALL write_output( dels, ktau, met, canopy, casaflux, casapool, & casamet, ssnow, & rad, bal, air, soil, veg, CSBOLTZ, CEMLEAF, CEMSOIL ) - + if (ktau_tot == kstart) call cable_output_write_parameters(kstart, patch, landpt, met) + call cable_output_update(ktau_tot, dels, leaps, start_year, met) + call cable_output_write(ktau_tot, dels, leaps, start_year, met, patch, landpt) END SELECT END IF ENDIF @@ -1053,10 +1099,15 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) CASE ('plum', 'cru', 'gswp', 'gswp3') CALL write_output( dels, ktau_tot, met, canopy, casaflux, casapool, casamet, & ssnow, rad, bal, air, soil, veg, CSBOLTZ, CEMLEAF, CEMSOIL ) + if (ktau_tot == kstart) call cable_output_write_parameters(kstart, patch, landpt, met) + call cable_output_update(ktau_tot, dels, leaps, start_year, met) + call cable_output_write(ktau_tot, dels, leaps, start_year, met, patch, landpt) CASE DEFAULT CALL write_output( dels, ktau, met, canopy, casaflux, casapool, casamet, & ssnow, rad, bal, air, soil, veg, CSBOLTZ, CEMLEAF, CEMSOIL ) - + if (ktau == kstart) call cable_output_write_parameters(kstart, patch, landpt, met) + call cable_output_update(ktau, dels, leaps, start_year, met) + call cable_output_write(ktau, dels, leaps, start_year, met, patch, landpt) END SELECT END IF @@ -1258,6 +1309,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) if(.not.l_landuse) then CALL create_restart( logn, dels, ktau, soil, veg, ssnow, & canopy, rough, rad, bgc, bal, met ) + call cable_output_write_restart(current_time=ktau * dels) endif IF (cable_user%CALL_climate) THEN @@ -1321,7 +1373,9 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU, mpi_grp_master) call landuse_deallocate_mp(cend(mland),ms,msn,nrb,mplant,mlitter,msoil,mwood,lucmp) ENDIF + if (.not. casaonly) call cable_output_mod_end() + call cable_netcdf_mod_end() ! Close met data input file: IF ( TRIM(cable_user%MetType) .NE. "gswp" .AND. & diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 0e1d1f101..70ab29bce 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -194,6 +194,12 @@ SUBROUTINE get_default_params(logn, vegparmnew, LUC_EXPT, mpi_grp) ! count to obtain 'landpt_global', 'max_vegpatches' and 'mp_global' CALL countPatch(nlon, nlat, npatch) + ! TODO(Sean): this a pretty odd place to call init_local_structure_variables, + ! it's definitely a sign we need to pull some things out of load_parameters. + ! This was put here because init_local_structure_variables requires to be + ! called after countPatch is called, as this initialises the rest of the + ! global structure variables, and before allocate_cable_vars, which require + ! the local structure variables to be initialised to allocate cable arrays. CALL init_local_structure_variables(mpi_grp) END SUBROUTINE get_default_params diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 14aa29449..790f825ae 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -85,6 +85,7 @@ MODULE cable_serial patch_type,landpt,& defaultLAI, sdoy, smoy, syear, timeunits, calendar, & NO_CHECK + use cable_io_vars_module, only: patch USE cable_io_decomp_mod, ONLY: io_decomp_t USE cable_io_decomp_mod, ONLY: cable_io_decomp_init USE casa_ncdf_module, ONLY: is_casa_time @@ -112,6 +113,16 @@ MODULE cable_serial ncid_wd,ncid_mask USE cable_output_module, ONLY: create_restart,open_output_file, & write_output,close_output_file + use cable_output_prototype_v2_mod, only: cable_output_mod_init + use cable_output_prototype_v2_mod, only: cable_output_mod_end + use cable_output_prototype_v2_mod, only: cable_output_register_output_variables + use cable_output_prototype_v2_mod, only: cable_output_profiles_init + use cable_output_prototype_v2_mod, only: cable_output_update + use cable_output_prototype_v2_mod, only: cable_output_write + use cable_output_prototype_v2_mod, only: cable_output_write_parameters + use cable_output_prototype_v2_mod, only: cable_output_write_restart + use cable_output_prototype_v2_mod, only: cable_output_core_outputs + use cable_netcdf_mod, only: cable_netcdf_mod_init, cable_netcdf_mod_end USE cable_checks_module, ONLY: constant_check_range USE cable_write_module, ONLY: nullify_write USE cable_IO_vars_module, ONLY: timeunits,calendar @@ -273,10 +284,13 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi real(r_2), dimension(:,:,:), allocatable, save :: patchfrac_new type(io_decomp_t) :: io_decomp + + integer :: start_year ! END header ! INISTUFF + call cable_netcdf_mod_init(mpi_grp) ! outer loop - spinup loop no. ktau_tot : ktau = 0 @@ -462,6 +476,12 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi call cable_io_decomp_init(io_decomp) + if (.not. casaonly) then + call cable_output_mod_init(io_decomp) + call cable_output_register_output_variables(cable_output_core_outputs(canopy, soil)) + call cable_output_profiles_init() + end if + ENDIF ! CALL 1 ! globally (WRT code) accessible kend through USE cable_common_module @@ -572,6 +592,17 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi IF ( CABLE_USER%POPLUC) CALL POPLUC_set_patchfrac(POPLUC,LUC_EXPT) ENDIF + ! TODO(Sean): this is a hack for determining if the current time step + ! is the last of the month. Better way to do this? + IF(ktau == 1) THEN + !MC - use met%year(1) instead of CABLE_USER%YearStart for non-GSWP forcing and leap years + IF ( TRIM(cable_user%MetType) .EQ. '' ) THEN + start_year = met%year(1) + ELSE + start_year = CABLE_USER%YearStart + ENDIF + END IF + IF ( .NOT. CASAONLY ) THEN ! Feedback prognostic vcmax and daily LAI from casaCNP to CABLE @@ -581,6 +612,13 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi IF (l_laiFeedbk.AND.icycle>0) veg%vlai(:) = casamet%glai(:) IF (.NOT. allocated(c1)) ALLOCATE( c1(mp,nrb), rhoch(mp,nrb), xk(mp,nrb) ) + + if (ktau > kstart .and. mod(ktau - kstart, ktauday) == 0) then + ! Reset daily aggregators if previous time step was the end of day + call canopy%tscrn_max_daily%reset() + call canopy%tscrn_min_daily%reset() + end if + ! Call land surface scheme for this timestep, all grid points: CALL cbm( ktau, dels, air, bgc, canopy, met, bal, & rad, rough, soil, ssnow, sum_flux, veg, climate, xk, c1, rhoch ) @@ -595,9 +633,8 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi ssnow%rnof2 = ssnow%rnof2*dels ssnow%runoff = ssnow%runoff*dels - - - + call canopy%tscrn_max_daily%accumulate() + call canopy%tscrn_min_daily%accumulate() ELSE IF ( IS_CASA_TIME("dread", yyyy, ktau, kstart, & koffset, kend, ktauday, logn) ) THEN ! CLN READ FROM FILE INSTEAD ! @@ -722,6 +759,9 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi CALL write_output( dels, ktau, met, canopy, casaflux, casapool, casamet, & ssnow, rad, bal, air, soil, veg, CSBOLTZ, CEMLEAF, CEMSOIL ) END SELECT + if (ktau == kstart) call cable_output_write_parameters(kstart, patch, landpt, met) + call cable_output_update(ktau, dels, leaps, start_year, met) + call cable_output_write(ktau, dels, leaps, start_year, met, patch, landpt) ENDIF @@ -1004,9 +1044,11 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi IF ( .NOT. CASAONLY.and. .not. l_landuse ) THEN ! Write restart file if requested: - IF(output%restart) & + IF(output%restart) then CALL create_restart( logn, dels, ktau, soil, veg, ssnow, & canopy, rough, rad, bgc, bal, met ) + call cable_output_write_restart(current_time=ktau * dels) + end if !mpidiff IF (cable_user%CALL_climate) & CALL WRITE_CLIMATE_RESTART_NC ( climate, ktauday ) @@ -1014,7 +1056,9 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi !--- LN ------------------------------------------[ ENDIF + if (.not. casaonly) call cable_output_mod_end() + call cable_netcdf_mod_end() IF ( TRIM(cable_user%MetType) .NE. "gswp" .AND. & TRIM(cable_user%MetType) .NE. "gswp3" .AND. & diff --git a/src/util/aggregator.F90 b/src/util/aggregator.F90 new file mode 100644 index 000000000..838c164a7 --- /dev/null +++ b/src/util/aggregator.F90 @@ -0,0 +1,614 @@ +module aggregator_mod + + use iso_fortran_env, only: int32, real32, real64 + use cable_abort_module, only: cable_abort + + implicit none + private + + public :: aggregator_t + public :: aggregator_int32_0d_t + public :: aggregator_int32_1d_t + public :: aggregator_int32_2d_t + public :: aggregator_int32_3d_t + public :: aggregator_real32_0d_t + public :: aggregator_real32_1d_t + public :: aggregator_real32_2d_t + public :: aggregator_real32_3d_t + public :: aggregator_real64_0d_t + public :: aggregator_real64_1d_t + public :: aggregator_real64_2d_t + public :: aggregator_real64_3d_t + public :: new_aggregator + + type, abstract :: aggregator_t + integer :: counter = 0 + procedure(accumulate_data), pointer :: accumulate + procedure(reset_data), pointer :: reset + contains + procedure :: init => aggregator_init + procedure :: set_method => aggregator_set_method + procedure :: rank => aggregator_rank + procedure :: shape => aggregator_shape + end type aggregator_t + + abstract interface + subroutine accumulate_data(this) + import aggregator_t + class(aggregator_t), intent(inout) :: this + end subroutine accumulate_data + subroutine reset_data(this) + import aggregator_t + class(aggregator_t), intent(inout) :: this + end subroutine reset_data + end interface + + type, extends(aggregator_t) :: aggregator_int32_0d_t + integer(kind=int32), allocatable :: aggregated_data + integer(kind=int32), pointer :: source_data => null() + end type aggregator_int32_0d_t + + type, extends(aggregator_t) :: aggregator_int32_1d_t + integer(kind=int32), dimension(:), allocatable :: aggregated_data + integer(kind=int32), dimension(:), pointer :: source_data => null() + end type aggregator_int32_1d_t + + type, extends(aggregator_t) :: aggregator_int32_2d_t + integer(kind=int32), dimension(:,:), allocatable :: aggregated_data + integer(kind=int32), dimension(:,:), pointer :: source_data => null() + end type aggregator_int32_2d_t + + type, extends(aggregator_t) :: aggregator_int32_3d_t + integer(kind=int32), dimension(:,:,:), allocatable :: aggregated_data + integer(kind=int32), dimension(:,:,:), pointer :: source_data => null() + end type aggregator_int32_3d_t + + type, extends(aggregator_t) :: aggregator_real32_0d_t + real(kind=real32), allocatable :: aggregated_data + real(kind=real32), pointer :: source_data => null() + end type aggregator_real32_0d_t + + type, extends(aggregator_t) :: aggregator_real32_1d_t + real(kind=real32), dimension(:), allocatable :: aggregated_data + real(kind=real32), dimension(:), pointer :: source_data => null() + end type aggregator_real32_1d_t + + type, extends(aggregator_t) :: aggregator_real32_2d_t + real(kind=real32), dimension(:,:), allocatable :: aggregated_data + real(kind=real32), dimension(:,:), pointer :: source_data => null() + end type aggregator_real32_2d_t + + type, extends(aggregator_t) :: aggregator_real32_3d_t + real(kind=real32), dimension(:,:,:), allocatable :: aggregated_data + real(kind=real32), dimension(:,:,:), pointer :: source_data => null() + end type aggregator_real32_3d_t + + type, extends(aggregator_t) :: aggregator_real64_0d_t + real(kind=real64), allocatable :: aggregated_data + real(kind=real64), pointer :: source_data => null() + end type aggregator_real64_0d_t + + type, extends(aggregator_t) :: aggregator_real64_1d_t + real(kind=real64), dimension(:), allocatable :: aggregated_data + real(kind=real64), dimension(:), pointer :: source_data => null() + end type aggregator_real64_1d_t + + type, extends(aggregator_t) :: aggregator_real64_2d_t + real(kind=real64), dimension(:,:), allocatable :: aggregated_data + real(kind=real64), dimension(:,:), pointer :: source_data => null() + end type aggregator_real64_2d_t + + type, extends(aggregator_t) :: aggregator_real64_3d_t + real(kind=real64), dimension(:,:,:), allocatable :: aggregated_data + real(kind=real64), dimension(:,:,:), pointer :: source_data => null() + end type aggregator_real64_3d_t + + interface new_aggregator + module procedure new_aggregator_int32_0d_t + module procedure new_aggregator_int32_1d_t + module procedure new_aggregator_int32_2d_t + module procedure new_aggregator_int32_3d_t + module procedure new_aggregator_real32_0d + module procedure new_aggregator_real32_1d + module procedure new_aggregator_real32_2d + module procedure new_aggregator_real32_3d + module procedure new_aggregator_real64_0d + module procedure new_aggregator_real64_1d + module procedure new_aggregator_real64_2d + module procedure new_aggregator_real64_3d + end interface + +contains + + subroutine aggregator_init(this, method) + class(aggregator_t), intent(inout) :: this + character(len=*), intent(in) :: method + + select type (this) + type is (aggregator_int32_0d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_int32_1d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_int32_2d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_int32_3d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real32_0d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real32_1d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real32_2d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real32_3d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real64_0d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real64_1d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real64_2d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + type is (aggregator_real64_3d_t) + if (.not. allocated(this%aggregated_data)) allocate(this%aggregated_data, mold=this%source_data) + end select + + call this%set_method(method) + + call this%reset() + + end subroutine aggregator_init + + subroutine aggregator_set_method(this, method) + class(aggregator_t), intent(inout) :: this + character(len=*), intent(in) :: method + + if (method == "mean") then + this%accumulate => mean_accumulate + this%reset => other_reset + elseif (method == "sum") then + this%accumulate => sum_accumulate + this%reset => other_reset + elseif (method == "point") then + this%accumulate => point_accumulate + this%reset => point_reset + elseif (method == "min") then + this%accumulate => min_accumulate + this%reset => min_reset + elseif (method == "max") then + this%accumulate => max_accumulate + this%reset => max_reset + else + call cable_abort("Aggregation method "//method//" is invalid.") + endif + + end subroutine aggregator_set_method + + integer function aggregator_rank(this) + class(aggregator_t), intent(in) :: this + + select type (this) + type is (aggregator_int32_0d_t) + aggregator_rank = 0 + type is (aggregator_int32_1d_t) + aggregator_rank = 1 + type is (aggregator_int32_2d_t) + aggregator_rank = 2 + type is (aggregator_int32_3d_t) + aggregator_rank = 3 + type is (aggregator_real32_0d_t) + aggregator_rank = 0 + type is (aggregator_real32_1d_t) + aggregator_rank = 1 + type is (aggregator_real32_2d_t) + aggregator_rank = 2 + type is (aggregator_real32_3d_t) + aggregator_rank = 3 + type is (aggregator_real64_0d_t) + aggregator_rank = 0 + type is (aggregator_real64_1d_t) + aggregator_rank = 1 + type is (aggregator_real64_2d_t) + aggregator_rank = 2 + type is (aggregator_real64_3d_t) + aggregator_rank = 3 + end select + + end function aggregator_rank + + function aggregator_shape(this) result(agg_shape) + class(aggregator_t), intent(in) :: this + integer, allocatable :: agg_shape(:) + + select type (this) + type is (aggregator_int32_0d_t) + agg_shape = shape(this%source_data) + type is (aggregator_int32_1d_t) + agg_shape = shape(this%source_data) + type is (aggregator_int32_2d_t) + agg_shape = shape(this%source_data) + type is (aggregator_int32_3d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real32_0d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real32_1d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real32_2d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real32_3d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real64_0d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real64_1d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real64_2d_t) + agg_shape = shape(this%source_data) + type is (aggregator_real64_3d_t) + agg_shape = shape(this%source_data) + end select + + end function aggregator_shape + + subroutine mean_accumulate(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_real32_0d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + type is (aggregator_real32_1d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + type is (aggregator_real32_2d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + type is (aggregator_real32_3d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + type is (aggregator_real64_0d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + type is (aggregator_real64_1d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + type is (aggregator_real64_2d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + type is (aggregator_real64_3d_t) + this%aggregated_data = this%aggregated_data + (this%source_data - this%aggregated_data) / (this%counter + 1) + end select + + this%counter = this%counter + 1 + + end subroutine mean_accumulate + + subroutine sum_accumulate(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_int32_0d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_int32_1d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_int32_2d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_int32_3d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real32_0d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real32_1d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real32_2d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real32_3d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real64_0d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real64_1d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real64_2d_t) + this%aggregated_data = this%aggregated_data + this%source_data + type is (aggregator_real64_3d_t) + this%aggregated_data = this%aggregated_data + this%source_data + end select + + this%counter = this%counter + 1 + + end subroutine sum_accumulate + + subroutine point_accumulate(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_int32_0d_t) + this%aggregated_data = this%source_data + type is (aggregator_int32_1d_t) + this%aggregated_data = this%source_data + type is (aggregator_int32_2d_t) + this%aggregated_data = this%source_data + type is (aggregator_int32_3d_t) + this%aggregated_data = this%source_data + type is (aggregator_real32_0d_t) + this%aggregated_data = this%source_data + type is (aggregator_real32_1d_t) + this%aggregated_data = this%source_data + type is (aggregator_real32_2d_t) + this%aggregated_data = this%source_data + type is (aggregator_real32_3d_t) + this%aggregated_data = this%source_data + type is (aggregator_real64_0d_t) + this%aggregated_data = this%source_data + type is (aggregator_real64_1d_t) + this%aggregated_data = this%source_data + type is (aggregator_real64_2d_t) + this%aggregated_data = this%source_data + type is (aggregator_real64_3d_t) + this%aggregated_data = this%source_data + end select + + this%counter = this%counter + 1 + + end subroutine point_accumulate + + subroutine min_accumulate(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_int32_0d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_int32_1d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_int32_2d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_int32_3d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real32_0d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real32_1d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real32_2d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real32_3d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real64_0d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real64_1d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real64_2d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + type is (aggregator_real64_3d_t) + this%aggregated_data = min(this%aggregated_data, this%source_data) + end select + + this%counter = this%counter + 1 + + end subroutine min_accumulate + + subroutine max_accumulate(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_int32_0d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_int32_1d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_int32_2d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_int32_3d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real32_0d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real32_1d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real32_2d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real32_3d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real64_0d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real64_1d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real64_2d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + type is (aggregator_real64_3d_t) + this%aggregated_data = max(this%aggregated_data, this%source_data) + end select + + this%counter = this%counter + 1 + + end subroutine max_accumulate + + subroutine point_reset(this) + class(aggregator_t), intent(inout) :: this + end subroutine point_reset + + subroutine min_reset(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_int32_0d_t) + this%aggregated_data = huge(int(0_int32)) + type is (aggregator_int32_1d_t) + this%aggregated_data = huge(int(0_int32)) + type is (aggregator_int32_2d_t) + this%aggregated_data = huge(int(0_int32)) + type is (aggregator_int32_3d_t) + this%aggregated_data = huge(int(0_int32)) + type is (aggregator_real32_0d_t) + this%aggregated_data = huge(real(0.0_real32)) + type is (aggregator_real32_1d_t) + this%aggregated_data = huge(real(0.0_real32)) + type is (aggregator_real32_2d_t) + this%aggregated_data = huge(real(0.0_real32)) + type is (aggregator_real32_3d_t) + this%aggregated_data = huge(real(0.0_real32)) + type is (aggregator_real64_0d_t) + this%aggregated_data = huge(real(0.0_real64)) + type is (aggregator_real64_1d_t) + this%aggregated_data = huge(real(0.0_real64)) + type is (aggregator_real64_2d_t) + this%aggregated_data = huge(real(0.0_real64)) + type is (aggregator_real64_3d_t) + this%aggregated_data = huge(real(0.0_real64)) + end select + + this%counter = 0 + + end subroutine min_reset + + subroutine max_reset(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_int32_0d_t) + this%aggregated_data = -huge(int(0_int32)) + type is (aggregator_int32_1d_t) + this%aggregated_data = -huge(int(0_int32)) + type is (aggregator_int32_2d_t) + this%aggregated_data = -huge(int(0_int32)) + type is (aggregator_int32_3d_t) + this%aggregated_data = -huge(int(0_int32)) + type is (aggregator_real32_0d_t) + this%aggregated_data = -huge(real(0.0_real32)) + type is (aggregator_real32_1d_t) + this%aggregated_data = -huge(real(0.0_real32)) + type is (aggregator_real32_2d_t) + this%aggregated_data = -huge(real(0.0_real32)) + type is (aggregator_real32_3d_t) + this%aggregated_data = -huge(real(0.0_real32)) + type is (aggregator_real64_0d_t) + this%aggregated_data = -huge(real(0.0_real64)) + type is (aggregator_real64_1d_t) + this%aggregated_data = -huge(real(0.0_real64)) + type is (aggregator_real64_2d_t) + this%aggregated_data = -huge(real(0.0_real64)) + type is (aggregator_real64_3d_t) + this%aggregated_data = -huge(real(0.0_real64)) + end select + + this%counter = 0 + + end subroutine max_reset + + subroutine other_reset(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + type is (aggregator_int32_0d_t) + this%aggregated_data = 0_int32 + type is (aggregator_int32_1d_t) + this%aggregated_data = 0_int32 + type is (aggregator_int32_2d_t) + this%aggregated_data = 0_int32 + type is (aggregator_int32_3d_t) + this%aggregated_data = 0_int32 + type is (aggregator_real32_0d_t) + this%aggregated_data = 0.0_real32 + type is (aggregator_real32_1d_t) + this%aggregated_data = 0.0_real32 + type is (aggregator_real32_2d_t) + this%aggregated_data = 0.0_real32 + type is (aggregator_real32_3d_t) + this%aggregated_data = 0.0_real32 + type is (aggregator_real64_0d_t) + this%aggregated_data = 0.0_real64 + type is (aggregator_real64_1d_t) + this%aggregated_data = 0.0_real64 + type is (aggregator_real64_2d_t) + this%aggregated_data = 0.0_real64 + type is (aggregator_real64_3d_t) + this%aggregated_data = 0.0_real64 + end select + + this%counter = 0 + + end subroutine other_reset + + function new_aggregator_int32_0d_t(source_data) result(agg) + integer(kind=int32), intent(inout), target :: source_data + type(aggregator_int32_0d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_int32_0d_t + + function new_aggregator_int32_1d_t(source_data) result(agg) + integer(kind=int32), dimension(:), intent(inout), target :: source_data + type(aggregator_int32_1d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_int32_1d_t + + function new_aggregator_int32_2d_t(source_data) result(agg) + integer(kind=int32), dimension(:,:), intent(inout), target :: source_data + type(aggregator_int32_2d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_int32_2d_t + + function new_aggregator_int32_3d_t(source_data) result(agg) + integer(kind=int32), dimension(:,:,:), intent(inout), target :: source_data + type(aggregator_int32_3d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_int32_3d_t + + function new_aggregator_real32_0d(source_data) result(agg) + real(kind=real32), intent(inout), target :: source_data + type(aggregator_real32_0d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real32_0d + + function new_aggregator_real32_1d(source_data) result(agg) + real(kind=real32), dimension(:), intent(inout), target :: source_data + type(aggregator_real32_1d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real32_1d + + function new_aggregator_real32_2d(source_data) result(agg) + real(kind=real32), dimension(:,:), intent(inout), target :: source_data + type(aggregator_real32_2d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real32_2d + + function new_aggregator_real32_3d(source_data) result(agg) + real(kind=real32), dimension(:,:,:), intent(inout), target :: source_data + type(aggregator_real32_3d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real32_3d + + function new_aggregator_real64_0d(source_data) result(agg) + real(kind=real64), intent(inout), target :: source_data + type(aggregator_real64_0d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real64_0d + + function new_aggregator_real64_1d(source_data) result(agg) + real(kind=real64), dimension(:), intent(inout), target :: source_data + type(aggregator_real64_1d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real64_1d + + function new_aggregator_real64_2d(source_data) result(agg) + real(kind=real64), dimension(:,:), intent(inout), target :: source_data + type(aggregator_real64_2d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real64_2d + + function new_aggregator_real64_3d(source_data) result(agg) + real(kind=real64), dimension(:,:,:), intent(inout), target :: source_data + type(aggregator_real64_3d_t) :: agg + + agg%source_data => source_data + + end function new_aggregator_real64_3d + +end module diff --git a/src/util/cable_enum.F90 b/src/util/cable_enum.F90 new file mode 100644 index 000000000..7bbb3c18d --- /dev/null +++ b/src/util/cable_enum.F90 @@ -0,0 +1,25 @@ +module cable_enum_mod + implicit none + + type :: cable_enum_t + integer :: value + contains + procedure, private :: cable_enum_eq + generic :: operator(==) => cable_enum_eq + procedure, private :: cable_enum_ne + generic :: operator(/=) => cable_enum_ne + end type + +contains + + elemental logical function cable_enum_eq(this, other) + class(cable_enum_t), intent(in) :: this, other + cable_enum_eq = (this%value == other%value) + end function + + elemental logical function cable_enum_ne(this, other) + class(cable_enum_t), intent(in) :: this, other + cable_enum_ne = (this%value /= other%value) + end function + +end module diff --git a/src/util/cable_grid_reductions.F90 b/src/util/cable_grid_reductions.F90 new file mode 100644 index 000000000..44ae1d6f9 --- /dev/null +++ b/src/util/cable_grid_reductions.F90 @@ -0,0 +1,282 @@ +module cable_grid_reductions_mod + + use iso_fortran_env, only: int32, real32, real64 + + use cable_io_vars_module, only: patch_type, land_type + + implicit none + private + + public :: grid_cell_average + public :: first_patch_in_grid_cell + + interface first_patch_in_grid_cell + module procedure first_patch_in_grid_cell_int32_1d + module procedure first_patch_in_grid_cell_int32_2d + module procedure first_patch_in_grid_cell_int32_3d + module procedure first_patch_in_grid_cell_real32_1d + module procedure first_patch_in_grid_cell_real32_2d + module procedure first_patch_in_grid_cell_real32_3d + module procedure first_patch_in_grid_cell_real64_1d + module procedure first_patch_in_grid_cell_real64_2d + module procedure first_patch_in_grid_cell_real64_3d + end interface + + interface grid_cell_average + module procedure grid_cell_average_real32_1d + module procedure grid_cell_average_real32_2d + module procedure grid_cell_average_real32_3d + module procedure grid_cell_average_real64_1d + module procedure grid_cell_average_real64_2d + module procedure grid_cell_average_real64_3d + end interface + +contains + + subroutine grid_cell_average_real32_1d(input_array, output_array, patch, landpt) + real(kind=real32), intent(in) :: input_array(:) + real(kind=real32), intent(out) :: output_array(:) + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, patch_index + + do land_index = 1, size(output_array) + output_array(land_index) = 0.0_real32 + do patch_index = landpt(land_index)%cstart, landpt(land_index)%cend + output_array(land_index) = output_array(land_index) + & + input_array(patch_index) * patch(patch_index)%frac + end do + end do + + end subroutine + + subroutine grid_cell_average_real32_2d(input_array, output_array, patch, landpt) + real(kind=real32), intent(in) :: input_array(:, :) + real(kind=real32), intent(out) :: output_array(:, :) + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, patch_index, j + + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j) = 0.0_real32 + do patch_index = landpt(land_index)%cstart, landpt(land_index)%cend + output_array(land_index, j) = ( & + output_array(land_index, j) + input_array(patch_index, j) * patch(patch_index)%frac & + ) + end do + end do + end do + + end subroutine + + subroutine grid_cell_average_real32_3d(input_array, output_array, patch, landpt) + real(kind=real32), intent(in) :: input_array(:, :, :) + real(kind=real32), intent(out) :: output_array(:, :, :) + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, patch_index, j, k + + do k = 1, size(output_array, 3) + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j, k) = 0.0_real32 + do patch_index = landpt(land_index)%cstart, landpt(land_index)%cend + output_array(land_index, j, k) = ( & + output_array(land_index, j, k) + & + input_array(patch_index, j, k) * patch(patch_index)%frac & + ) + end do + end do + end do + end do + + end subroutine + + subroutine grid_cell_average_real64_1d(input_array, output_array, patch, landpt) + real(kind=real64), intent(in) :: input_array(:) + real(kind=real64), intent(out) :: output_array(:) + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, patch_index + + do land_index = 1, size(output_array) + output_array(land_index) = 0.0_real64 + do patch_index = landpt(land_index)%cstart, landpt(land_index)%cend + output_array(land_index) = output_array(land_index) + & + input_array(patch_index) * patch(patch_index)%frac + end do + end do + + end subroutine + + subroutine grid_cell_average_real64_2d(input_array, output_array, patch, landpt) + real(kind=real64), intent(in) :: input_array(:, :) + real(kind=real64), intent(out) :: output_array(:, :) + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, patch_index, j + + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j) = 0.0_real64 + do patch_index = landpt(land_index)%cstart, landpt(land_index)%cend + output_array(land_index, j) = ( & + output_array(land_index, j) + input_array(patch_index, j) * patch(patch_index)%frac & + ) + end do + end do + end do + + end subroutine + + subroutine grid_cell_average_real64_3d(input_array, output_array, patch, landpt) + real(kind=real64), intent(in) :: input_array(:, :, :) + real(kind=real64), intent(out) :: output_array(:, :, :) + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, patch_index, j, k + + do k = 1, size(output_array, 3) + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j, k) = 0.0_real64 + do patch_index = landpt(land_index)%cstart, landpt(land_index)%cend + output_array(land_index, j, k) = ( & + output_array(land_index, j, k) + & + input_array(patch_index, j, k) * patch(patch_index)%frac & + ) + end do + end do + end do + end do + + end subroutine + + subroutine first_patch_in_grid_cell_int32_1d(input_array, output_array, landpt) + integer(kind=int32), intent(in) :: input_array(:) + integer(kind=int32), intent(out) :: output_array(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index + + do land_index = 1, size(output_array) + output_array(land_index) = input_array(landpt(land_index)%cstart) + end do + + end subroutine + + subroutine first_patch_in_grid_cell_int32_2d(input_array, output_array, landpt) + integer(kind=int32), intent(in) :: input_array(:, :) + integer(kind=int32), intent(out) :: output_array(:, :) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, j + + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j) = input_array(landpt(land_index)%cstart, j) + end do + end do + + end subroutine + + subroutine first_patch_in_grid_cell_int32_3d(input_array, output_array, landpt) + integer(kind=int32), intent(in) :: input_array(:, :, :) + integer(kind=int32), intent(out) :: output_array(:, :, :) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, j, k + + do k = 1, size(output_array, 3) + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j, k) = input_array(landpt(land_index)%cstart, j, k) + end do + end do + end do + + end subroutine + + subroutine first_patch_in_grid_cell_real32_1d(input_array, output_array, landpt) + real(kind=real32), intent(in) :: input_array(:) + real(kind=real32), intent(out) :: output_array(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index + + do land_index = 1, size(output_array) + output_array(land_index) = input_array(landpt(land_index)%cstart) + end do + + end subroutine + + subroutine first_patch_in_grid_cell_real32_2d(input_array, output_array, landpt) + real(kind=real32), intent(in) :: input_array(:, :) + real(kind=real32), intent(out) :: output_array(:, :) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, j + + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j) = input_array(landpt(land_index)%cstart, j) + end do + end do + + end subroutine + + subroutine first_patch_in_grid_cell_real32_3d(input_array, output_array, landpt) + real(kind=real32), intent(in) :: input_array(:, :, :) + real(kind=real32), intent(out) :: output_array(:, :, :) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, j, k + + do k = 1, size(output_array, 3) + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j, k) = input_array(landpt(land_index)%cstart, j, k) + end do + end do + end do + + end subroutine + + subroutine first_patch_in_grid_cell_real64_1d(input_array, output_array, landpt) + real(kind=real64), intent(in) :: input_array(:) + real(kind=real64), intent(out) :: output_array(:) + type(land_type), intent(in) :: landpt(:) + integer :: land_index + + do land_index = 1, size(output_array) + output_array(land_index) = input_array(landpt(land_index)%cstart) + end do + + end subroutine + + subroutine first_patch_in_grid_cell_real64_2d(input_array, output_array, landpt) + real(kind=real64), intent(in) :: input_array(:, :) + real(kind=real64), intent(out) :: output_array(:, :) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, j + + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j) = input_array(landpt(land_index)%cstart, j) + end do + end do + + end subroutine + + subroutine first_patch_in_grid_cell_real64_3d(input_array, output_array, landpt) + real(kind=real64), intent(in) :: input_array(:, :, :) + real(kind=real64), intent(out) :: output_array(:, :, :) + type(land_type), intent(in) :: landpt(:) + integer :: land_index, j, k + + do k = 1, size(output_array, 3) + do j = 1, size(output_array, 2) + do land_index = 1, size(output_array, 1) + output_array(land_index, j, k) = input_array(landpt(land_index)%cstart, j, k) + end do + end do + end do + + end subroutine + +end module diff --git a/src/util/cable_timing_utils.F90 b/src/util/cable_timing_utils.F90 new file mode 100644 index 000000000..e738035eb --- /dev/null +++ b/src/util/cable_timing_utils.F90 @@ -0,0 +1,64 @@ +module cable_timing_utils_mod + use cable_abort_module, only: cable_abort + use cable_common_module, only: is_leapyear, current_year => CurYear + implicit none + private + + public :: time_step_matches + + integer, parameter :: seconds_per_hour = 3600 + integer, parameter :: hours_per_day = 24 + integer, parameter :: months_in_year = 12 + integer, parameter, dimension(months_in_year) :: & + daysm = [31,28,31,30,31,30,31,31,30,31,30,31], & + daysml = [31,29,31,30,31,30,31,31,30,31,30,31], & + lastday = [31,59,90,120,151,181,212,243,273,304,334,365], & + lastdayl = [31,60,91,121,152,182,213,244,274,305,335,366] + +contains + + function time_step_matches(dels, ktau, frequency, leaps, start_year) result(match) + real, intent(in) :: dels !! Model time step in seconds + integer, intent(in) :: ktau !! Current time step index + character(len=*), intent(in) :: frequency !! Frequency string: 'all', 'daily', 'monthly' + logical, intent(in) :: leaps !! Are we using leap years? + integer, intent(in) :: start_year !! Start year of the simulation + logical :: match + integer :: i + integer :: time_steps_per_interval + integer :: interval_in_hours + integer :: last_day_of_month_in_accumulated_days(months_in_year) ! TODO(Sean): better variable name? + + select case (frequency) + case ('user') + read(frequency(5:7), *) interval_in_hours + time_steps_per_interval = seconds_per_hour * interval_in_hours / int(dels) + match = mod(ktau, time_steps_per_interval) == 0 + case ('all') + match = .true. + case ('daily') + time_steps_per_interval = seconds_per_hour * hours_per_day / int(dels) + match = mod(ktau, time_steps_per_interval) == 0 + case ('monthly') + ! TODO(Sean): is there a better algorithm for monthly matching that doesn't involve looping over years? + last_day_of_month_in_accumulated_days = 0 + do i = start_year, current_year - 1 + if (leaps .and. is_leapyear(i)) then + last_day_of_month_in_accumulated_days = last_day_of_month_in_accumulated_days + 366 + else + last_day_of_month_in_accumulated_days = last_day_of_month_in_accumulated_days + 365 + end if + end do + if (leaps .and. is_leapyear(current_year)) then + last_day_of_month_in_accumulated_days = last_day_of_month_in_accumulated_days + lastdayl + else + last_day_of_month_in_accumulated_days = last_day_of_month_in_accumulated_days + lastday + end if + match = any(int(real(last_day_of_month_in_accumulated_days) * hours_per_day * seconds_per_hour / dels) == ktau) + case default + call cable_abort('Error: unknown frequency "' // trim(adjustl(frequency)) // '"', __FILE__, __LINE__) + end select + + end function + +end module diff --git a/src/util/io/output/cable_output.F90 b/src/util/io/output/cable_output.F90 new file mode 100644 index 000000000..5f3f24b57 --- /dev/null +++ b/src/util/io/output/cable_output.F90 @@ -0,0 +1,27 @@ +module cable_output_prototype_v2_mod + + use cable_output_core_mod, only: cable_output_mod_init + use cable_output_core_mod, only: cable_output_mod_end + use cable_output_core_mod, only: cable_output_register_output_variables + use cable_output_core_mod, only: cable_output_profiles_init + use cable_output_core_mod, only: cable_output_update + use cable_output_core_mod, only: cable_output_write + use cable_output_core_mod, only: cable_output_write_parameters + use cable_output_core_mod, only: cable_output_write_restart + + use cable_output_types_mod, only: cable_output_attribute_t + use cable_output_types_mod, only: cable_output_variable_t + + use cable_output_definitions_mod, only: cable_output_core_outputs => core_outputs + + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_PATCH + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_SOIL + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_SNOW + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_RAD + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_PLANTCARBON + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_SOILCARBON + + implicit none + public + +end module diff --git a/src/util/io/output/cable_output_core.F90 b/src/util/io/output/cable_output_core.F90 new file mode 100644 index 000000000..c8af63505 --- /dev/null +++ b/src/util/io/output/cable_output_core.F90 @@ -0,0 +1,793 @@ +module cable_output_core_mod + + use iso_fortran_env, only: int32, real32, real64 + + use cable_def_types_mod, only: met_type + + use cable_io_vars_module, only: patch_type + use cable_io_vars_module, only: land_type + use cable_io_vars_module, only: output + use cable_io_vars_module, only: check + use cable_io_vars_module, only: ON_TIMESTEP + use cable_io_vars_module, only: ON_WRITE + + use aggregator_mod, only: aggregator_int32_0d_t + use aggregator_mod, only: aggregator_int32_1d_t + use aggregator_mod, only: aggregator_int32_2d_t + use aggregator_mod, only: aggregator_int32_3d_t + use aggregator_mod, only: aggregator_real32_0d_t + use aggregator_mod, only: aggregator_real32_1d_t + use aggregator_mod, only: aggregator_real32_2d_t + use aggregator_mod, only: aggregator_real32_3d_t + use aggregator_mod, only: aggregator_real64_0d_t + use aggregator_mod, only: aggregator_real64_1d_t + use aggregator_mod, only: aggregator_real64_2d_t + use aggregator_mod, only: aggregator_real64_3d_t + + use cable_netcdf_mod, only: cable_netcdf_decomp_t + use cable_netcdf_mod, only: cable_netcdf_file_t + use cable_netcdf_mod, only: cable_netcdf_create_file + use cable_netcdf_mod, only: CABLE_NETCDF_INT + use cable_netcdf_mod, only: CABLE_NETCDF_FLOAT + use cable_netcdf_mod, only: CABLE_NETCDF_DOUBLE + use cable_netcdf_mod, only: CABLE_NETCDF_IOTYPE_CLASSIC + + use cable_abort_module, only: cable_abort + + use cable_checks_module, only: check_range + + use cable_io_decomp_mod, only: io_decomp_t + + use cable_timing_utils_mod, only: time_step_matches + + use cable_grid_reductions_mod, only: grid_cell_average + use cable_grid_reductions_mod, only: first_patch_in_grid_cell + + use cable_output_types_mod, only: cable_output_variable_t + use cable_output_types_mod, only: cable_output_profile_t + use cable_output_types_mod, only: FILL_VALUE_INT32 + use cable_output_types_mod, only: FILL_VALUE_REAL32 + use cable_output_types_mod, only: FILL_VALUE_REAL64 + + use cable_output_utils_mod, only: init_decomp_pointers + use cable_output_utils_mod, only: allocate_grid_reduction_buffers + use cable_output_utils_mod, only: deallocate_grid_reduction_buffers + use cable_output_utils_mod, only: check_invalid_frequency + use cable_output_utils_mod, only: dim_size + use cable_output_utils_mod, only: define_variables + use cable_output_utils_mod, only: set_global_attributes + use cable_output_utils_mod, only: associate_decomp_int32 + use cable_output_utils_mod, only: associate_decomp_real32 + use cable_output_utils_mod, only: associate_decomp_real64 + use cable_output_utils_mod, only: associate_temp_buffer_int32 + use cable_output_utils_mod, only: associate_temp_buffer_real32 + use cable_output_utils_mod, only: associate_temp_buffer_real64 + + use cable_output_definitions_mod, only: coordinate_variables + + implicit none + private + + public :: cable_output_mod_init + public :: cable_output_mod_end + public :: cable_output_register_output_variables + public :: cable_output_profiles_init + public :: cable_output_update + public :: cable_output_write + public :: cable_output_write_parameters + public :: cable_output_write_restart + + type(cable_output_profile_t), allocatable :: global_profile + + type(cable_output_variable_t), allocatable :: registered_output_variables(:) + +contains + + subroutine cable_output_mod_init(io_decomp) + class(io_decomp_t), intent(in) :: io_decomp + class(cable_netcdf_file_t), allocatable :: output_file + + call init_decomp_pointers(io_decomp) + call allocate_grid_reduction_buffers() + + end subroutine + + subroutine cable_output_register_output_variables(output_variables) + type(cable_output_variable_t), dimension(:), intent(in) :: output_variables + integer :: i + + do i = 1, size(output_variables) + associate(output_var => output_variables(i)) + if (all(output_var%reduction_method /= [character(32) :: "none", "grid_cell_average", "first_patch_in_grid_cell"])) then + call cable_abort("Invalid reduction method for variable " // trim(output_var%name), __FILE__, __LINE__) + end if + if (all(output_var%aggregation_method /= [character(32) :: "point", "mean", "max", "min", "sum"])) then + call cable_abort("Invalid aggregation method for variable " // trim(output_var%name), __FILE__, __LINE__) + end if + if (all(output_var%var_type /= [CABLE_NETCDF_INT, CABLE_NETCDF_FLOAT, CABLE_NETCDF_DOUBLE])) then + call cable_abort("Invalid variable type for variable " // trim(output_var%name), __FILE__, __LINE__) + end if + if (count(output_var%name == output_variables(:)%name) > 1) then + call cable_abort("Duplicate variable name found: " // trim(output_var%name), __FILE__, __LINE__) + end if + if (( & + .not. allocated(output_var%data_shape) .and. output_var%aggregator%rank() /= 0 & + ) .or. ( & + allocated(output_var%data_shape) .and. any(dim_size(output_var%data_shape) /= output_var%aggregator%shape()) & + )) then + call cable_abort("Data shape does not match aggregator shape for variable " // trim(output_var%name), __FILE__, __LINE__) + end if + if (output_var%reduction_method /= "none" .and. .not. output_var%distributed) then + call cable_abort("Grid cell reductions require distributed output for variable " // trim(output_var%name), __FILE__, __LINE__) + end if + end associate + end do + + registered_output_variables = output_variables + + end subroutine cable_output_register_output_variables + + subroutine cable_output_mod_end() + + if (allocated(global_profile%output_file)) call global_profile%output_file%close() + + deallocate(global_profile) + + call deallocate_grid_reduction_buffers() + + end subroutine + + subroutine cable_output_profiles_init() + class(cable_netcdf_file_t), allocatable :: output_file + integer :: i + + allocate(cable_output_profile_t::global_profile) + + global_profile%sampling_frequency = output%averaging + + global_profile%file_name = "test_output.nc" ! TODO(Sean): use filename from namelist + + global_profile%output_variables = [ & + coordinate_variables(), & + pack(registered_output_variables, registered_output_variables(:)%active) & + ] + + do i = 1, size(global_profile%output_variables) + associate(output_var => global_profile%output_variables(i)) + call check_invalid_frequency( & + sampling_frequency=global_profile%sampling_frequency, & + accumulation_frequency=output_var%accumulation_frequency, & + var_name=output_var%name, & + file_name=global_profile%file_name & + ) + end associate + end do + + global_profile%output_file = cable_netcdf_create_file(global_profile%file_name, iotype=CABLE_NETCDF_IOTYPE_CLASSIC) + + call define_variables(global_profile%output_file, global_profile%output_variables) + + call set_global_attributes(global_profile) + + call global_profile%output_file%end_def() + + do i = 1, size(global_profile%output_variables) + associate(output_variable => global_profile%output_variables(i)) + call output_variable%aggregator%init(method=output_variable%aggregation_method) + end associate + end do + + end subroutine + + subroutine cable_output_write_parameters(time_index, patch, landpt, met) + integer, intent(in) :: time_index + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + type(met_type), intent(in) :: met + + integer :: i + + do i = 1, size(global_profile%output_variables) + associate(output_variable => global_profile%output_variables(i)) + if (.not. output_variable%parameter) cycle + call check_variable_range(output_variable, time_index, met) + call output_variable%aggregator%accumulate() + call write_variable(output_variable, global_profile%output_file, patch, landpt) + call output_variable%aggregator%reset() + end associate + end do + + end subroutine cable_output_write_parameters + + subroutine cable_output_update(time_index, dels, leaps, start_year, met) + integer, intent(in) :: time_index + real, intent(in) :: dels + logical, intent(in) :: leaps + integer, intent(in) :: start_year + type(met_type), intent(in) :: met + + real :: current_time + integer :: i + + if (check%ranges == ON_TIMESTEP) then + do i = 1, size(global_profile%output_variables) + call check_variable_range(global_profile%output_variables(i), time_index, met) + end do + end if + + do i = 1, size(global_profile%output_variables) + associate(output_variable => global_profile%output_variables(i)) + if (time_step_matches(dels, time_index, output_variable%accumulation_frequency, leaps, start_year)) then + call output_variable%aggregator%accumulate() + end if + end associate + end do + + end subroutine cable_output_update + + subroutine cable_output_write(time_index, dels, leaps, start_year, met, patch, landpt) + integer, intent(in) :: time_index + real, intent(in) :: dels + logical, intent(in) :: leaps + integer, intent(in) :: start_year + type(met_type), intent(in) :: met + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + + real :: current_time + integer :: i + + if (time_step_matches(dels, time_index, global_profile%sampling_frequency, leaps, start_year)) then + + do i = 1, size(global_profile%output_variables) + associate(output_variable => global_profile%output_variables(i)) + if (output_variable%parameter) cycle + if (check%ranges == ON_WRITE) call check_variable_range(output_variable, time_index, met) + call write_variable(output_variable, global_profile%output_file, patch, landpt, global_profile%frame + 1) + call output_variable%aggregator%reset() + end associate + end do + + current_time = time_index * dels + + if (global_profile%sampling_frequency == "all") then + call global_profile%output_file%put_var("time", current_time, start=[global_profile%frame + 1]) + else + call global_profile%output_file%put_var("time", (current_time + global_profile%previous_write_time) / 2.0, start=[global_profile%frame + 1]) + end if + + call global_profile%output_file%put_var("time_bnds", [global_profile%previous_write_time, current_time], start=[1, global_profile%frame + 1]) + + global_profile%previous_write_time = current_time + global_profile%frame = global_profile%frame + 1 + + end if + + end subroutine cable_output_write + + subroutine cable_output_write_restart(current_time) + real, intent(in) :: current_time !! Current simulation time + class(cable_netcdf_file_t), allocatable :: restart_output_file + + type(cable_output_variable_t), allocatable :: restart_variables(:) + integer :: i + + restart_variables = [ & + coordinate_variables(restart=.true.), & + pack(registered_output_variables, registered_output_variables(:)%restart) & + ] + + restart_output_file = cable_netcdf_create_file("test_restart.nc", iotype=CABLE_NETCDF_IOTYPE_CLASSIC) ! TODO(Sean): use filename from namelist + + call define_variables(restart_output_file, restart_variables, restart=.true.) + + call restart_output_file%end_def() + + call restart_output_file%put_var("time", [current_time]) + + do i = 1, size(restart_variables) + call write_variable(restart_variables(i), restart_output_file, restart=.true.) + end do + + call restart_output_file%close() + + end subroutine cable_output_write_restart + + subroutine check_variable_range(output_variable, time_index, met) + type(cable_output_variable_t), intent(in) :: output_variable + integer, intent(in) :: time_index + type(met_type), intent(in) :: met + + select type (aggregator => output_variable%aggregator) + type is (aggregator_int32_0d_t) + ! TODO(Sean): implement range checking for integer types + type is (aggregator_int32_1d_t) + ! TODO(Sean): implement range checking for integer types + type is (aggregator_int32_2d_t) + ! TODO(Sean): implement range checking for integer types + type is (aggregator_int32_3d_t) + ! TODO(Sean): implement range checking for integer types + type is (aggregator_real32_0d_t) + ! TODO(Sean): implement range checking for scalars + type is (aggregator_real32_1d_t) + call check_range(output_variable%name, aggregator%source_data, output_variable%range, time_index, met) + type is (aggregator_real32_2d_t) + call check_range(output_variable%name, aggregator%source_data, output_variable%range, time_index, met) + type is (aggregator_real32_3d_t) + call check_range(output_variable%name, aggregator%source_data, output_variable%range, time_index, met) + type is (aggregator_real64_0d_t) + ! TODO(Sean): implement range checking for double precision types + type is (aggregator_real64_1d_t) + ! TODO(Sean): implement range checking for double precision types + type is (aggregator_real64_2d_t) + ! TODO(Sean): implement range checking for double precision types + type is (aggregator_real64_3d_t) + ! TODO(Sean): implement range checking for double precision types + class default + call cable_abort("Unexpected aggregator type", __FILE__, __LINE__) + end select + + end subroutine check_variable_range + + subroutine write_variable(output_variable, output_file, patch, landpt, frame, restart) + type(cable_output_variable_t), intent(inout), target :: output_variable + class(cable_netcdf_file_t), intent(inout) :: output_file + type(patch_type), intent(in), optional :: patch(:) + type(land_type), intent(in), optional :: landpt(:) + integer, intent(in), optional :: frame + logical, intent(in), optional :: restart + + class(cable_netcdf_decomp_t), pointer :: decomp + integer :: i, ndims + logical :: restart_local + + integer(kind=int32), pointer :: write_buffer_int32_0d + integer(kind=int32), pointer :: write_buffer_int32_1d(:) + integer(kind=int32), pointer :: write_buffer_int32_2d(:, :) + integer(kind=int32), pointer :: write_buffer_int32_3d(:, :, :) + real(kind=real32), pointer :: write_buffer_real32_0d + real(kind=real32), pointer :: write_buffer_real32_1d(:) + real(kind=real32), pointer :: write_buffer_real32_2d(:, :) + real(kind=real32), pointer :: write_buffer_real32_3d(:, :, :) + real(kind=real64), pointer :: write_buffer_real64_0d + real(kind=real64), pointer :: write_buffer_real64_1d(:) + real(kind=real64), pointer :: write_buffer_real64_2d(:, :) + real(kind=real64), pointer :: write_buffer_real64_3d(:, :, :) + + decomp => null() + + write_buffer_int32_0d => null() + write_buffer_int32_1d => null() + write_buffer_int32_2d => null() + write_buffer_int32_3d => null() + write_buffer_real32_0d => null() + write_buffer_real32_1d => null() + write_buffer_real32_2d => null() + write_buffer_real32_3d => null() + write_buffer_real64_0d => null() + write_buffer_real64_1d => null() + write_buffer_real64_2d => null() + write_buffer_real64_3d => null() + + restart_local = .false. + if (present(restart)) restart_local = restart + + if (.not. restart_local .and. output_variable%reduction_method /= "none") then + if (.not. present(patch) .or. .not. present(landpt)) then + call cable_abort("Optional arguments patch and landpt must be present for grid reductions", __FILE__, __LINE__) + end if + end if + + select type (aggregator => output_variable%aggregator) + type is (aggregator_int32_0d_t) + if (output_variable%reduction_method /= "none") then + call cable_abort("Grid cell reductions are not supported for scalar variables", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call cable_abort("Distributed writes are not supported for scalar variables", __FILE__, __LINE__) + end if + write_buffer_int32_0d => aggregator%aggregated_data + if (restart_local) write_buffer_int32_0d => aggregator%source_data + if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_0d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_0d) + end if + type is (aggregator_int32_1d_t) + if (restart_local) then + write_buffer_int32_1d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_int32_1d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call cable_abort("Reduction method grid_cell_average is not supported for integer variables", __FILE__, __LINE__) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_int32(output_variable, temp_buffer_int32_1d=write_buffer_int32_1d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_int32_1d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_int32(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_int32_1d, & + decomp=decomp, & + fill_value=FILL_VALUE_INT32, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_1d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_1d) + end if + type is (aggregator_int32_2d_t) + if (restart_local) then + write_buffer_int32_2d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_int32_2d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call cable_abort("Reduction method grid_cell_average is not supported for integer variables", __FILE__, __LINE__) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_int32(output_variable, temp_buffer_int32_2d=write_buffer_int32_2d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_int32_2d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_int32(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_int32_2d, & + decomp=decomp, & + fill_value=FILL_VALUE_INT32, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_2d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_2d) + end if + type is (aggregator_int32_3d_t) + if (restart_local) then + write_buffer_int32_3d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_int32_3d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call cable_abort("Reduction method grid_cell_average is not supported for integer variables", __FILE__, __LINE__) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_int32(output_variable, temp_buffer_int32_3d=write_buffer_int32_3d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_int32_3d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_int32(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_int32_3d, & + decomp=decomp, & + fill_value=FILL_VALUE_INT32, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_3d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_int32_3d) + end if + type is (aggregator_real32_0d_t) + if (output_variable%reduction_method /= "none") then + call cable_abort("Grid cell reductions are not supported for scalar variables", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call cable_abort("Distributed writes are not supported for scalar variables", __FILE__, __LINE__) + end if + write_buffer_real32_0d => aggregator%aggregated_data + if (restart_local) write_buffer_real32_0d => aggregator%source_data + if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_0d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_0d) + end if + type is (aggregator_real32_1d_t) + if (restart_local) then + write_buffer_real32_1d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_real32_1d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call associate_temp_buffer_real32(output_variable, temp_buffer_real32_1d=write_buffer_real32_1d) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real32_1d, & + landpt=landpt, & + patch=patch) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_real32(output_variable, temp_buffer_real32_1d=write_buffer_real32_1d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real32_1d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_real32(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_real32_1d, & + decomp=decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_1d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_1d) + end if + type is (aggregator_real32_2d_t) + if (restart_local) then + write_buffer_real32_2d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_real32_2d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call associate_temp_buffer_real32(output_variable, temp_buffer_real32_2d=write_buffer_real32_2d) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real32_2d, & + landpt=landpt, & + patch=patch) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_real32(output_variable, temp_buffer_real32_2d=write_buffer_real32_2d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real32_2d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_real32(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_real32_2d, & + decomp=decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_2d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_2d) + end if + type is (aggregator_real32_3d_t) + if (restart_local) then + write_buffer_real32_3d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_real32_3d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call associate_temp_buffer_real32(output_variable, temp_buffer_real32_3d=write_buffer_real32_3d) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real32_3d, & + landpt=landpt, & + patch=patch) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_real32(output_variable, temp_buffer_real32_3d=write_buffer_real32_3d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real32_3d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_real32(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_real32_3d, & + decomp=decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_3d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real32_3d) + end if + type is (aggregator_real64_0d_t) + if (output_variable%reduction_method /= "none") then + call cable_abort("Grid cell reductions are not supported for scalar variables", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call cable_abort("Distributed writes are not supported for scalar variables", __FILE__, __LINE__) + end if + write_buffer_real64_0d => aggregator%aggregated_data + if (restart_local) write_buffer_real64_0d => aggregator%source_data + if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_0d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_0d) + end if + type is (aggregator_real64_1d_t) + if (restart_local) then + write_buffer_real64_1d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_real64_1d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call associate_temp_buffer_real64(output_variable, temp_buffer_real64_1d=write_buffer_real64_1d) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real64_1d, & + landpt=landpt, & + patch=patch) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_real64(output_variable, temp_buffer_real64_1d=write_buffer_real64_1d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real64_1d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_real64(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_real64_1d, & + decomp=decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_1d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_1d) + end if + type is (aggregator_real64_2d_t) + if (restart_local) then + write_buffer_real64_2d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_real64_2d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call associate_temp_buffer_real64(output_variable, temp_buffer_real64_2d=write_buffer_real64_2d) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real64_2d, & + landpt=landpt, & + patch=patch) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_real64(output_variable, temp_buffer_real64_2d=write_buffer_real64_2d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real64_2d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_real64(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_real64_2d, & + decomp=decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_2d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_2d) + end if + type is (aggregator_real64_3d_t) + if (restart_local) then + write_buffer_real64_3d => aggregator%source_data + else if (output_variable%reduction_method == "none") then + write_buffer_real64_3d => aggregator%aggregated_data + else if (output_variable%reduction_method == "grid_cell_average") then + call associate_temp_buffer_real64(output_variable, temp_buffer_real64_3d=write_buffer_real64_3d) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real64_3d, & + landpt=landpt, & + patch=patch) + else if (output_variable%reduction_method == "first_patch_in_grid_cell") then + call associate_temp_buffer_real64(output_variable, temp_buffer_real64_3d=write_buffer_real64_3d) + call first_patch_in_grid_cell( & + input_array=aggregator%aggregated_data, & + output_array=write_buffer_real64_3d, & + landpt=landpt) + else + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end if + if (output_variable%distributed) then + call associate_decomp_real64(output_variable, decomp, restart=restart_local) + call output_file%write_darray( & + var_name=output_variable%name, & + values=write_buffer_real64_3d, & + decomp=decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=frame) + else if (present(frame)) then + call output_file%inq_var_ndims(output_variable%name, ndims) + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_3d, & + start=[(1, i = 1, ndims - 1), frame]) + else + call output_file%put_var( & + var_name=output_variable%name, & + values=write_buffer_real64_3d) + end if + class default + call cable_abort("Unexpected aggregator type", __FILE__, __LINE__) + end select + + end subroutine write_variable + +end module diff --git a/src/util/io/output/cable_output_definitions.F90 b/src/util/io/output/cable_output_definitions.F90 new file mode 100644 index 000000000..985c18178 --- /dev/null +++ b/src/util/io/output/cable_output_definitions.F90 @@ -0,0 +1,339 @@ +module cable_output_definitions_mod + use cable_def_types_mod, only: canopy_type + use cable_def_types_mod, only: soil_parameter_type + use cable_def_types_mod, only: mvtype, mstype + + use cable_netcdf_mod, only: CABLE_NETCDF_INT + use cable_netcdf_mod, only: CABLE_NETCDF_FLOAT + + use aggregator_mod, only: new_aggregator + + use cable_io_vars_module, only: output, patchout + use cable_io_vars_module, only: landpt_global + use cable_io_vars_module, only: patch + use cable_io_vars_module, only: lat_all, lon_all + use cable_io_vars_module, only: latitude, longitude + use cable_io_vars_module, only: metgrid + + use cable_output_types_mod, only: cable_output_variable_t + use cable_output_types_mod, only: attribute => cable_output_attribute_t + use cable_output_types_mod, only: DIM_PATCH => CABLE_OUTPUT_DIM_PATCH + use cable_output_types_mod, only: DIM_SOIL => CABLE_OUTPUT_DIM_SOIL + use cable_output_types_mod, only: DIM_RAD => CABLE_OUTPUT_DIM_RAD + use cable_output_types_mod, only: DIM_LAND => CABLE_OUTPUT_DIM_LAND + use cable_output_types_mod, only: DIM_LAND_GLOBAL => CABLE_OUTPUT_DIM_LAND_GLOBAL + use cable_output_types_mod, only: DIM_X => CABLE_OUTPUT_DIM_X + use cable_output_types_mod, only: DIM_Y => CABLE_OUTPUT_DIM_Y + + use cable_output_utils_mod, only: requires_x_y_output_grid + use cable_output_utils_mod, only: requires_land_output_grid + + use cable_checks_module, only: ranges + + use cable_abort_module, only: cable_abort + + implicit none + private + + public :: core_outputs + public :: coordinate_variables + +contains + + function core_outputs(canopy, soil) result(output_variables) + type(canopy_type), intent(inout) :: canopy + type(soil_parameter_type), intent(in) :: soil + + type(cable_output_variable_t), allocatable :: output_variables(:) + + output_variables = [ & + cable_output_variable_t( & + name="isoil", & + data_shape=[DIM_PATCH], & + var_type=CABLE_NETCDF_INT, & + range=ranges%isoil, & + active=output%isoil, & + patchout=output%patch .or. patchout%isoil, & + reduction_method="first_patch_in_grid_cell", & + aggregation_method="point", & + parameter=.true., & + restart=.true., & + aggregator=new_aggregator(soil%isoilm), & + metadata=[ & + attribute("units", "-"), & + attribute("long_name", "Soil type") & + ] & + ), & + cable_output_variable_t( & + name="swilt", & + data_shape=[DIM_PATCH], & + var_type=CABLE_NETCDF_FLOAT, & + range=ranges%swilt, & + active=output%swilt, & + patchout=output%patch .or. patchout%swilt, & + reduction_method="grid_cell_average", & + aggregation_method="point", & + parameter=.true., & + aggregator=new_aggregator(soil%swilt), & + metadata=[ & + attribute("units", "1"), & + attribute("long_name", "") & + ] & + ), & + cable_output_variable_t( & + name="albsoil", & + data_shape=[DIM_PATCH, DIM_RAD], & + var_type=CABLE_NETCDF_FLOAT, & + range=ranges%albsoil, & + active=output%albsoil, & + patchout=output%patch .or. patchout%albsoil, & + reduction_method="first_patch_in_grid_cell", & + aggregation_method="point", & + parameter=.true., & + aggregator=new_aggregator(soil%albsoil), & + metadata=[ & + attribute("units", "1"), & + attribute("long_name", "") & + ] & + ), & + cable_output_variable_t( & + name="Qh", & + data_shape=[DIM_PATCH], & + var_type=CABLE_NETCDF_FLOAT, & + range=ranges%Qh, & + active=output%Qh, & + patchout=output%patch .or. patchout%Qh, & + reduction_method="grid_cell_average", & + aggregation_method="mean", & + aggregator=new_aggregator(canopy%fh), & + metadata=[ & + attribute("units", "W/m^2"), & + attribute("long_name", "Surface sensible heat flux") & + ] & + ), & + cable_output_variable_t( & + name="Tmx", & + data_shape=[DIM_PATCH], & + var_type=CABLE_NETCDF_FLOAT, & + active=output%Tex .and. output%averaging == "monthly", & + patchout=output%patch .or. patchout%Tex, & + reduction_method="grid_cell_average", & + aggregation_method="mean", & + range=ranges%Tscrn, & + accumulation_frequency="daily", & + aggregator=new_aggregator(canopy%tscrn_max_daily%aggregated_data), & + metadata=[ & + attribute("units", "oC"), & + attribute("long_name", "averaged daily maximum screen-level T") & + ] & + ), & + cable_output_variable_t( & + name="nap", & + data_shape=[DIM_LAND_GLOBAL], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-huge(0.0), huge(0.0)], & + active=.false., & + restart=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(landpt_global(:)%nap), & + metadata=[ & + attribute("units", ""), & + attribute("long_name", "") & + ] & + ), & + cable_output_variable_t( & + name="patchfrac", & + data_shape=[DIM_PATCH], & + var_type=CABLE_NETCDF_FLOAT, & + range=[0.0, 1.0], & + active=.false., & + restart=.true., & + distributed=.true., & + aggregation_method="point", & + aggregator=new_aggregator(patch(:)%frac), & + metadata=[ & + attribute("units", ""), & + attribute("long_name", "Fraction of vegetated grid cell area occupied by a vegetation/soil patch") & + ] & + ), & + cable_output_variable_t( & + name="mvtype", & + var_type=CABLE_NETCDF_FLOAT, & + range=[-huge(0.0), huge(0.0)], & + active=.false., & + restart=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(mvtype), & + metadata=[ & + attribute("units", ""), & + attribute("long_name", "Number of vegetation types") & + ] & + ), & + cable_output_variable_t( & + name="mstype", & + var_type=CABLE_NETCDF_FLOAT, & + range=[-huge(0.0), huge(0.0)], & + active=.false., & + restart=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(mstype), & + metadata=[ & + attribute("units", ""), & + attribute("long_name", "Number of soil types") & + ] & + ) & + ] + + end function core_outputs + + function coordinate_variables(restart) result(output_variables) + logical, intent(in), optional :: restart + + type(cable_output_variable_t), allocatable :: output_variables(:) + + if (present(restart)) then; if (restart) then + output_variables = [ & + cable_output_variable_t( & + name="latitude", & + data_shape=[DIM_LAND_GLOBAL], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-90.0, 90.0], & + active=.false., & + restart=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(latitude), & + metadata=[ & + attribute("units", "degrees_north"), & + attribute("long_name", "") & + ] & + ), & + cable_output_variable_t( & + name="longitude", & + data_shape=[DIM_LAND_GLOBAL], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-huge(0.0), huge(0.0)], & ! TODO(Sean): this depends on the met forcing input? + active=.false., & + restart=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(longitude), & + metadata=[ & + attribute("units", "degrees_east"), & + attribute("long_name", "") & + ] & + ) & + ] + return + end if; end if + + if (requires_x_y_output_grid(output%grid, metgrid)) then + output_variables = [ & + ! Define latitude and longitude variable (ALMA): + cable_output_variable_t( & + name="latitude", & + data_shape=[DIM_X, DIM_Y], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-90.0, 90.0], & + active=.true., & + parameter=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(lat_all), & + metadata=[ & + attribute("units", "degrees_north"), & + attribute("long_name", "") & + ] & + ), & + cable_output_variable_t( & + name="longitude", & + data_shape=[DIM_X, DIM_Y], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-huge(0.0), huge(0.0)], & ! TODO(Sean): this depends on the met forcing input? + active=.true., & + parameter=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(lon_all), & + metadata=[ & + attribute("units", "degrees_east"), & + attribute("long_name", "") & + ] & + ), & + ! Write "cordinate variables" to enable reading by GrADS: + cable_output_variable_t( & + name="x", & + data_shape=[DIM_X], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-huge(0.0), huge(0.0)], & ! TODO(Sean): this depends on the met forcing input? + active=.true., & + parameter=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(lon_all(:, 1)), & + metadata=[ & + attribute("units", "degrees_east"), & + attribute("long_name", ""), & + attribute("comment", "x coordinate variable for GrADS compatibility") & + ] & + ), & + cable_output_variable_t( & + name="y", & + data_shape=[DIM_Y], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-90.0, 90.0], & ! TODO(Sean): this depends on the met forcing input? + active=.true., & + parameter=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(lat_all(1, :)), & + metadata=[ & + attribute("units", "degrees_north"), & + attribute("long_name", ""), & + attribute("comment", "y coordinate variable for GrADS compatibility") & + ] & + ) & + ] + else if (requires_land_output_grid(output%grid, metgrid)) then + output_variables = [ & + cable_output_variable_t( & + name="local_lat", & + data_shape=[DIM_LAND_GLOBAL], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-90.0, 90.0], & + active=requires_land_output_grid(output%grid, metgrid), & + parameter=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(latitude), & + metadata=[ & + attribute("units", "degrees_north"), & + attribute("long_name", "") & + ] & + ), & + cable_output_variable_t( & + name="local_lon", & + data_shape=[DIM_LAND_GLOBAL], & + var_type=CABLE_NETCDF_FLOAT, & + range=[-huge(0.0), huge(0.0)], & ! TODO(Sean): this depends on the met forcing input? + active=requires_land_output_grid(output%grid, metgrid), & + parameter=.true., & + distributed=.false., & + aggregation_method="point", & + aggregator=new_aggregator(longitude), & + metadata=[ & + attribute("units", "degrees_east"), & + attribute("long_name", "") & + ] & + ) & + ] + else + call cable_abort("Unable to determine coordinate variables for output grid.", __FILE__, __LINE__) + end if + + end function + +end module diff --git a/src/util/io/output/cable_output_types.F90 b/src/util/io/output/cable_output_types.F90 new file mode 100644 index 000000000..4d0904777 --- /dev/null +++ b/src/util/io/output/cable_output_types.F90 @@ -0,0 +1,64 @@ +module cable_output_types_mod + + use iso_fortran_env, only: int32, real32, real64 + + use aggregator_mod, only: aggregator_t + + use cable_netcdf_mod, only: cable_netcdf_file_t + + use cable_enum_mod, only: cable_enum_t + + implicit none + private + + type, extends(cable_enum_t), public :: cable_output_dim_t + end type + + type, public :: cable_output_attribute_t + character(len=64) :: name + character(len=256) :: value + end type + + type, public :: cable_output_variable_t + character(len=64) :: name + character(len=64) :: accumulation_frequency = "all" + character(len=64) :: reduction_method = "none" + character(len=64) :: aggregation_method + logical :: active + logical :: parameter = .false. + logical :: distributed = .true. + logical :: restart = .false. + logical :: patchout = .false. + integer :: var_type + real, dimension(2) :: range + type(cable_output_dim_t), allocatable :: data_shape(:) + class(aggregator_t), allocatable :: aggregator + type(cable_output_attribute_t), allocatable :: metadata(:) + end type + + type, public :: cable_output_profile_t + real :: previous_write_time = 0.0 + integer :: frame = 0 + character(len=64) :: sampling_frequency + character(len=256) :: file_name + class(cable_netcdf_file_t), allocatable :: output_file + type(cable_output_variable_t), allocatable :: output_variables(:) + type(cable_output_attribute_t), allocatable :: metadata(:) + end type + + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_PATCH = cable_output_dim_t(0) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_SOIL = cable_output_dim_t(1) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_SNOW = cable_output_dim_t(2) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_RAD = cable_output_dim_t(3) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_PLANTCARBON = cable_output_dim_t(4) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_SOILCARBON = cable_output_dim_t(5) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_LAND = cable_output_dim_t(6) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_LAND_GLOBAL = cable_output_dim_t(7) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_X = cable_output_dim_t(8) + type(cable_output_dim_t), parameter, public :: CABLE_OUTPUT_DIM_Y = cable_output_dim_t(9) + + integer(kind=int32), parameter, public :: FILL_VALUE_INT32 = -9999999_int32 + real(kind=real32), parameter, public :: FILL_VALUE_REAL32 = -1.0e+33_real32 + real(kind=real64), parameter, public :: FILL_VALUE_REAL64 = -1.0e+33_real64 + +end module diff --git a/src/util/io/output/cable_output_utils.F90 b/src/util/io/output/cable_output_utils.F90 new file mode 100644 index 000000000..5363494fa --- /dev/null +++ b/src/util/io/output/cable_output_utils.F90 @@ -0,0 +1,988 @@ +module cable_output_utils_mod + + use iso_fortran_env, only: int32, real32, real64 + + use cable_common_module, only: filename + + use cable_def_types_mod, only: mp + use cable_def_types_mod, only: mp_global + use cable_def_types_mod, only: mland + use cable_def_types_mod, only: mland_global + use cable_def_types_mod, only: ms + use cable_def_types_mod, only: msn + use cable_def_types_mod, only: nrb + use cable_def_types_mod, only: ncs + use cable_def_types_mod, only: ncp + + use cable_io_vars_module, only: output + use cable_io_vars_module, only: metgrid + use cable_io_vars_module, only: xdimsize + use cable_io_vars_module, only: ydimsize + use cable_io_vars_module, only: max_vegpatches + use cable_io_vars_module, only: timeunits + use cable_io_vars_module, only: time_coord + use cable_io_vars_module, only: calendar + + use cable_abort_module, only: cable_abort + + use cable_netcdf_mod, only: cable_netcdf_decomp_t + use cable_netcdf_mod, only: cable_netcdf_file_t + use cable_netcdf_mod, only: MAX_LEN_DIM => CABLE_NETCDF_MAX_STR_LEN_DIM + use cable_netcdf_mod, only: CABLE_NETCDF_UNLIMITED + use cable_netcdf_mod, only: CABLE_NETCDF_INT + use cable_netcdf_mod, only: CABLE_NETCDF_FLOAT + use cable_netcdf_mod, only: CABLE_NETCDF_DOUBLE + + use cable_io_decomp_mod, only: io_decomp_t + + use cable_output_types_mod, only: cable_output_dim_t + use cable_output_types_mod, only: cable_output_variable_t + use cable_output_types_mod, only: cable_output_profile_t + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_PATCH + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_SOIL + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_SNOW + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_RAD + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_PLANTCARBON + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_SOILCARBON + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_LAND + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_LAND_GLOBAL + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_X + use cable_output_types_mod, only: CABLE_OUTPUT_DIM_Y + use cable_output_types_mod, only: FILL_VALUE_INT32 + use cable_output_types_mod, only: FILL_VALUE_REAL32 + use cable_output_types_mod, only: FILL_VALUE_REAL64 + + implicit none + private + + public :: init_decomp_pointers + public :: allocate_grid_reduction_buffers + public :: deallocate_grid_reduction_buffers + public :: requires_x_y_output_grid + public :: requires_land_output_grid + public :: check_invalid_frequency + public :: dim_size + public :: infer_dim_names + public :: define_variables + public :: set_global_attributes + public :: associate_decomp_int32 + public :: associate_decomp_real32 + public :: associate_decomp_real64 + public :: associate_temp_buffer_int32 + public :: associate_temp_buffer_real32 + public :: associate_temp_buffer_real64 + + ! Decomposition pointers for each variable class and data type + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_soil_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_soil_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_soil_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_snow_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_snow_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_snow_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_rad_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_rad_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_rad_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_plantcarbon_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_plantcarbon_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_plantcarbon_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_soilcarbon_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_soilcarbon_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_soilcarbon_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_soil_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_soil_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_soil_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_snow_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_snow_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_snow_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_rad_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_rad_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_rad_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_plantcarbon_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_plantcarbon_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_plantcarbon_real64 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_soilcarbon_int32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_soilcarbon_real32 + class(cable_netcdf_decomp_t), pointer :: output_decomp_base_patch_soilcarbon_real64 + + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_int32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_real32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_real64 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_soil_int32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_soil_real32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_soil_real64 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_snow_int32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_snow_real32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_snow_real64 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_rad_int32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_rad_real32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_rad_real64 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_plantcarbon_int32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_plantcarbon_real32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_plantcarbon_real64 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_soilcarbon_int32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_soilcarbon_real32 + class(cable_netcdf_decomp_t), pointer :: restart_decomp_patch_soilcarbon_real64 + + ! Temporary buffers for computing grid-cell averages for each variable class + integer(kind=int32), allocatable, target :: temp_buffer_land_int32(:) + real(kind=real32), allocatable, target :: temp_buffer_land_real32(:) + real(kind=real64), allocatable, target :: temp_buffer_land_real64(:) + integer(kind=int32), allocatable, target :: temp_buffer_land_soil_int32(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_soil_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_soil_real64(:, :) + integer(kind=int32), allocatable, target :: temp_buffer_land_snow_int32(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_snow_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_snow_real64(:, :) + integer(kind=int32), allocatable, target :: temp_buffer_land_rad_int32(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_rad_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_rad_real64(:, :) + integer(kind=int32), allocatable, target :: temp_buffer_land_plantcarbon_int32(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_plantcarbon_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_plantcarbon_real64(:, :) + integer(kind=int32), allocatable, target :: temp_buffer_land_soilcarbon_int32(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_soilcarbon_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_soilcarbon_real64(:, :) + +contains + + subroutine init_decomp_pointers(io_decomp) + type(io_decomp_t), intent(in), target :: io_decomp + + if (requires_x_y_output_grid(output%grid, metGrid)) then + output_decomp_base_int32 => io_decomp%land_to_x_y_int32 + output_decomp_base_real32 => io_decomp%land_to_x_y_real32 + output_decomp_base_real64 => io_decomp%land_to_x_y_real64 + output_decomp_base_soil_int32 => io_decomp%land_soil_to_x_y_soil_int32 + output_decomp_base_soil_real32 => io_decomp%land_soil_to_x_y_soil_real32 + output_decomp_base_soil_real64 => io_decomp%land_soil_to_x_y_soil_real64 + output_decomp_base_snow_int32 => io_decomp%land_snow_to_x_y_snow_int32 + output_decomp_base_snow_real32 => io_decomp%land_snow_to_x_y_snow_real32 + output_decomp_base_snow_real64 => io_decomp%land_snow_to_x_y_snow_real64 + output_decomp_base_rad_int32 => io_decomp%land_rad_to_x_y_rad_int32 + output_decomp_base_rad_real32 => io_decomp%land_rad_to_x_y_rad_real32 + output_decomp_base_rad_real64 => io_decomp%land_rad_to_x_y_rad_real64 + output_decomp_base_plantcarbon_int32 => io_decomp%land_plantcarbon_to_x_y_plantcarbon_int32 + output_decomp_base_plantcarbon_real32 => io_decomp%land_plantcarbon_to_x_y_plantcarbon_real32 + output_decomp_base_plantcarbon_real64 => io_decomp%land_plantcarbon_to_x_y_plantcarbon_real64 + output_decomp_base_soilcarbon_int32 => io_decomp%land_soilcarbon_to_x_y_soilcarbon_int32 + output_decomp_base_soilcarbon_real32 => io_decomp%land_soilcarbon_to_x_y_soilcarbon_real32 + output_decomp_base_soilcarbon_real64 => io_decomp%land_soilcarbon_to_x_y_soilcarbon_real64 + output_decomp_base_patch_int32 => io_decomp%patch_to_x_y_patch_int32 + output_decomp_base_patch_real32 => io_decomp%patch_to_x_y_patch_real32 + output_decomp_base_patch_real64 => io_decomp%patch_to_x_y_patch_real64 + output_decomp_base_patch_soil_int32 => io_decomp%patch_soil_to_x_y_patch_soil_int32 + output_decomp_base_patch_soil_real32 => io_decomp%patch_soil_to_x_y_patch_soil_real32 + output_decomp_base_patch_soil_real64 => io_decomp%patch_soil_to_x_y_patch_soil_real64 + output_decomp_base_patch_snow_int32 => io_decomp%patch_snow_to_x_y_patch_snow_int32 + output_decomp_base_patch_snow_real32 => io_decomp%patch_snow_to_x_y_patch_snow_real32 + output_decomp_base_patch_snow_real64 => io_decomp%patch_snow_to_x_y_patch_snow_real64 + output_decomp_base_patch_rad_int32 => io_decomp%patch_rad_to_x_y_patch_rad_int32 + output_decomp_base_patch_rad_real32 => io_decomp%patch_rad_to_x_y_patch_rad_real32 + output_decomp_base_patch_rad_real64 => io_decomp%patch_rad_to_x_y_patch_rad_real64 + output_decomp_base_patch_plantcarbon_int32 => io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_int32 + output_decomp_base_patch_plantcarbon_real32 => io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_real32 + output_decomp_base_patch_plantcarbon_real64 => io_decomp%patch_plantcarbon_to_x_y_patch_plantcarbon_real64 + output_decomp_base_patch_soilcarbon_int32 => io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_int32 + output_decomp_base_patch_soilcarbon_real32 => io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_real32 + output_decomp_base_patch_soilcarbon_real64 => io_decomp%patch_soilcarbon_to_x_y_patch_soilcarbon_real64 + else if (requires_land_output_grid(output%grid, metGrid)) then + output_decomp_base_int32 => io_decomp%land_to_land_int32 + output_decomp_base_real32 => io_decomp%land_to_land_real32 + output_decomp_base_real64 => io_decomp%land_to_land_real64 + output_decomp_base_soil_int32 => io_decomp%land_soil_to_land_soil_int32 + output_decomp_base_soil_real32 => io_decomp%land_soil_to_land_soil_real32 + output_decomp_base_soil_real64 => io_decomp%land_soil_to_land_soil_real64 + output_decomp_base_snow_int32 => io_decomp%land_snow_to_land_snow_int32 + output_decomp_base_snow_real32 => io_decomp%land_snow_to_land_snow_real32 + output_decomp_base_snow_real64 => io_decomp%land_snow_to_land_snow_real64 + output_decomp_base_rad_int32 => io_decomp%land_rad_to_land_rad_int32 + output_decomp_base_rad_real32 => io_decomp%land_rad_to_land_rad_real32 + output_decomp_base_rad_real64 => io_decomp%land_rad_to_land_rad_real64 + output_decomp_base_plantcarbon_int32 => io_decomp%land_plantcarbon_to_land_plantcarbon_int32 + output_decomp_base_plantcarbon_real32 => io_decomp%land_plantcarbon_to_land_plantcarbon_real32 + output_decomp_base_plantcarbon_real64 => io_decomp%land_plantcarbon_to_land_plantcarbon_real64 + output_decomp_base_soilcarbon_int32 => io_decomp%land_soilcarbon_to_land_soilcarbon_int32 + output_decomp_base_soilcarbon_real32 => io_decomp%land_soilcarbon_to_land_soilcarbon_real32 + output_decomp_base_soilcarbon_real64 => io_decomp%land_soilcarbon_to_land_soilcarbon_real64 + output_decomp_base_patch_int32 => io_decomp%patch_to_land_patch_int32 + output_decomp_base_patch_real32 => io_decomp%patch_to_land_patch_real32 + output_decomp_base_patch_real64 => io_decomp%patch_to_land_patch_real64 + output_decomp_base_patch_soil_int32 => io_decomp%patch_soil_to_land_patch_soil_int32 + output_decomp_base_patch_soil_real32 => io_decomp%patch_soil_to_land_patch_soil_real32 + output_decomp_base_patch_soil_real64 => io_decomp%patch_soil_to_land_patch_soil_real64 + output_decomp_base_patch_snow_int32 => io_decomp%patch_snow_to_land_patch_snow_int32 + output_decomp_base_patch_snow_real32 => io_decomp%patch_snow_to_land_patch_snow_real32 + output_decomp_base_patch_snow_real64 => io_decomp%patch_snow_to_land_patch_snow_real64 + output_decomp_base_patch_rad_int32 => io_decomp%patch_rad_to_land_patch_rad_int32 + output_decomp_base_patch_rad_real32 => io_decomp%patch_rad_to_land_patch_rad_real32 + output_decomp_base_patch_rad_real64 => io_decomp%patch_rad_to_land_patch_rad_real64 + output_decomp_base_patch_plantcarbon_int32 => io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_int32 + output_decomp_base_patch_plantcarbon_real32 => io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_real32 + output_decomp_base_patch_plantcarbon_real64 => io_decomp%patch_plantcarbon_to_land_patch_plantcarbon_real64 + output_decomp_base_patch_soilcarbon_int32 => io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_int32 + output_decomp_base_patch_soilcarbon_real32 => io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_real32 + output_decomp_base_patch_soilcarbon_real64 => io_decomp%patch_soilcarbon_to_land_patch_soilcarbon_real64 + else + call cable_abort("Error: Unable to determine output grid type", __FILE__, __LINE__) + end if + + restart_decomp_patch_int32 => io_decomp%patch_to_patch_int32 + restart_decomp_patch_real32 => io_decomp%patch_to_patch_real32 + restart_decomp_patch_real64 => io_decomp%patch_to_patch_real64 + restart_decomp_patch_soil_int32 => io_decomp%patch_soil_to_patch_soil_int32 + restart_decomp_patch_soil_real32 => io_decomp%patch_soil_to_patch_soil_real32 + restart_decomp_patch_soil_real64 => io_decomp%patch_soil_to_patch_soil_real64 + restart_decomp_patch_snow_int32 => io_decomp%patch_snow_to_patch_snow_int32 + restart_decomp_patch_snow_real32 => io_decomp%patch_snow_to_patch_snow_real32 + restart_decomp_patch_snow_real64 => io_decomp%patch_snow_to_patch_snow_real64 + restart_decomp_patch_rad_int32 => io_decomp%patch_rad_to_patch_rad_int32 + restart_decomp_patch_rad_real32 => io_decomp%patch_rad_to_patch_rad_real32 + restart_decomp_patch_rad_real64 => io_decomp%patch_rad_to_patch_rad_real64 + restart_decomp_patch_plantcarbon_int32 => io_decomp%patch_plantcarbon_to_patch_plantcarbon_int32 + restart_decomp_patch_plantcarbon_real32 => io_decomp%patch_plantcarbon_to_patch_plantcarbon_real32 + restart_decomp_patch_plantcarbon_real64 => io_decomp%patch_plantcarbon_to_patch_plantcarbon_real64 + restart_decomp_patch_soilcarbon_int32 => io_decomp%patch_soilcarbon_to_patch_soilcarbon_int32 + restart_decomp_patch_soilcarbon_real32 => io_decomp%patch_soilcarbon_to_patch_soilcarbon_real32 + restart_decomp_patch_soilcarbon_real64 => io_decomp%patch_soilcarbon_to_patch_soilcarbon_real64 + + end subroutine + + subroutine allocate_grid_reduction_buffers() + + allocate(temp_buffer_land_int32(mland)) + allocate(temp_buffer_land_real32(mland)) + allocate(temp_buffer_land_real64(mland)) + allocate(temp_buffer_land_soil_int32(mland, ms)) + allocate(temp_buffer_land_soil_real32(mland, ms)) + allocate(temp_buffer_land_soil_real64(mland, ms)) + allocate(temp_buffer_land_snow_int32(mland, msn)) + allocate(temp_buffer_land_snow_real32(mland, msn)) + allocate(temp_buffer_land_snow_real64(mland, msn)) + allocate(temp_buffer_land_rad_int32(mland, nrb)) + allocate(temp_buffer_land_rad_real32(mland, nrb)) + allocate(temp_buffer_land_rad_real64(mland, nrb)) + allocate(temp_buffer_land_plantcarbon_int32(mland, ncp)) + allocate(temp_buffer_land_plantcarbon_real32(mland, ncp)) + allocate(temp_buffer_land_plantcarbon_real64(mland, ncp)) + allocate(temp_buffer_land_soilcarbon_int32(mland, ncs)) + allocate(temp_buffer_land_soilcarbon_real32(mland, ncs)) + allocate(temp_buffer_land_soilcarbon_real64(mland, ncs)) + + end subroutine + + subroutine deallocate_grid_reduction_buffers() + + deallocate(temp_buffer_land_int32) + deallocate(temp_buffer_land_real32) + deallocate(temp_buffer_land_real64) + deallocate(temp_buffer_land_soil_int32) + deallocate(temp_buffer_land_soil_real32) + deallocate(temp_buffer_land_soil_real64) + deallocate(temp_buffer_land_snow_int32) + deallocate(temp_buffer_land_snow_real32) + deallocate(temp_buffer_land_snow_real64) + deallocate(temp_buffer_land_rad_int32) + deallocate(temp_buffer_land_rad_real32) + deallocate(temp_buffer_land_rad_real64) + deallocate(temp_buffer_land_plantcarbon_int32) + deallocate(temp_buffer_land_plantcarbon_real32) + deallocate(temp_buffer_land_plantcarbon_real64) + deallocate(temp_buffer_land_soilcarbon_int32) + deallocate(temp_buffer_land_soilcarbon_real32) + deallocate(temp_buffer_land_soilcarbon_real64) + + end subroutine + + logical function requires_x_y_output_grid(output_grid, met_grid) + character(len=*), intent(in) :: output_grid + character(len=*), intent(in) :: met_grid + requires_x_y_output_grid = (( & + output_grid == "default" .AND. met_grid == "mask" & + ) .OR. ( & + output_grid == "mask" .OR. output_grid == "ALMA" & + )) + end function + + logical function requires_land_output_grid(output_grid, met_grid) + character(len=*), intent(in) :: output_grid + character(len=*), intent(in) :: met_grid + requires_land_output_grid = ( & + output_grid == "land" .OR. (output_grid == "default" .AND. met_grid == "land") & + ) + end function + + logical function data_shape_eq(shape1, shape2) + type(cable_output_dim_t), dimension(:), intent(in) :: shape1, shape2 + data_shape_eq = size(shape1) == size(shape2) .and. all(shape1 == shape2) + end function + + elemental integer function dim_size(dim) + type(cable_output_dim_t), intent(in) :: dim + + select case (dim%value) + case (CABLE_OUTPUT_DIM_PATCH%value) + dim_size = mp + case (CABLE_OUTPUT_DIM_SOIL%value) + dim_size = ms + case (CABLE_OUTPUT_DIM_SNOW%value) + dim_size = msn + case (CABLE_OUTPUT_DIM_RAD%value) + dim_size = nrb + case (CABLE_OUTPUT_DIM_PLANTCARBON%value) + dim_size = ncp + case (CABLE_OUTPUT_DIM_SOILCARBON%value) + dim_size = ncs + case (CABLE_OUTPUT_DIM_LAND%value) + dim_size = mland + case (CABLE_OUTPUT_DIM_LAND_GLOBAL%value) + dim_size = mland_global + case (CABLE_OUTPUT_DIM_X%value) + dim_size = xdimsize + case (CABLE_OUTPUT_DIM_Y%value) + dim_size = ydimsize + case default + dim_size = -1 ! Unknown dimension + end select + + end function + + subroutine check_invalid_frequency(sampling_frequency, accumulation_frequency, var_name, file_name) + character(len=*), intent(in) :: sampling_frequency + character(len=*), intent(in) :: accumulation_frequency + character(len=*), intent(in) :: var_name + character(len=*), intent(in) :: file_name + + integer :: sampling_period_in_hours, accumulation_period_in_hours + + character(len=256) :: err_message + + err_message = ( & + "Invalid combination of sampling frequency '" // sampling_frequency // & + "' with accumulation frequency '" // accumulation_frequency // "' for variable '" // & + var_name // "' in file '" // file_name // "'" & + ) + + select case (sampling_frequency) + case ("all") + if (accumulation_frequency /= "all") call cable_abort(err_message, __FILE__, __LINE__) + case ("user") + read(sampling_frequency(5:7), *) sampling_period_in_hours + if (accumulation_frequency == "user") then + read(accumulation_frequency(5:7), *) accumulation_period_in_hours + if (sampling_period_in_hours < accumulation_period_in_hours) then + call cable_abort(err_message, __FILE__, __LINE__) + end if + else if (accumulation_frequency /= "all") then + call cable_abort(err_message, __FILE__, __LINE__) + end if + case ("daily") + if (.not. any(accumulation_frequency == ["all ", "daily", "user "])) then + call cable_abort(err_message, __FILE__, __LINE__) + end if + case ("monthly") + if (.not. any(accumulation_frequency == ["all ", "daily ", "user ", "monthly"])) then + call cable_abort(err_message, __FILE__, __LINE__) + end if + case default + call cable_abort("Invalid sampling frequency '" // sampling_frequency // & + "' for variable '" // var_name // "' in file '" // file_name // "'", __FILE__, __LINE__) + end select + + end subroutine check_invalid_frequency + + function infer_dim_names(output_variable, restart) result(dim_names) + type(cable_output_variable_t), intent(in) :: output_variable + logical, intent(in), optional :: restart + + character(MAX_LEN_DIM), allocatable :: dim_names(:) + logical :: restart_local + integer :: j + + restart_local = .false. + if (present(restart)) restart_local = restart + + allocate(dim_names(0)) + if (allocated(output_variable%data_shape)) then + do j = 1, size(output_variable%data_shape) + select case (output_variable%data_shape(j)%value) + case (CABLE_OUTPUT_DIM_PATCH%value) + if (restart_local) then + dim_names = [dim_names, "mp"] + else if (requires_land_output_grid(output%grid, metgrid)) then + if (output_variable%reduction_method == "none") then + dim_names = [dim_names, "land", "patch"] + else + dim_names = [dim_names, "land"] + end if + else if (requires_x_y_output_grid(output%grid, metgrid)) then + if (output_variable%reduction_method == "none") then + dim_names = [dim_names, "x", "y", "patch"] + else + dim_names = [dim_names, "x", "y"] + end if + end if + case (CABLE_OUTPUT_DIM_LAND%value) + if (restart_local) then + dim_names = [dim_names, "mland"] + else if (requires_land_output_grid(output%grid, metgrid)) then + dim_names = [dim_names, "land"] + else if (requires_x_y_output_grid(output%grid, metgrid)) then + dim_names = [dim_names, "x", "y"] + end if + case (CABLE_OUTPUT_DIM_LAND_GLOBAL%value) + if (restart_local) then + dim_names = [dim_names, "mland"] + else + dim_names = [dim_names, "land"] + end if + case (CABLE_OUTPUT_DIM_SOIL%value) + dim_names = [dim_names, "soil"] + case (CABLE_OUTPUT_DIM_SNOW%value) + dim_names = [dim_names, "snow"] + case (CABLE_OUTPUT_DIM_RAD%value) + dim_names = [dim_names, "rad"] + case (CABLE_OUTPUT_DIM_PLANTCARBON%value) + if (restart_local) then + dim_names = [dim_names, "plant_carbon_pools"] + else + dim_names = [dim_names, "plantcarbon"] + end if + case (CABLE_OUTPUT_DIM_SOILCARBON%value) + if (restart_local) then + dim_names = [dim_names, "soil_carbon_pools"] + else + dim_names = [dim_names, "soilcarbon"] + end if + case (CABLE_OUTPUT_DIM_X%value) + dim_names = [dim_names, "x"] + case (CABLE_OUTPUT_DIM_Y%value) + dim_names = [dim_names, "y"] + case default + call cable_abort("Unexpected data shape for variable " // output_variable%name, __FILE__, __LINE__) + end select + end do + end if + + if (.not. restart_local .and. .not. output_variable%parameter) dim_names = [dim_names, "time"] + + end function + + function infer_cell_methods(output_variable) result(cell_methods) + type(cable_output_variable_t), intent(in) :: output_variable + character(len=256) :: cell_methods + character(len=256) :: cell_methods_time, cell_methods_area + + if (.not. output_variable%parameter) then + select case (output_variable%aggregation_method) + case ("point") + cell_methods_time = "time: point" + case ("mean") + cell_methods_time = "time: mean" + case ("sum") + cell_methods_time = "time: sum" + case ("min") + cell_methods_time = "time: minimum" + case ("max") + cell_methods_time = "time: maximum" + case default + call cable_abort("Unexpected aggregation method '" // output_variable%aggregation_method // & + "' for variable '" // output_variable%name // "'", __FILE__, __LINE__) + end select + else + cell_methods_time = "" + end if + + select case (output_variable%reduction_method) + case ("none") + ! TODO(Sean): the cell_method for this case should be `area: point where + ! pft` where `pft` is a string-valued auxiliary coordinate variable + ! describing the labels for all patches: + cell_methods_area = "" + case ("first_patch_in_grid_cell") + ! TODO(Sean): the cell_method for this case should be `area: point where + ! ` where `` is the area type of the first patch in + ! the grid cell: + cell_methods_area = "" + case ("grid_cell_average") + cell_methods_area = "area: mean where land" + case default + call cable_abort("Unexpected reduction method '" // output_variable%reduction_method // & + "' for variable '" // output_variable%name // "'", __FILE__, __LINE__) + end select + + cell_methods = adjustl(trim(cell_methods_time) // " " // trim(cell_methods_area)) + + end function + + subroutine define_variables(output_file, output_variables, restart) + class(cable_netcdf_file_t), intent(inout) :: output_file + type(cable_output_variable_t), intent(in) :: output_variables(:) + logical, intent(in), optional :: restart + + integer :: i, j + logical :: restart_local + + character(MAX_LEN_DIM), allocatable :: required_dimensions(:), dim_names(:) + + restart_local = .false. + if (present(restart)) restart_local = restart + + do i = 1, size(output_variables) + associate(output_var => output_variables(i)) + if (restart_local .and. .not. output_var%restart) cycle + if (.not. allocated(output_var%data_shape)) cycle + dim_names = infer_dim_names(output_var, restart_local) + if (.not. allocated(required_dimensions)) then + required_dimensions = dim_names + else + required_dimensions = [ & + required_dimensions, & + pack(dim_names, [( & + .not. any(dim_names(j) == required_dimensions), & + j = 1, & + size(dim_names) & + )]) & + ] + end if + end associate + end do + + do i = 1, size(required_dimensions) + select case (required_dimensions(i)) + case ("mp") + call output_file%def_dims(["mp"], [mp_global]) + case ("mland") + call output_file%def_dims(["mland"], [mland_global]) + case ("land") + call output_file%def_dims(["land"], [mland_global]) + case ("x") + call output_file%def_dims(["x"], [xdimsize]) + case ("y") + call output_file%def_dims(["y"], [ydimsize]) + case ("patch") + call output_file%def_dims(["patch"], [max_vegpatches]) + case ("soil") + call output_file%def_dims(["soil"], [ms]) + case ("rad") + call output_file%def_dims(["rad"], [nrb]) + case ("soil_carbon_pools") + call output_file%def_dims(["soil_carbon_pools"], [ncs]) + case ("soilcarbon") + call output_file%def_dims(["soilcarbon"], [ncs]) + case ("plant_carbon_pools") + call output_file%def_dims(["plant_carbon_pools"], [ncp]) + case ("plantcarbon") + call output_file%def_dims(["plantcarbon"], [ncp]) + case ("time") + ! time dimension defined separately below + case default + call cable_abort("Unexpected dimension name '" // required_dimensions(i) // "'", __FILE__, __LINE__) + end select + end do + + if (restart_local) then + call output_file%def_dims(["time"], [1]) + else + call output_file%def_dims(["time"], [CABLE_NETCDF_UNLIMITED]) + end if + + call output_file%def_var("time", ["time"], CABLE_NETCDF_DOUBLE) + call output_file%put_att("time", "units", timeunits) + call output_file%put_att("time", "coordinate", time_coord) + call output_file%put_att("time", "calendar", calendar) + + if (.not. restart_local) then + call output_file%def_dims(["nv"], [2]) + call output_file%def_var("time_bnds", ["nv ", "time"], CABLE_NETCDF_DOUBLE) + call output_file%put_att("time", "bounds", "time_bnds") + end if + + do i = 1, size(output_variables) + associate(output_var => output_variables(i)) + if (restart_local .and. .not. output_var%restart) cycle + call output_file%def_var( & + var_name=output_var%name, & + dim_names=infer_dim_names(output_var, restart_local), & + type=output_var%var_type & + ) + if (allocated(output_var%metadata)) then + do j = 1, size(output_var%metadata) + call output_file%put_att(output_var%name, output_var%metadata(j)%name, output_var%metadata(j)%value) + end do + end if + select case (output_var%var_type) + case (CABLE_NETCDF_INT) + call output_file%put_att(output_var%name, "_FillValue", FILL_VALUE_INT32) + call output_file%put_att(output_var%name, "missing_value", FILL_VALUE_INT32) + case (CABLE_NETCDF_FLOAT) + call output_file%put_att(output_var%name, "_FillValue", FILL_VALUE_REAL32) + call output_file%put_att(output_var%name, "missing_value", FILL_VALUE_REAL32) + case (CABLE_NETCDF_DOUBLE) + call output_file%put_att(output_var%name, "_FillValue", FILL_VALUE_REAL64) + call output_file%put_att(output_var%name, "missing_value", FILL_VALUE_REAL64) + end select + call output_file%put_att(output_var%name, "cell_methods", infer_cell_methods(output_var)) + end associate + end do + + end subroutine define_variables + + subroutine set_global_attributes(output_profile) + type(cable_output_profile_t), intent(inout) :: output_profile + + character(32) :: todaydate, nowtime + integer :: i + + if (allocated(output_profile%metadata)) then + do i = 1, size(output_profile%metadata) + call output_profile%output_file%put_att( & + att_name=output_profile%metadata(i)%name, & + att_value=output_profile%metadata(i)%value & + ) + end do + end if + + call date_and_time(todaydate, nowtime) + todaydate = todaydate(1:4) // "/" // todaydate(5:6) // "/" // todaydate(7:8) + nowtime = nowtime(1:2) // ":" // nowtime(3:4) // ":" // nowtime(5:6) + call output_profile%output_file%put_att("Production", trim(todaydate) // " at " // trim(nowtime)) + call output_profile%output_file%put_att("Source", "CABLE LSM output file") + call output_profile%output_file%put_att("CABLE_input_file", trim(filename%met)) + + select case (output_profile%sampling_frequency) + case ("user") + call output_profile%output_file%put_att("Output_averaging", TRIM(output_profile%sampling_frequency(5:7)) // "-hourly output") + case ("all") + call output_profile%output_file%put_att("Output_averaging", "all timesteps recorded") + case ("daily") + call output_profile%output_file%put_att("Output_averaging", "daily") + case ("monthly") + call output_profile%output_file%put_att("Output_averaging", "monthly") + case default + call cable_abort("Invalid sampling frequency '" // output_profile%sampling_frequency // "'", __FILE__, __LINE__) + end select + + end subroutine set_global_attributes + + subroutine associate_decomp_int32(output_var, decomp, restart) + type(cable_output_variable_t), intent(in) :: output_var + class(cable_netcdf_decomp_t), pointer, intent(inout) :: decomp + logical, intent(in), optional :: restart + + if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_int32 + else + decomp => output_decomp_base_int32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_int32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOIL])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_soil_int32 + else + decomp => output_decomp_base_soil_int32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_soil_int32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SNOW])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_snow_int32 + else + decomp => output_decomp_base_snow_int32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_snow_int32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_rad_int32 + else + decomp => output_decomp_base_rad_int32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_rad_int32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_PLANTCARBON])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_plantcarbon_int32 + else + decomp => output_decomp_base_plantcarbon_int32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_plantcarbon_int32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOILCARBON])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_soilcarbon_int32 + else + decomp => output_decomp_base_soilcarbon_int32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_soilcarbon_int32 + end if + else + call cable_abort("Unsupported data shape for output variable " // output_var%name, __FILE__, __LINE__) + end if + + end subroutine associate_decomp_int32 + + subroutine associate_decomp_real32(output_var, decomp, restart) + type(cable_output_variable_t), intent(in) :: output_var + class(cable_netcdf_decomp_t), pointer, intent(inout) :: decomp + logical, intent(in), optional :: restart + + if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_real32 + else + decomp => output_decomp_base_real32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_real32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOIL])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_soil_real32 + else + decomp => output_decomp_base_soil_real32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_soil_real32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SNOW])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_snow_real32 + else + decomp => output_decomp_base_snow_real32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_snow_real32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_rad_real32 + else + decomp => output_decomp_base_rad_real32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_rad_real32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_PLANTCARBON])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_plantcarbon_real32 + else + decomp => output_decomp_base_plantcarbon_real32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_plantcarbon_real32 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOILCARBON])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_soilcarbon_real32 + else + decomp => output_decomp_base_soilcarbon_real32 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_soilcarbon_real32 + end if + else + call cable_abort("Unsupported data shape for output variable " // output_var%name, __FILE__, __LINE__) + end if + + end subroutine associate_decomp_real32 + + subroutine associate_decomp_real64(output_var, decomp, restart) + type(cable_output_variable_t), intent(in) :: output_var + class(cable_netcdf_decomp_t), pointer, intent(inout) :: decomp + logical, intent(in), optional :: restart + + if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_real64 + else + decomp => output_decomp_base_real64 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_real64 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOIL])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_soil_real64 + else + decomp => output_decomp_base_soil_real64 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_soil_real64 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SNOW])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_snow_real64 + else + decomp => output_decomp_base_snow_real64 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_snow_real64 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_rad_real64 + else + decomp => output_decomp_base_rad_real64 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_rad_real64 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_PLANTCARBON])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_plantcarbon_real64 + else + decomp => output_decomp_base_plantcarbon_real64 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_plantcarbon_real64 + end if + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOILCARBON])) then + if (output_var%reduction_method == "none") then + decomp => output_decomp_base_patch_soilcarbon_real64 + else + decomp => output_decomp_base_soilcarbon_real64 + end if + if (present(restart)) then + if (restart) decomp => restart_decomp_patch_soilcarbon_real64 + end if + else + call cable_abort("Unsupported data shape for output variable " // output_var%name, __FILE__, __LINE__) + end if + + end subroutine associate_decomp_real64 + + subroutine associate_temp_buffer_int32(output_var, temp_buffer_int32_1d, temp_buffer_int32_2d, temp_buffer_int32_3d) + type(cable_output_variable_t), intent(inout) :: output_var + integer(kind=int32), pointer, intent(inout), optional :: temp_buffer_int32_1d(:) + integer(kind=int32), pointer, intent(inout), optional :: temp_buffer_int32_2d(:,:) + integer(kind=int32), pointer, intent(inout), optional :: temp_buffer_int32_3d(:,:,:) + + if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH])) then + if (.not. present(temp_buffer_int32_1d)) call cable_abort( & + "temp_buffer_int32_1d must be provided for 1D data shape", __FILE__, __LINE__) + temp_buffer_int32_1d => temp_buffer_land_int32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOIL])) then + if (.not. present(temp_buffer_int32_2d)) call cable_abort( & + "temp_buffer_int32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_int32_2d => temp_buffer_land_soil_int32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (.not. present(temp_buffer_int32_2d)) call cable_abort( & + "temp_buffer_int32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_int32_2d => temp_buffer_land_rad_int32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SNOW])) then + if (.not. present(temp_buffer_int32_2d)) call cable_abort( & + "temp_buffer_int32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_int32_2d => temp_buffer_land_snow_int32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (.not. present(temp_buffer_int32_2d)) call cable_abort( & + "temp_buffer_int32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_int32_2d => temp_buffer_land_rad_int32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_PLANTCARBON])) then + if (.not. present(temp_buffer_int32_2d)) call cable_abort( & + "temp_buffer_int32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_int32_2d => temp_buffer_land_plantcarbon_int32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOILCARBON])) then + if (.not. present(temp_buffer_int32_2d)) call cable_abort( & + "temp_buffer_int32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_int32_2d => temp_buffer_land_soilcarbon_int32 + else + call cable_abort("Unexpected source data shape for grid reduction", __FILE__, __LINE__) + end if + + end subroutine associate_temp_buffer_int32 + + subroutine associate_temp_buffer_real32(output_var, temp_buffer_real32_1d, temp_buffer_real32_2d, temp_buffer_real32_3d) + type(cable_output_variable_t), intent(inout) :: output_var + real(kind=real32), pointer, intent(inout), optional :: temp_buffer_real32_1d(:) + real(kind=real32), pointer, intent(inout), optional :: temp_buffer_real32_2d(:,:) + real(kind=real32), pointer, intent(inout), optional :: temp_buffer_real32_3d(:,:,:) + + if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH])) then + if (.not. present(temp_buffer_real32_1d)) call cable_abort( & + "temp_buffer_real32_1d must be provided for 1D data shape", __FILE__, __LINE__) + temp_buffer_real32_1d => temp_buffer_land_real32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOIL])) then + if (.not. present(temp_buffer_real32_2d)) call cable_abort( & + "temp_buffer_real32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real32_2d => temp_buffer_land_soil_real32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (.not. present(temp_buffer_real32_2d)) call cable_abort( & + "temp_buffer_real32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real32_2d => temp_buffer_land_rad_real32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SNOW])) then + if (.not. present(temp_buffer_real32_2d)) call cable_abort( & + "temp_buffer_real32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real32_2d => temp_buffer_land_snow_real32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (.not. present(temp_buffer_real32_2d)) call cable_abort( & + "temp_buffer_real32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real32_2d => temp_buffer_land_rad_real32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_PLANTCARBON])) then + if (.not. present(temp_buffer_real32_2d)) call cable_abort( & + "temp_buffer_real32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real32_2d => temp_buffer_land_plantcarbon_real32 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOILCARBON])) then + if (.not. present(temp_buffer_real32_2d)) call cable_abort( & + "temp_buffer_real32_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real32_2d => temp_buffer_land_soilcarbon_real32 + else + call cable_abort("Unexpected source data shape for grid reduction", __FILE__, __LINE__) + end if + + end subroutine associate_temp_buffer_real32 + + subroutine associate_temp_buffer_real64(output_var, temp_buffer_real64_1d, temp_buffer_real64_2d, temp_buffer_real64_3d) + type(cable_output_variable_t), intent(inout) :: output_var + real(kind=real64), pointer, intent(inout), optional :: temp_buffer_real64_1d(:) + real(kind=real64), pointer, intent(inout), optional :: temp_buffer_real64_2d(:,:) + real(kind=real64), pointer, intent(inout), optional :: temp_buffer_real64_3d(:,:,:) + + if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH])) then + if (.not. present(temp_buffer_real64_1d)) call cable_abort( & + "temp_buffer_real64_1d must be provided for 1D data shape", __FILE__, __LINE__) + temp_buffer_real64_1d => temp_buffer_land_real64 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOIL])) then + if (.not. present(temp_buffer_real64_2d)) call cable_abort( & + "temp_buffer_real64_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real64_2d => temp_buffer_land_soil_real64 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (.not. present(temp_buffer_real64_2d)) call cable_abort( & + "temp_buffer_real64_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real64_2d => temp_buffer_land_rad_real64 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SNOW])) then + if (.not. present(temp_buffer_real64_2d)) call cable_abort( & + "temp_buffer_real64_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real64_2d => temp_buffer_land_snow_real64 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_RAD])) then + if (.not. present(temp_buffer_real64_2d)) call cable_abort( & + "temp_buffer_real64_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real64_2d => temp_buffer_land_rad_real64 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_PLANTCARBON])) then + if (.not. present(temp_buffer_real64_2d)) call cable_abort( & + "temp_buffer_real64_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real64_2d => temp_buffer_land_plantcarbon_real64 + else if (data_shape_eq(output_var%data_shape, [CABLE_OUTPUT_DIM_PATCH, CABLE_OUTPUT_DIM_SOILCARBON])) then + if (.not. present(temp_buffer_real64_2d)) call cable_abort( & + "temp_buffer_real64_2d must be provided for 2D data shape", __FILE__, __LINE__) + temp_buffer_real64_2d => temp_buffer_land_soilcarbon_real64 + else + call cable_abort("Unexpected source data shape for grid reduction", __FILE__, __LINE__) + end if + + end subroutine associate_temp_buffer_real64 + +end module diff --git a/src/util/netcdf/cable_netcdf.F90 b/src/util/netcdf/cable_netcdf.F90 index 086c9c307..500368249 100644 --- a/src/util/netcdf/cable_netcdf.F90 +++ b/src/util/netcdf/cable_netcdf.F90 @@ -30,7 +30,14 @@ module cable_netcdf_mod CABLE_NETCDF_MAX_STR_LEN_VAR, & CABLE_NETCDF_MAX_STR_LEN_DIM, & CABLE_NETCDF_MAX_RANK, & - CABLE_NETCDF_UNLIMITED + CABLE_NETCDF_UNLIMITED, & + CABLE_NETCDF_IOTYPE_CLASSIC, & + CABLE_NETCDF_IOTYPE_NETCDF4C, & + CABLE_NETCDF_IOTYPE_NETCDF4P, & + CABLE_NETCDF_MODE_CLOBBER, & + CABLE_NETCDF_MODE_NOCLOBBER, & + CABLE_NETCDF_MODE_WRITE, & + CABLE_NETCDF_MODE_NOWRITE enum, bind(c) enumerator :: & @@ -39,6 +46,21 @@ module cable_netcdf_mod CABLE_NETCDF_DOUBLE end enum + enum, bind(c) + enumerator :: & + CABLE_NETCDF_IOTYPE_CLASSIC, & + CABLE_NETCDF_IOTYPE_NETCDF4C, & + CABLE_NETCDF_IOTYPE_NETCDF4P + end enum + + enum, bind(c) + enumerator :: & + CABLE_NETCDF_MODE_CLOBBER, & + CABLE_NETCDF_MODE_NOCLOBBER, & + CABLE_NETCDF_MODE_WRITE, & + CABLE_NETCDF_MODE_NOWRITE + end enum + integer, parameter :: CABLE_NETCDF_MAX_STR_LEN_FILE = 200 integer, parameter :: CABLE_NETCDF_MAX_STR_LEN_VAR = 80 integer, parameter :: CABLE_NETCDF_MAX_STR_LEN_DIM = 20 @@ -55,6 +77,7 @@ module cable_netcdf_mod contains procedure(cable_netcdf_file_close), deferred :: close procedure(cable_netcdf_file_end_def), deferred :: end_def + procedure(cable_netcdf_file_redef), deferred :: redef procedure(cable_netcdf_file_sync), deferred :: sync procedure(cable_netcdf_file_def_dims), deferred :: def_dims procedure(cable_netcdf_file_def_var), deferred :: def_var @@ -81,6 +104,7 @@ module cable_netcdf_mod get_att_global_string, get_att_global_int32, get_att_global_real32, get_att_global_real64, & get_att_var_string, get_att_var_int32, get_att_var_real32, get_att_var_real64 procedure(cable_netcdf_file_inq_dim_len), deferred :: inq_dim_len + procedure(cable_netcdf_file_inq_var_ndims), deferred :: inq_var_ndims procedure(cable_netcdf_file_put_var_int32_0d), deferred :: put_var_int32_0d procedure(cable_netcdf_file_put_var_int32_1d), deferred :: put_var_int32_1d procedure(cable_netcdf_file_put_var_int32_2d), deferred :: put_var_int32_2d @@ -150,6 +174,10 @@ subroutine cable_netcdf_file_end_def(this) import cable_netcdf_file_t class(cable_netcdf_file_t), intent(inout) :: this end subroutine + subroutine cable_netcdf_file_redef(this) + import cable_netcdf_file_t + class(cable_netcdf_file_t), intent(inout) :: this + end subroutine subroutine cable_netcdf_file_sync(this) import cable_netcdf_file_t class(cable_netcdf_file_t), intent(inout) :: this @@ -163,7 +191,8 @@ subroutine cable_netcdf_file_def_dims(this, dim_names, dim_lens) subroutine cable_netcdf_file_def_var(this, var_name, dim_names, type) import cable_netcdf_file_t class(cable_netcdf_file_t), intent(inout) :: this - character(len=*), intent(in) :: var_name, dim_names(:) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: dim_names(:) integer, intent(in) :: type end subroutine subroutine cable_netcdf_file_put_att_global_string(this, att_name, att_value) @@ -266,6 +295,12 @@ subroutine cable_netcdf_file_inq_dim_len(this, dim_name, dim_len) character(len=*), intent(in) :: dim_name integer, intent(out) :: dim_len end subroutine + subroutine cable_netcdf_file_inq_var_ndims(this, var_name, ndims) + import cable_netcdf_file_t + class(cable_netcdf_file_t), intent(inout) :: this + character(len=*), intent(in) :: var_name + integer, intent(out) :: ndims + end subroutine subroutine cable_netcdf_file_put_var_int32_0d(this, var_name, values, start, count) import cable_netcdf_file_t, CABLE_NETCDF_INT32_KIND class(cable_netcdf_file_t), intent(inout) :: this @@ -607,16 +642,20 @@ subroutine cable_netcdf_io_finalise(this) import cable_netcdf_io_t class(cable_netcdf_io_t), intent(inout) :: this end subroutine - function cable_netcdf_io_create_file(this, path) result(file) + function cable_netcdf_io_create_file(this, path, iotype, mode) result(file) import cable_netcdf_io_t, cable_netcdf_file_t class(cable_netcdf_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file end function - function cable_netcdf_io_open_file(this, path) result(file) + function cable_netcdf_io_open_file(this, path, iotype, mode) result(file) import cable_netcdf_io_t, cable_netcdf_file_t class(cable_netcdf_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file end function function cable_netcdf_io_create_decomp(this, compmap, dims, type) result(decomp) @@ -634,12 +673,16 @@ module subroutine cable_netcdf_mod_init(mpi_grp) end subroutine module subroutine cable_netcdf_mod_end() end subroutine - module function cable_netcdf_create_file(path) result(file) + module function cable_netcdf_create_file(path, iotype, mode) result(file) character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file end function - module function cable_netcdf_open_file(path) result(file) + module function cable_netcdf_open_file(path, iotype, mode) result(file) character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file end function module function cable_netcdf_create_decomp(compmap, dims, type) result(decomp) diff --git a/src/util/netcdf/cable_netcdf_internal.F90 b/src/util/netcdf/cable_netcdf_internal.F90 index 1129c60bb..44616872b 100644 --- a/src/util/netcdf/cable_netcdf_internal.F90 +++ b/src/util/netcdf/cable_netcdf_internal.F90 @@ -19,16 +19,20 @@ module subroutine cable_netcdf_mod_end() call cable_netcdf_io_handler%finalise() end subroutine - module function cable_netcdf_create_file(path) result(file) + module function cable_netcdf_create_file(path, iotype, mode) result(file) character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file - file = cable_netcdf_io_handler%create_file(path) + file = cable_netcdf_io_handler%create_file(path, iotype, mode) end function - module function cable_netcdf_open_file(path) result(file) + module function cable_netcdf_open_file(path, iotype, mode) result(file) character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file - file = cable_netcdf_io_handler%open_file(path) + file = cable_netcdf_io_handler%open_file(path, iotype, mode) end function module function cable_netcdf_create_decomp(compmap, dims, type) result(decomp) diff --git a/src/util/netcdf/cable_netcdf_stub_types.F90 b/src/util/netcdf/cable_netcdf_stub_types.F90 index 6f89dc8f4..25a06b0e1 100644 --- a/src/util/netcdf/cable_netcdf_stub_types.F90 +++ b/src/util/netcdf/cable_netcdf_stub_types.F90 @@ -26,6 +26,7 @@ module cable_netcdf_stub_types_mod contains procedure :: close => cable_netcdf_stub_file_close procedure :: end_def => cable_netcdf_stub_file_end_def + procedure :: redef => cable_netcdf_stub_file_redef procedure :: sync => cable_netcdf_stub_file_sync procedure :: def_dims => cable_netcdf_stub_file_def_dims procedure :: def_var => cable_netcdf_stub_file_def_var @@ -46,6 +47,7 @@ module cable_netcdf_stub_types_mod procedure :: get_att_var_real32 => cable_netcdf_stub_file_get_att_var_real32 procedure :: get_att_var_real64 => cable_netcdf_stub_file_get_att_var_real64 procedure :: inq_dim_len => cable_netcdf_stub_file_inq_dim_len + procedure :: inq_var_ndims => cable_netcdf_stub_file_inq_var_ndims procedure :: put_var_int32_0d => cable_netcdf_stub_file_put_var_int32_0d procedure :: put_var_int32_1d => cable_netcdf_stub_file_put_var_int32_1d procedure :: put_var_int32_2d => cable_netcdf_stub_file_put_var_int32_2d @@ -100,16 +102,20 @@ subroutine cable_netcdf_stub_io_finalise(this) class(cable_netcdf_stub_io_t), intent(inout) :: this end subroutine - function cable_netcdf_stub_io_create_file(this, path) result(file) + function cable_netcdf_stub_io_create_file(this, path, iotype, mode) result(file) class(cable_netcdf_stub_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file file = cable_netcdf_stub_file_t() end function - function cable_netcdf_stub_io_open_file(this, path) result(file) + function cable_netcdf_stub_io_open_file(this, path, iotype, mode) result(file) class(cable_netcdf_stub_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file file = cable_netcdf_stub_file_t() end function @@ -130,6 +136,10 @@ subroutine cable_netcdf_stub_file_end_def(this) class(cable_netcdf_stub_file_t), intent(inout) :: this end subroutine + subroutine cable_netcdf_stub_file_redef(this) + class(cable_netcdf_stub_file_t), intent(inout) :: this + end subroutine + subroutine cable_netcdf_stub_file_sync(this) class(cable_netcdf_stub_file_t), intent(inout) :: this end subroutine @@ -142,7 +152,8 @@ subroutine cable_netcdf_stub_file_def_dims(this, dim_names, dim_lens) subroutine cable_netcdf_stub_file_def_var(this, var_name, dim_names, type) class(cable_netcdf_stub_file_t), intent(inout) :: this - character(len=*), intent(in) :: var_name, dim_names(:) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: dim_names(:) integer, intent(in) :: type end subroutine @@ -255,6 +266,13 @@ subroutine cable_netcdf_stub_file_inq_dim_len(this, dim_name, dim_len) dim_len = 0 end subroutine + subroutine cable_netcdf_stub_file_inq_var_ndims(this, var_name, ndims) + class(cable_netcdf_stub_file_t), intent(inout) :: this + character(len=*), intent(in) :: var_name + integer, intent(out) :: ndims + ndims = 0 + end subroutine + subroutine cable_netcdf_stub_file_put_var_int32_0d(this, var_name, values, start, count) class(cable_netcdf_stub_file_t), intent(inout) :: this character(len=*), intent(in) :: var_name diff --git a/src/util/netcdf/nf90/cable_netcdf_nf90.F90 b/src/util/netcdf/nf90/cable_netcdf_nf90.F90 index dfa8a08d4..b50b30605 100644 --- a/src/util/netcdf/nf90/cable_netcdf_nf90.F90 +++ b/src/util/netcdf/nf90/cable_netcdf_nf90.F90 @@ -23,8 +23,14 @@ module cable_netcdf_nf90_mod use netcdf, only: nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_enddef + use netcdf, only: nf90_redef use netcdf, only: NF90_NOERR use netcdf, only: NF90_NETCDF4 + use netcdf, only: NF90_CLASSIC_MODEL + use netcdf, only: NF90_CLOBBER + use netcdf, only: NF90_NOCLOBBER + use netcdf, only: NF90_WRITE + use netcdf, only: NF90_NOWRITE use netcdf, only: NF90_UNLIMITED use netcdf, only: NF90_INT use netcdf, only: NF90_FLOAT @@ -55,6 +61,7 @@ module cable_netcdf_nf90_mod contains procedure :: close => cable_netcdf_nf90_file_close procedure :: end_def => cable_netcdf_nf90_file_end_def + procedure :: redef => cable_netcdf_nf90_file_redef procedure :: sync => cable_netcdf_nf90_file_sync procedure :: def_dims => cable_netcdf_nf90_file_def_dims procedure :: def_var => cable_netcdf_nf90_file_def_var @@ -75,6 +82,7 @@ module cable_netcdf_nf90_mod procedure :: get_att_var_real32 => cable_netcdf_nf90_file_get_att_var_real32 procedure :: get_att_var_real64 => cable_netcdf_nf90_file_get_att_var_real64 procedure :: inq_dim_len => cable_netcdf_nf90_file_inq_dim_len + procedure :: inq_var_ndims => cable_netcdf_nf90_file_inq_var_ndims procedure :: put_var_int32_0d => cable_netcdf_nf90_file_put_var_int32_0d procedure :: put_var_int32_1d => cable_netcdf_nf90_file_put_var_int32_1d procedure :: put_var_int32_2d => cable_netcdf_nf90_file_put_var_int32_2d @@ -136,6 +144,34 @@ function type_nf90(type) end select end function type_nf90 + function cmode_nf90(iotype, mode) + integer, intent(in) :: iotype + integer, intent(in), optional :: mode + integer :: cmode_nf90 + select case(iotype) + case (CABLE_NETCDF_IOTYPE_CLASSIC) + cmode_nf90 = NF90_CLASSIC_MODEL + case (CABLE_NETCDF_IOTYPE_NETCDF4C, CABLE_NETCDF_IOTYPE_NETCDF4P) + cmode_nf90 = NF90_NETCDF4 + case default + call cable_abort("Error: iotype not supported", __FILE__, __LINE__) + end select + if (present(mode)) then + select case(mode) + case (CABLE_NETCDF_MODE_NOCLOBBER) + cmode_nf90 = ior(cmode_nf90, NF90_NOCLOBBER) + case (CABLE_NETCDF_MODE_CLOBBER) + cmode_nf90 = ior(cmode_nf90, NF90_CLOBBER) + case (CABLE_NETCDF_MODE_WRITE) + cmode_nf90 = ior(cmode_nf90, NF90_WRITE) + case (CABLE_NETCDF_MODE_NOWRITE) + cmode_nf90 = ior(cmode_nf90, NF90_NOWRITE) + case default + call cable_abort("Error: mode not supported", __FILE__, __LINE__) + end select + end if + end function cmode_nf90 + subroutine check_nf90(status) integer, intent ( in) :: status if(status /= NF90_NOERR) then @@ -151,21 +187,25 @@ subroutine cable_netcdf_nf90_io_finalise(this) class(cable_netcdf_nf90_io_t), intent(inout) :: this end subroutine - function cable_netcdf_nf90_io_create_file(this, path) result(file) + function cable_netcdf_nf90_io_create_file(this, path, iotype, mode) result(file) class(cable_netcdf_nf90_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file integer :: ncid - call check_nf90(nf90_create(path, NF90_NETCDF4, ncid)) + call check_nf90(nf90_create(path, cmode_nf90(iotype, mode), ncid)) file = cable_netcdf_nf90_file_t(ncid=ncid) end function - function cable_netcdf_nf90_io_open_file(this, path) result(file) + function cable_netcdf_nf90_io_open_file(this, path, iotype, mode) result(file) class(cable_netcdf_nf90_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file integer :: ncid - call check_nf90(nf90_open(path, NF90_NETCDF4, ncid)) + call check_nf90(nf90_open(path, cmode_nf90(iotype, mode), ncid)) file = cable_netcdf_nf90_file_t(ncid=ncid) end function @@ -187,6 +227,11 @@ subroutine cable_netcdf_nf90_file_end_def(this) call check_nf90(nf90_enddef(this%ncid)) end subroutine + subroutine cable_netcdf_nf90_file_redef(this) + class(cable_netcdf_nf90_file_t), intent(inout) :: this + call check_nf90(nf90_redef(this%ncid)) + end subroutine + subroutine cable_netcdf_nf90_file_sync(this) class(cable_netcdf_nf90_file_t), intent(inout) :: this call check_nf90(nf90_sync(this%ncid)) @@ -208,15 +253,23 @@ subroutine cable_netcdf_nf90_file_def_dims(this, dim_names, dim_lens) subroutine cable_netcdf_nf90_file_def_var(this, var_name, dim_names, type) class(cable_netcdf_nf90_file_t), intent(inout) :: this - character(len=*), intent(in) :: var_name, dim_names(:) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: dim_names(:) integer, intent(in) :: type integer, allocatable :: dimids(:) integer :: i, tmp + + if (.not. present(dim_names)) then + call check_nf90(nf90_def_var(this%ncid, var_name, type_nf90(type), tmp)) + return + end if + allocate(dimids(size(dim_names))) do i = 1, size(dimids) call check_nf90(nf90_inq_dimid(this%ncid, dim_names(i), dimids(i))) end do call check_nf90(nf90_def_var(this%ncid, var_name, type_nf90(type), dimids, tmp)) + end subroutine subroutine cable_netcdf_nf90_file_put_att_global_string(this, att_name, att_value) @@ -354,6 +407,15 @@ subroutine cable_netcdf_nf90_file_inq_dim_len(this, dim_name, dim_len) call check_nf90(nf90_inquire_dimension(this%ncid, dimid, len=dim_len)) end subroutine + subroutine cable_netcdf_nf90_file_inq_var_ndims(this, var_name, ndims) + class(cable_netcdf_nf90_file_t), intent(inout) :: this + character(len=*), intent(in) :: var_name + integer, intent(out) :: ndims + integer :: varid + call check_nf90(nf90_inq_varid(this%ncid, var_name, varid)) + call check_nf90(nf90_inquire_variable(this%ncid, varid, ndims=ndims)) + end subroutine + subroutine cable_netcdf_nf90_file_put_var_int32_0d(this, var_name, values, start, count) class(cable_netcdf_nf90_file_t), intent(inout) :: this character(len=*), intent(in) :: var_name diff --git a/src/util/netcdf/pio/cable_netcdf_pio.F90 b/src/util/netcdf/pio/cable_netcdf_pio.F90 index 3a2cacb20..550b5b537 100644 --- a/src/util/netcdf/pio/cable_netcdf_pio.F90 +++ b/src/util/netcdf/pio/cable_netcdf_pio.F90 @@ -25,8 +25,10 @@ module cable_netcdf_pio_mod use pio, only: pio_read_darray use pio, only: pio_strerror use pio, only: pio_enddef + use pio, only: pio_redef use pio, only: pio_inq_dimid use pio, only: pio_inquire_dimension + use pio, only: pio_inquire_variable use pio, only: pio_inq_varid use pio, only: pio_finalize use pio, only: PIO_MAX_NAME @@ -35,8 +37,13 @@ module cable_netcdf_pio_mod use pio, only: PIO_REAL use pio, only: PIO_DOUBLE use pio, only: PIO_REARR_BOX + use pio, only: PIO_IOTYPE_NETCDF use pio, only: PIO_IOTYPE_NETCDF4C + use pio, only: PIO_IOTYPE_NETCDF4P use pio, only: PIO_CLOBBER + use pio, only: PIO_NOCLOBBER + use pio, only: PIO_WRITE + use pio, only: PIO_NOWRITE use pio, only: PIO_UNLIMITED use pio, only: PIO_NOERR use pio, only: PIO_GLOBAL @@ -70,6 +77,7 @@ module cable_netcdf_pio_mod contains procedure :: close => cable_netcdf_pio_file_close procedure :: end_def => cable_netcdf_pio_file_end_def + procedure :: redef => cable_netcdf_pio_file_redef procedure :: sync => cable_netcdf_pio_file_sync procedure :: def_dims => cable_netcdf_pio_file_def_dims procedure :: def_var => cable_netcdf_pio_file_def_var @@ -90,6 +98,7 @@ module cable_netcdf_pio_mod procedure :: get_att_var_real32 => cable_netcdf_pio_file_get_att_var_real32 procedure :: get_att_var_real64 => cable_netcdf_pio_file_get_att_var_real64 procedure :: inq_dim_len => cable_netcdf_pio_file_inq_dim_len + procedure :: inq_var_ndims => cable_netcdf_pio_file_inq_var_ndims procedure :: put_var_int32_0d => cable_netcdf_pio_file_put_var_int32_0d procedure :: put_var_int32_1d => cable_netcdf_pio_file_put_var_int32_1d procedure :: put_var_int32_2d => cable_netcdf_pio_file_put_var_int32_2d @@ -151,6 +160,45 @@ function type_pio(basetype) end select end function type_pio + function iotype_pio(iotype) + integer, intent(in) :: iotype + integer :: iotype_pio + select case(iotype) + case(CABLE_NETCDF_IOTYPE_CLASSIC) + iotype_pio = PIO_IOTYPE_NETCDF + case(CABLE_NETCDF_IOTYPE_NETCDF4C) + iotype_pio = PIO_IOTYPE_NETCDF4C + case(CABLE_NETCDF_IOTYPE_NETCDF4P) + iotype_pio = PIO_IOTYPE_NETCDF4P + case default + call cable_abort("cable_netcdf_pio_mod: Error: iotype not supported") + end select + end function iotype_pio + + function mode_pio(mode) + integer, intent(in), optional :: mode + integer :: mode_pio + + if (.not. present(mode)) then + mode_pio = PIO_WRITE + return + end if + + select case(mode) + case(CABLE_NETCDF_MODE_CLOBBER) + mode_pio = PIO_CLOBBER + case(CABLE_NETCDF_MODE_NOCLOBBER) + mode_pio = PIO_NOCLOBBER + case(CABLE_NETCDF_MODE_WRITE) + mode_pio = PIO_WRITE + case(CABLE_NETCDF_MODE_NOWRITE) + mode_pio = PIO_NOWRITE + case default + call cable_abort("Error: mode not supported", __FILE__, __LINE__) + end select + + end function mode_pio + subroutine check_pio(status) integer, intent(in) :: status integer :: strerror_status @@ -207,21 +255,25 @@ subroutine cable_netcdf_pio_io_finalise(this) end subroutine - function cable_netcdf_pio_io_create_file(this, path) result(file) + function cable_netcdf_pio_io_create_file(this, path, iotype, mode) result(file) class(cable_netcdf_pio_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file type(pio_file_desc_t) :: pio_file_desc - call check_pio(pio_createfile(this%pio_iosystem_desc, pio_file_desc, PIO_IOTYPE_NETCDF4C, path, PIO_CLOBBER)) + call check_pio(pio_createfile(this%pio_iosystem_desc, pio_file_desc, iotype_pio(iotype), path, mode_pio(mode))) file = cable_netcdf_pio_file_t(pio_file_desc) end function - function cable_netcdf_pio_io_open_file(this, path) result(file) + function cable_netcdf_pio_io_open_file(this, path, iotype, mode) result(file) class(cable_netcdf_pio_io_t), intent(inout) :: this character(len=*), intent(in) :: path + integer, intent(in) :: iotype + integer, intent(in), optional :: mode class(cable_netcdf_file_t), allocatable :: file type(pio_file_desc_t) :: pio_file_desc - call check_pio(pio_openfile(this%pio_iosystem_desc, pio_file_desc, PIO_IOTYPE_NETCDF4C, path)) + call check_pio(pio_openfile(this%pio_iosystem_desc, pio_file_desc, iotype_pio(iotype), path, mode_pio(mode))) file = cable_netcdf_pio_file_t(pio_file_desc) end function @@ -251,6 +303,11 @@ subroutine cable_netcdf_pio_file_end_def(this) call check_pio(pio_enddef(this%pio_file_desc)) end subroutine + subroutine cable_netcdf_pio_file_redef(this) + class(cable_netcdf_pio_file_t), intent(inout) :: this + call check_pio(pio_redef(this%pio_file_desc)) + end subroutine + subroutine cable_netcdf_pio_file_sync(this) class(cable_netcdf_pio_file_t), intent(inout) :: this call pio_syncfile(this%pio_file_desc) @@ -272,16 +329,24 @@ subroutine cable_netcdf_pio_file_def_dims(this, dim_names, dim_lens) subroutine cable_netcdf_pio_file_def_var(this, var_name, dim_names, type) class(cable_netcdf_pio_file_t), intent(inout) :: this - character(len=*), intent(in) :: var_name, dim_names(:) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: dim_names(:) integer, intent(in) :: type integer, allocatable :: dimids(:) integer :: i type(pio_var_desc_t) :: tmp + + if (.not. present(dim_names)) then + call check_pio(pio_def_var(this%pio_file_desc, var_name, type_pio(type), tmp)) + return + end if + allocate(dimids(size(dim_names))) do i = 1, size(dimids) call check_pio(pio_inq_dimid(this%pio_file_desc, dim_names(i), dimids(i))) end do call check_pio(pio_def_var(this%pio_file_desc, var_name, type_pio(type), dimids, tmp)) + end subroutine subroutine cable_netcdf_pio_file_put_att_global_string(this, att_name, att_value) @@ -419,6 +484,15 @@ subroutine cable_netcdf_pio_file_inq_dim_len(this, dim_name, dim_len) call check_pio(pio_inquire_dimension(this%pio_file_desc, dimid, len=dim_len)) end subroutine + subroutine cable_netcdf_pio_file_inq_var_ndims(this, var_name, ndims) + class(cable_netcdf_pio_file_t), intent(inout) :: this + character(len=*), intent(in) :: var_name + integer, intent(out) :: ndims + integer :: varid + call check_pio(pio_inq_varid(this%pio_file_desc, var_name, varid)) + call check_pio(pio_inquire_variable(this%pio_file_desc, varid, ndims=ndims)) + end subroutine + subroutine cable_netcdf_pio_file_put_var_int32_0d(this, var_name, values, start, count) class(cable_netcdf_pio_file_t), intent(inout) :: this character(len=*), intent(in) :: var_name