diff --git a/CMakeLists.txt b/CMakeLists.txt index 9618bd2ea..428995298 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -286,10 +286,14 @@ else() src/offline/cable_mpi.F90 src/offline/cable_namelist_input.F90 src/offline/cable_output.F90 + src/offline/cable_output_prototype_v2.F90 + src/offline/cable_output_definitions.F90 src/offline/cable_parameters.F90 src/offline/cable_pft_params.F90 src/offline/cable_plume_mip.F90 src/offline/cable_read.F90 + src/offline/cable_restart.F90 + src/offline/cable_restart_write.F90 src/offline/cable_site.F90 src/offline/cable_serial.F90 src/offline/cable_soil_params.F90 @@ -299,9 +303,13 @@ else() src/offline/cbl_model_driver_offline.F90 src/offline/landuse_inout.F90 src/offline/spincasacnp.F90 + src/util/aggregator.F90 + src/util/aggregator_types.F90 src/util/cable_climate_type_mod.F90 src/util/masks_cbl.F90 src/util/cable_array_utils.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 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..2f81fb7a5 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, method="max"); CALL var%tscrn_max_daily%init() + var%tscrn_min_daily = new_aggregator(source_data=var%tscrn, method="min"); CALL var%tscrn_min_daily%init() + 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..c6e18cc89 100644 --- a/src/offline/cable_io_decomp.F90 +++ b/src/offline/cable_io_decomp.F90 @@ -12,7 +12,6 @@ module cable_io_decomp_mod use cable_io_vars_module, only: landpt use cable_io_vars_module, only: max_vegpatches use cable_io_vars_module, only: land_decomp_start - use cable_io_vars_module, only: patch_decomp_start use cable_io_vars_module, only: output use cable_io_vars_module, only: metGrid @@ -74,25 +73,6 @@ module cable_io_decomp_mod class(cable_netcdf_decomp_t), allocatable :: patch_soilcarbon_to_land_patch_soilcarbon_real32 class(cable_netcdf_decomp_t), allocatable :: patch_soilcarbon_to_land_patch_soilcarbon_real64 - class(cable_netcdf_decomp_t), allocatable :: patch_to_patch_int32 - class(cable_netcdf_decomp_t), allocatable :: patch_to_patch_real32 - class(cable_netcdf_decomp_t), allocatable :: patch_to_patch_real64 - class(cable_netcdf_decomp_t), allocatable :: patch_soil_to_patch_soil_int32 - class(cable_netcdf_decomp_t), allocatable :: patch_soil_to_patch_soil_real32 - class(cable_netcdf_decomp_t), allocatable :: patch_soil_to_patch_soil_real64 - class(cable_netcdf_decomp_t), allocatable :: patch_snow_to_patch_snow_int32 - class(cable_netcdf_decomp_t), allocatable :: patch_snow_to_patch_snow_real32 - class(cable_netcdf_decomp_t), allocatable :: patch_snow_to_patch_snow_real64 - class(cable_netcdf_decomp_t), allocatable :: patch_rad_to_patch_rad_int32 - class(cable_netcdf_decomp_t), allocatable :: patch_rad_to_patch_rad_real32 - class(cable_netcdf_decomp_t), allocatable :: patch_rad_to_patch_rad_real64 - class(cable_netcdf_decomp_t), allocatable :: patch_plantcarbon_to_patch_plantcarbon_int32 - class(cable_netcdf_decomp_t), allocatable :: patch_plantcarbon_to_patch_plantcarbon_real32 - class(cable_netcdf_decomp_t), allocatable :: patch_plantcarbon_to_patch_plantcarbon_real64 - class(cable_netcdf_decomp_t), allocatable :: patch_soilcarbon_to_patch_soilcarbon_int32 - class(cable_netcdf_decomp_t), allocatable :: patch_soilcarbon_to_patch_soilcarbon_real32 - class(cable_netcdf_decomp_t), allocatable :: patch_soilcarbon_to_patch_soilcarbon_real64 - class(cable_netcdf_decomp_t), allocatable :: land_to_x_y_int32 class(cable_netcdf_decomp_t), allocatable :: land_to_x_y_real32 class(cable_netcdf_decomp_t), allocatable :: land_to_x_y_real64 @@ -175,12 +155,6 @@ subroutine cable_io_decomp_init(io_decomp) type(dim_spec_t), allocatable :: var_shape_land_patch_rad(:) type(dim_spec_t), allocatable :: var_shape_land_patch_plantcarbon(:) type(dim_spec_t), allocatable :: var_shape_land_patch_soilcarbon(:) - type(dim_spec_t), allocatable :: var_shape_patch(:) - type(dim_spec_t), allocatable :: var_shape_patch_soil(:) - type(dim_spec_t), allocatable :: var_shape_patch_snow(:) - type(dim_spec_t), allocatable :: var_shape_patch_rad(:) - type(dim_spec_t), allocatable :: var_shape_patch_plantcarbon(:) - type(dim_spec_t), allocatable :: var_shape_patch_soilcarbon(:) logical :: requires_land_output_grid, requires_x_y_output_grid @@ -221,12 +195,6 @@ subroutine cable_io_decomp_init(io_decomp) var_shape_land_patch_rad = [dim_spec_t('land', mland_global), dim_spec_t('patch', max_vegpatches), dim_spec_t('rad', nrb)] var_shape_land_patch_plantcarbon = [dim_spec_t('land', mland_global), dim_spec_t('patch', max_vegpatches), dim_spec_t('plantcarbon', ncp)] var_shape_land_patch_soilcarbon = [dim_spec_t('land', mland_global), dim_spec_t('patch', max_vegpatches), dim_spec_t('soilcarbon', ncs)] - var_shape_patch = [dim_spec_t('patch', mp_global)] - var_shape_patch_soil = [dim_spec_t('patch', mp_global), dim_spec_t('soil', ms)] - var_shape_patch_snow = [dim_spec_t('patch', mp_global), dim_spec_t('snow', msn)] - var_shape_patch_rad = [dim_spec_t('patch', mp_global), dim_spec_t('rad', nrb)] - var_shape_patch_plantcarbon = [dim_spec_t('patch', mp_global), dim_spec_t('plantcarbon', ncp)] - var_shape_patch_soilcarbon = [dim_spec_t('patch', mp_global), dim_spec_t('soilcarbon', ncs)] io_decomp%land_to_x_y_int32 = io_decomp_land_to_x_y(land_x, land_y, mem_shape_land, var_shape_x_y, CABLE_NETCDF_INT) io_decomp%land_to_x_y_real32 = io_decomp_land_to_x_y(land_x, land_y, mem_shape_land, var_shape_x_y, CABLE_NETCDF_FLOAT) @@ -269,59 +237,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_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) - io_decomp%patch_to_patch_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch, var_shape_patch, CABLE_NETCDF_DOUBLE) - io_decomp%patch_soil_to_patch_soil_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soil, var_shape_patch_soil, CABLE_NETCDF_INT) - io_decomp%patch_soil_to_patch_soil_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soil, var_shape_patch_soil, CABLE_NETCDF_FLOAT) - io_decomp%patch_soil_to_patch_soil_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soil, var_shape_patch_soil, CABLE_NETCDF_DOUBLE) - io_decomp%patch_snow_to_patch_snow_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_snow, var_shape_patch_snow, CABLE_NETCDF_INT) - io_decomp%patch_snow_to_patch_snow_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_snow, var_shape_patch_snow, CABLE_NETCDF_FLOAT) - io_decomp%patch_snow_to_patch_snow_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_snow, var_shape_patch_snow, CABLE_NETCDF_DOUBLE) - io_decomp%patch_rad_to_patch_rad_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_rad, var_shape_patch_rad, CABLE_NETCDF_INT) - io_decomp%patch_rad_to_patch_rad_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_rad, var_shape_patch_rad, CABLE_NETCDF_FLOAT) - io_decomp%patch_rad_to_patch_rad_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_rad, var_shape_patch_rad, CABLE_NETCDF_DOUBLE) - io_decomp%patch_plantcarbon_to_patch_plantcarbon_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_plantcarbon, var_shape_patch_plantcarbon, CABLE_NETCDF_INT) - io_decomp%patch_plantcarbon_to_patch_plantcarbon_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_plantcarbon, var_shape_patch_plantcarbon, CABLE_NETCDF_FLOAT) - io_decomp%patch_plantcarbon_to_patch_plantcarbon_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_plantcarbon, var_shape_patch_plantcarbon, CABLE_NETCDF_DOUBLE) - io_decomp%patch_soilcarbon_to_patch_soilcarbon_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soilcarbon, var_shape_patch_soilcarbon, CABLE_NETCDF_INT) - io_decomp%patch_soilcarbon_to_patch_soilcarbon_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soilcarbon, var_shape_patch_soilcarbon, CABLE_NETCDF_FLOAT) - io_decomp%patch_soilcarbon_to_patch_soilcarbon_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soilcarbon, var_shape_patch_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) end subroutine diff --git a/src/offline/cable_output.F90 b/src/offline/cable_output.F90 index ed402645a..fab45e0b0 100644 --- a/src/offline/cable_output.F90 +++ b/src/offline/cable_output.F90 @@ -286,7 +286,7 @@ SUBROUTINE open_output_file(dels, soil, veg, bgc, rough, met) out_settings = output_par_settings_type(met=met, restart=.FALSE.) ! Create output file: - ok = NF90_CREATE(filename%out, NF90_CLOBBER, ncid_out) + ok = NF90_CREATE(filename%out, IOR(NF90_CLOBBER, NF90_NETCDF4), ncid_out) IF(ok /= NF90_NOERR) CALL nc_abort(ok, 'Error creating output file ' & //TRIM(filename%out)// ' (SUBROUTINE open_output_file)') ! Put the file in define mode: diff --git a/src/offline/cable_output_definitions.F90 b/src/offline/cable_output_definitions.F90 new file mode 100644 index 000000000..44cacf5ab --- /dev/null +++ b/src/offline/cable_output_definitions.F90 @@ -0,0 +1,297 @@ +module cable_output_definitions_mod + use cable_abort_module, only: cable_abort + + use cable_def_types_mod, only: canopy_type + use cable_def_types_mod, only: soil_parameter_type + + use cable_io_vars_module, only: metGrid + + use cable_netcdf_mod, only: CABLE_NETCDF_FLOAT + + use aggregator_mod, only: new_aggregator + + use cable_netcdf_mod, only: cable_netcdf_decomp_t + use cable_netcdf_mod, only: MAX_LEN_DIM => CABLE_NETCDF_MAX_STR_LEN_DIM + + use cable_io_decomp_mod, only: io_decomp_t + + use cable_output_prototype_v2_mod, only: requires_x_y_output_grid + use cable_output_prototype_v2_mod, only: requires_land_output_grid + use cable_output_prototype_v2_mod, only: output, patchout + use cable_output_prototype_v2_mod, only: cable_output_add_variable + + use cable_checks_module, only: ranges + + implicit none + private + + public :: cable_output_definitions_set + +contains + + subroutine cable_output_definitions_set(io_decomp, canopy, soil) + class(io_decomp_t), intent(in), target :: io_decomp + type(canopy_type), intent(inout) :: canopy + type(soil_parameter_type), intent(in) :: soil + + character(len=MAX_LEN_DIM), allocatable :: base_dims(:) + + 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 + + if (requires_x_y_output_grid(output%grid, metGrid)) then + base_dims = ["x", "y"] + 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 + base_dims = ["land"] + 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 + + call cable_output_add_variable( & + name="swilt", & + dims=[base_dims, "patch"], & + var_type=CABLE_NETCDF_FLOAT, & + units="1", & + long_name="", & + range=ranges%swilt, & + active=output%swilt .and. (output%patch .OR. patchout%swilt), & + parameter=.true., & + decomp=output_decomp_base_patch_real32, & + aggregator=new_aggregator( & + source_data=soil%swilt, & + method="point" & + ) & + ) + + call cable_output_add_variable( & + name="swilt", & + dims=[base_dims], & + var_type=CABLE_NETCDF_FLOAT, & + units="1", & + long_name="", & + range=ranges%swilt, & + active=output%swilt .and. .not. (output%patch .OR. patchout%swilt), & + parameter=.true., & + reduction_method="grid_cell_average", & + decomp=output_decomp_base_real32, & + aggregator=new_aggregator( & + source_data=soil%swilt, & + method="point" & + ) & + ) + + + call cable_output_add_variable( & + name="albsoil", & + dims=[base_dims, "patch", "rad"], & + var_type=CABLE_NETCDF_FLOAT, & + units="1", & + long_name="", & + range=ranges%albsoil, & + active=output%albsoil .and. (output%patch .OR. patchout%albsoil), & + parameter=.true., & + decomp=output_decomp_base_patch_rad_real32, & + aggregator=new_aggregator( & + source_data=soil%albsoil, & + method="point" & + ) & + ) + + call cable_output_add_variable( & + name="albsoil", & + dims=[base_dims, "rad"], & + var_type=CABLE_NETCDF_FLOAT, & + units="1", & + long_name="", & + range=ranges%albsoil, & + active=output%albsoil .and. .not. (output%patch .OR. patchout%albsoil), & + parameter=.true., & + reduction_method="grid_cell_average", & + decomp=output_decomp_base_rad_real32, & + aggregator=new_aggregator( & + source_data=soil%albsoil, & + method="point" & + ) & + ) + + call cable_output_add_variable( & + name="Qh", & + dims=[base_dims, "patch", "time"], & + var_type=CABLE_NETCDF_FLOAT, & + units="W/m^2", & + long_name="Surface sensible heat flux", & + range=ranges%Qh, & + active=output%Qh .and. (output%patch .OR. patchout%Qh), & + decomp=output_decomp_base_patch_real32, & + aggregator=new_aggregator( & + source_data=canopy%fh, & + method="mean" & + ) & + ) + + call cable_output_add_variable( & + name="Qh", & + dims=[base_dims, "time"], & + var_type=CABLE_NETCDF_FLOAT, & + units="W/m^2", & + long_name="Surface sensible heat flux", & + range=ranges%Qh, & + active=output%Qh .and. .not. (output%patch .OR. patchout%Qh), & + reduction_method="grid_cell_average", & + decomp=output_decomp_base_real32, & + aggregator=new_aggregator( & + source_data=canopy%fh, & + method="mean" & + ) & + ) + + call cable_output_add_variable( & + name="Tmx", & + dims=[base_dims, "patch", "time"], & + var_type=CABLE_NETCDF_FLOAT, & + units="oC", & + long_name="averaged daily maximum screen-level T", & + active=( & + output%Tex .and. & + output%averaging == "monthly" .and. & + (output%patch .OR. patchout%Tex) & + ), & + decomp=output_decomp_base_patch_real32, & + range=ranges%Tscrn, & + accumulation_frequency="daily", & + aggregator=new_aggregator( & + source_data=canopy%tscrn_max_daily%aggregated_data, & + method="mean" & + ) & + ) + + call cable_output_add_variable( & + name="Tmx", & + dims=[base_dims, "time"], & + var_type=CABLE_NETCDF_FLOAT, & + units="oC", & + long_name="averaged daily maximum screen-level T", & + active=( & + output%Tex .and. & + output%averaging == "monthly" .and. & + .not. (output%patch .OR. patchout%Tex) & + ), & + reduction_method="grid_cell_average", & + decomp=output_decomp_base_real32, & + range=ranges%Tscrn, & + accumulation_frequency="daily", & + aggregator=new_aggregator( & + source_data=canopy%tscrn_max_daily%aggregated_data, & + method="mean" & + ) & + ) + + end subroutine cable_output_definitions_set + +end module diff --git a/src/offline/cable_output_prototype_v2.F90 b/src/offline/cable_output_prototype_v2.F90 new file mode 100644 index 000000000..03a9c4215 --- /dev/null +++ b/src/offline/cable_output_prototype_v2.F90 @@ -0,0 +1,732 @@ +module cable_output_prototype_v2_mod + use iso_fortran_env, only: int32, real32, real64 + + use cable_def_types_mod, only: mp, mp_global + use cable_def_types_mod, only: mland, 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_def_types_mod, only: met_type + + use cable_abort_module, only: cable_abort + + use cable_io_vars_module, only: metGrid, patch_type, land_type, xdimsize, ydimsize, max_vegpatches + use cable_io_vars_module, only: timeunits, calendar, time_coord + use cable_io_vars_module, only: check, ON_TIMESTEP, ON_WRITE + + use cable_checks_module, only: check_range + + use cable_io_vars_module, only: output, patchout + + use cable_timing_utils_mod, only: time_step_matches + + use aggregator_mod, only: aggregator_mod_init + use aggregator_mod, only: aggregator_mod_end + use aggregator_mod, only: aggregator_t + use aggregator_mod, only: aggregator_handle_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_1d_t + use aggregator_mod, only: aggregator_real32_2d_t + use aggregator_mod, only: aggregator_real32_3d_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 aggregator_mod, only: store_aggregator + + use cable_netcdf_mod, only: cable_netcdf_file_t + use cable_netcdf_mod, only: cable_netcdf_decomp_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_UNLIMITED + use cable_netcdf_mod, only: CABLE_NETCDF_IOTYPE_NETCDF4P + use cable_netcdf_mod, only: MAX_LEN_VAR => CABLE_NETCDF_MAX_STR_LEN_VAR + use cable_netcdf_mod, only: MAX_LEN_DIM => CABLE_NETCDF_MAX_STR_LEN_DIM + + use cable_grid_reductions_mod, only: grid_cell_average + + implicit none + private + + public :: cable_output_mod_init + public :: cable_output_mod_end + public :: cable_output_add_variable + public :: cable_output_commit + public :: cable_output_update + public :: cable_output_write_parameters + public :: output + public :: patchout + public :: requires_x_y_output_grid + public :: requires_land_output_grid + + integer(kind=int32), parameter :: FILL_VALUE_INT32 = -9999_int32 + real(kind=real32), parameter :: FILL_VALUE_REAL32 = -1.0e+33_real32 + real(kind=real64), parameter :: FILL_VALUE_REAL64 = -1.0e+33_real64 + + character(len=*), parameter :: DEFAULT_ACCUMULATION_FREQUENCY = "all" + + type cable_output_variable_t + character(len=MAX_LEN_VAR) :: name + character(len=MAX_LEN_DIM), allocatable :: dims(:) + integer :: var_type + character(len=50) :: units + character(len=100) :: long_name + character(len=100) :: cell_methods + character(len=10) :: accumulation_frequency + logical :: active + logical :: parameter + character(len=50) :: reduction_method + real, dimension(2) :: range + type(aggregator_handle_t) :: aggregator_handle + class(cable_netcdf_decomp_t), pointer :: decomp => null() + real(kind=real32), pointer :: temp_buffer_real32_1d(:) => null() + real(kind=real32), pointer :: temp_buffer_real32_2d(:, :) => null() + real(kind=real32), pointer :: temp_buffer_real32_3d(:, :, :) => null() + real(kind=real64), pointer :: temp_buffer_real64_1d(:) => null() + real(kind=real64), pointer :: temp_buffer_real64_2d(:, :) => null() + real(kind=real64), pointer :: temp_buffer_real64_3d(:, :, :) => null() + end type + + type cable_output_profile_t + real :: previous_write_time = 0.0 + integer :: frame = 0 + class(cable_netcdf_file_t), allocatable :: output_file + type(cable_output_variable_t), allocatable :: output_variables(:) + + type(aggregator_handle_t), allocatable :: aggregators_accumulate_time_step(:) + type(aggregator_handle_t), allocatable :: aggregators_accumulate_daily(:) + end type + + ! Temporary buffers for computing grid-cell averages for each variable class + real(kind=real32), allocatable, target :: temp_buffer_land_real32(:) + real(kind=real64), allocatable, target :: temp_buffer_land_real64(:) + real(kind=real32), allocatable, target :: temp_buffer_land_soil_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_soil_real64(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_snow_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_snow_real64(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_rad_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_rad_real64(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_plantcarbon_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_plantcarbon_real64(:, :) + real(kind=real32), allocatable, target :: temp_buffer_land_soilcarbon_real32(:, :) + real(kind=real64), allocatable, target :: temp_buffer_land_soilcarbon_real64(:, :) + + ! TODO(Sean): once cable_write.F90 is removed, move the output_inclusion_type + ! from cable_io_vars_module to here (as this would no longer introduce a cyclic + ! module dependency). Then uncomment declarations below: + ! type(output_inclusion_t) :: output + ! type(output_inclusion_t) :: patchout ! do we want patch-specific info + + type(cable_output_profile_t), allocatable :: global_profile + +contains + + 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 check_invalid_frequency(sampling_frequency, accumulation_frequency) + character(len=*), intent(in) :: sampling_frequency + character(len=*), intent(in) :: accumulation_frequency + + check_invalid_frequency = .false. + + ! TODO(Sean): return true if sampling frequency is greater than accumulation frequency + + end function + + subroutine cable_output_mod_init() + class(cable_netcdf_file_t), allocatable :: output_file + + ! Initialize temporary buffers for grid-cell averaging + allocate(temp_buffer_land_real32(mland)) + allocate(temp_buffer_land_real64(mland)) + allocate(temp_buffer_land_soil_real32(mland, ms)) + allocate(temp_buffer_land_soil_real64(mland, ms)) + allocate(temp_buffer_land_snow_real32(mland, msn)) + allocate(temp_buffer_land_snow_real64(mland, msn)) + allocate(temp_buffer_land_rad_real32(mland, nrb)) + allocate(temp_buffer_land_rad_real64(mland, nrb)) + allocate(temp_buffer_land_plantcarbon_real32(mland, ncp)) + allocate(temp_buffer_land_plantcarbon_real64(mland, ncp)) + allocate(temp_buffer_land_soilcarbon_real32(mland, ncs)) + allocate(temp_buffer_land_soilcarbon_real64(mland, ncs)) + + call aggregator_mod_init() + + allocate(cable_output_profile_t::global_profile) + + end subroutine + + subroutine cable_output_mod_end() + + if (allocated(global_profile%output_file)) call global_profile%output_file%close() + + deallocate(global_profile) + + call aggregator_mod_end() + + deallocate(temp_buffer_land_real32) + deallocate(temp_buffer_land_real64) + deallocate(temp_buffer_land_soil_real32) + deallocate(temp_buffer_land_soil_real64) + deallocate(temp_buffer_land_snow_real32) + deallocate(temp_buffer_land_snow_real64) + deallocate(temp_buffer_land_rad_real32) + deallocate(temp_buffer_land_rad_real64) + deallocate(temp_buffer_land_plantcarbon_real32) + deallocate(temp_buffer_land_plantcarbon_real64) + deallocate(temp_buffer_land_soilcarbon_real32) + deallocate(temp_buffer_land_soilcarbon_real64) + + end subroutine + + subroutine cable_output_add_variable( & + name, dims, var_type, units, long_name, active, reduction_method, & + decomp, range, accumulation_frequency, aggregator, parameter & + ) + character(len=*), intent(in) :: name + character(len=*), dimension(:), intent(in) :: dims + integer, intent(in) :: var_type + character(len=*), intent(in) :: units + character(len=*), intent(in) :: long_name + logical, intent(in) :: active + character(len=*), intent(in), optional :: reduction_method + class(cable_netcdf_decomp_t), intent(in), target :: decomp + real, dimension(2), intent(in) :: range + character(len=*), intent(in), optional :: accumulation_frequency + class(aggregator_t), intent(in) :: aggregator + logical, intent(in), optional :: parameter + + type(cable_output_variable_t) :: output_var + + if (present(reduction_method)) then + select case (reduction_method) + case ("none") + ! No additional checks needed for 'none' reduction + case ("grid_cell_average") + select type (aggregator) + type is (aggregator_real32_1d_t) + if (size(aggregator%source_data, 1) /= mp) call cable_abort("Incompatible source data size for grid reduction", __FILE__, __LINE__) + type is (aggregator_real32_2d_t) + if (size(aggregator%source_data, 1) /= mp) call cable_abort("Incompatible source data size for grid reduction", __FILE__, __LINE__) + type is (aggregator_real32_3d_t) + if (size(aggregator%source_data, 1) /= mp) call cable_abort("Incompatible source data size for grid reduction", __FILE__, __LINE__) + type is (aggregator_real64_1d_t) + if (size(aggregator%source_data, 1) /= mp) call cable_abort("Incompatible source data size for grid reduction", __FILE__, __LINE__) + type is (aggregator_real64_2d_t) + if (size(aggregator%source_data, 1) /= mp) call cable_abort("Incompatible source data size for grid reduction", __FILE__, __LINE__) + type is (aggregator_real64_3d_t) + if (size(aggregator%source_data, 1) /= mp) call cable_abort("Incompatible source data size for grid reduction", __FILE__, __LINE__) + class default + call cable_abort("Unexpected aggregator type", __FILE__, __LINE__) + end select + case default + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end select + end if + + ! TODO(Sean): determine cell_methods based on reduction and aggregation method + + output_var%name = name + output_var%dims = dims + output_var%units = units + output_var%long_name = long_name + output_var%active = active + output_var%range = range + output_var%decomp => decomp + output_var%var_type = var_type + output_var%aggregator_handle = store_aggregator(aggregator) + + if (present(parameter)) then + output_var%parameter = parameter + else + output_var%parameter = .false. + end if + + if (present(reduction_method)) then + output_var%reduction_method = reduction_method + else + output_var%reduction_method = "none" + end if + + if (present(accumulation_frequency)) then + output_var%accumulation_frequency = accumulation_frequency + else + output_var%accumulation_frequency = DEFAULT_ACCUMULATION_FREQUENCY + end if + + if (reduction_method /= "none") then + select type(aggregator) + type is (aggregator_real32_1d_t) + if (all(shape(aggregator%source_data) == [mp])) then + output_var%temp_buffer_real32_1d => temp_buffer_land_real32 + else + call cable_abort("Unexpected source data shape for grid reduction", __FILE__, __LINE__) + end if + type is (aggregator_real64_1d_t) + if (all(shape(aggregator%source_data) == [mp])) then + output_var%temp_buffer_real64_1d => temp_buffer_land_real64 + else + call cable_abort("Unexpected source data shape for grid reduction", __FILE__, __LINE__) + end if + type is (aggregator_real32_2d_t) + if (all(shape(aggregator%source_data) == [mp, ms])) then + output_var%temp_buffer_real32_2d => temp_buffer_land_soil_real32 + else if (all(shape(aggregator%source_data) == [mp, nrb])) then + output_var%temp_buffer_real32_2d => temp_buffer_land_rad_real32 + else if (all(shape(aggregator%source_data) == [mp, msn])) then + output_var%temp_buffer_real32_2d => temp_buffer_land_snow_real32 + else if (all(shape(aggregator%source_data) == [mp, nrb])) then + output_var%temp_buffer_real32_2d => temp_buffer_land_rad_real32 + else if (all(shape(aggregator%source_data) == [mp, ncp])) then + output_var%temp_buffer_real32_2d => temp_buffer_land_plantcarbon_real32 + else if (all(shape(aggregator%source_data) == [mp, ncs])) then + output_var%temp_buffer_real32_2d => temp_buffer_land_soilcarbon_real32 + else + call cable_abort("Unexpected source data shape for grid reduction", __FILE__, __LINE__) + end if + type is (aggregator_real64_2d_t) + if (all(shape(aggregator%source_data) == [mp, ms])) then + output_var%temp_buffer_real64_2d => temp_buffer_land_soil_real64 + else if (all(shape(aggregator%source_data) == [mp, nrb])) then + output_var%temp_buffer_real64_2d => temp_buffer_land_rad_real64 + else if (all(shape(aggregator%source_data) == [mp, msn])) then + output_var%temp_buffer_real64_2d => temp_buffer_land_snow_real64 + else if (all(shape(aggregator%source_data) == [mp, nrb])) then + output_var%temp_buffer_real64_2d => temp_buffer_land_rad_real64 + else if (all(shape(aggregator%source_data) == [mp, ncp])) then + output_var%temp_buffer_real64_2d => temp_buffer_land_plantcarbon_real64 + else if (all(shape(aggregator%source_data) == [mp, ncs])) then + output_var%temp_buffer_real64_2d => temp_buffer_land_soilcarbon_real64 + else + call cable_abort("Unexpected source data shape for grid reduction", __FILE__, __LINE__) + end if + class default + call cable_abort("Unexpected aggregator type", __FILE__, __LINE__) + end select + end if + + if (.not. allocated(global_profile%output_variables)) then + global_profile%output_variables = [output_var] + else + global_profile%output_variables = [global_profile%output_variables, output_var] + end if + + end subroutine cable_output_add_variable + + subroutine cable_output_commit() + class(cable_netcdf_file_t), allocatable :: output_file + integer :: i + + output_file = cable_netcdf_create_file("test_output.nc", iotype=CABLE_NETCDF_IOTYPE_NETCDF4P) ! TODO(Sean): use filename from namelist + + call output_file%def_dims(["x", "y"], [xdimsize, ydimsize]) + call output_file%def_dims(["patch"], [max_vegpatches]) + call output_file%def_dims(["soil"], [ms]) + call output_file%def_dims(["rad"], [nrb]) + call output_file%def_dims(["soil_carbon_pools"], [ncs]) + call output_file%def_dims(["plant_carbon_pools"], [ncp]) + call output_file%def_dims(["time"], [CABLE_NETCDF_UNLIMITED]) + call output_file%def_dims(["nv"], [2]) + + if (requires_x_y_output_grid(output%grid, metgrid)) then + call output_file%def_dims(["z"], [1]) ! Atmospheric 'z' dim of size 1 to comply with ALMA grid type + else if (requires_land_output_grid(output%grid, metgrid)) then + call output_file%def_dims(["land"], [mland_global]) + call output_file%def_var("local_lat", ["land"], CABLE_NETCDF_FLOAT) + call output_file%put_att("local_lat", "units", "degrees_north") + call output_file%def_var("local_lon", ["land"], CABLE_NETCDF_FLOAT) + call output_file%put_att("local_lon", "units", "degrees_east") + else + call cable_abort("Error: Unable to determine output grid type", __FILE__, __LINE__) + 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) + call output_file%put_att("time", "bounds", "time_bnds") + call output_file%def_var("time_bnds", ["nv", "time"], CABLE_NETCDF_DOUBLE) + + ! Define latitude and longitude variable (ALMA): + call output_file%def_var("latitude", ["x", "y"], CABLE_NETCDF_FLOAT) + call output_file%put_att("latitude", "units", "degrees_north") + call output_file%def_var("longitude", ["x", "y"], CABLE_NETCDF_FLOAT) + call output_file%put_att("longitude", "units", "degrees_east") + + ! Write "cordinate variables" to enable reading by GrADS: + call output_file%def_var("x", ["x"], CABLE_NETCDF_FLOAT) + call output_file%put_att("x", "units", "degrees_east") + call output_file%put_att("x", "comment", "x coordinate variable for GrADS compatibility") + call output_file%def_var("y", ["y"], CABLE_NETCDF_FLOAT) + call output_file%put_att("y", "units", "degrees_north") + call output_file%put_att("y", "comment", "y coordinate variable for GrADS compatibility") + + ! TODO(Sean): define remaining coordinate variables + + ! TODO(Sean): add global attributes + + global_profile%output_variables = pack(global_profile%output_variables, global_profile%output_variables(:)%active) + + do i = 1, size(global_profile%output_variables) + associate(output_var => global_profile%output_variables(i)) + if (check_invalid_frequency( & + sampling_frequency=output%averaging, & + accumulation_frequency=output_var%accumulation_frequency & + )) then + call cable_abort("Sampling frequency and accumulation frequency are incompatible", __FILE__, __LINE__) + end if + end associate + end do + + do i = 1, size(global_profile%output_variables) + associate(output_var => global_profile%output_variables(i)) + call output_file%def_var( & + var_name=output_var%name, & + dim_names=output_var%dims, & + type=output_var%var_type & + ) + call output_file%put_att(output_var%name, 'units', output_var%units) + call output_file%put_att(output_var%name, 'long_name', output_var%long_name) + 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 + ! TODO(Sean): set cell_methods attribute + end associate + end do + + call output_file%end_def() + + global_profile%output_file = output_file + + ! Initialise aggregators and accumulation lists + + do i = 1, size(global_profile%output_variables) + associate(output_var => global_profile%output_variables(i)) + call output_var%aggregator_handle%init() + if (output_var%parameter) cycle ! Register only time-varying variables for accumulation + select case(output_var%accumulation_frequency) + case("all") + if (.not. allocated(global_profile%aggregators_accumulate_time_step)) then + global_profile%aggregators_accumulate_time_step = [output_var%aggregator_handle] + else + global_profile%aggregators_accumulate_time_step = [global_profile%aggregators_accumulate_time_step, output_var%aggregator_handle] + end if + case("daily") + if (.not. allocated(global_profile%aggregators_accumulate_daily)) then + global_profile%aggregators_accumulate_daily = [output_var%aggregator_handle] + else + global_profile%aggregators_accumulate_daily = [global_profile%aggregators_accumulate_daily, output_var%aggregator_handle] + end if + case default + call cable_abort("Invalid accumulation frequency", __FILE__, __LINE__) + end select + 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_handle%accumulate() + select case (output_variable%reduction_method) + case ("grid_cell_average") + call write_variable_grid_cell_average(output_variable, global_profile%output_file, patch, landpt) + case ("none") + call write_variable(output_variable, global_profile%output_file) + case default + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end select + call output_variable%aggregator_handle%reset() + end associate + end do + + end subroutine cable_output_write_parameters + + subroutine cable_output_update(time_index, dels, leaps, start_year, patch, landpt, met) + integer, intent(in) :: time_index + real, intent(in) :: dels + logical, intent(in) :: leaps + integer, intent(in) :: start_year + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + 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%aggregators_accumulate_time_step) + associate(aggregator_handle => global_profile%aggregators_accumulate_time_step(i)) + call aggregator_handle%accumulate() + end associate + end do + + if (time_step_matches(dels, time_index, "daily", leaps, start_year)) then + do i = 1, size(global_profile%aggregators_accumulate_daily) + associate(aggregator_handle => global_profile%aggregators_accumulate_daily(i)) + call aggregator_handle%accumulate() + end associate + end do + end if + + if (time_step_matches(dels, time_index, output%averaging, 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) + select case (output_variable%reduction_method) + case ("grid_cell_average") + call write_variable_grid_cell_average(output_variable, global_profile%output_file, patch, landpt, global_profile%frame + 1) + case ("none") + call write_variable(output_variable, global_profile%output_file, global_profile%frame + 1) + case default + call cable_abort("Invalid reduction method", __FILE__, __LINE__) + end select + call output_variable%aggregator_handle%reset() + end associate + end do + + current_time = time_index * dels + call global_profile%output_file%put_var("time", (current_time + global_profile%previous_write_time) / 2.0, start=[global_profile%frame + 1]) + 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_update + + 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_handle%aggregator) + 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_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_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, time_index) + type(cable_output_variable_t), intent(inout) :: output_variable + class(cable_netcdf_file_t), intent(inout) :: output_file + integer, intent(in), optional :: time_index + + select type (aggregator => output_variable%aggregator_handle%aggregator) + type is (aggregator_int32_1d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + frame=time_index) + type is (aggregator_int32_2d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + frame=time_index) + type is (aggregator_int32_3d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + frame=time_index) + type is (aggregator_real32_1d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=time_index) + type is (aggregator_real32_2d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=time_index) + type is (aggregator_real32_3d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=time_index) + type is (aggregator_real64_1d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=time_index) + type is (aggregator_real64_2d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=time_index) + type is (aggregator_real64_3d_t) + call output_file%write_darray( & + var_name=output_variable%name, & + values=aggregator%aggregated_data, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=time_index) + class default + call cable_abort("Unexpected aggregator type", __FILE__, __LINE__) + end select + + end subroutine write_variable + + subroutine write_variable_grid_cell_average(output_variable, output_file, patch, landpt, time_index) + type(cable_output_variable_t), intent(inout) :: output_variable + class(cable_netcdf_file_t), intent(inout) :: output_file + type(patch_type), intent(in) :: patch(:) + type(land_type), intent(in) :: landpt(:) + integer, intent(in), optional :: time_index + + select type (aggregator => output_variable%aggregator_handle%aggregator) + type is (aggregator_real32_1d_t) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=output_variable%temp_buffer_real32_1d, & + landpt=landpt, & + patch=patch) + call output_file%write_darray( & + var_name=output_variable%name, & + values=output_variable%temp_buffer_real32_1d, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=time_index) + type is (aggregator_real32_2d_t) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=output_variable%temp_buffer_real32_2d, & + landpt=landpt, & + patch=patch) + call output_file%write_darray( & + var_name=output_variable%name, & + values=output_variable%temp_buffer_real32_2d, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=time_index) + type is (aggregator_real32_3d_t) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=output_variable%temp_buffer_real32_3d, & + landpt=landpt, & + patch=patch) + call output_file%write_darray( & + var_name=output_variable%name, & + values=output_variable%temp_buffer_real32_3d, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL32, & + frame=time_index) + type is (aggregator_real64_1d_t) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=output_variable%temp_buffer_real64_1d, & + landpt=landpt, & + patch=patch) + call output_file%write_darray( & + var_name=output_variable%name, & + values=output_variable%temp_buffer_real64_1d, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=time_index) + type is (aggregator_real64_2d_t) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=output_variable%temp_buffer_real64_2d, & + landpt=landpt, & + patch=patch) + call output_file%write_darray( & + var_name=output_variable%name, & + values=output_variable%temp_buffer_real64_2d, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=time_index) + type is (aggregator_real64_3d_t) + call grid_cell_average( & + input_array=aggregator%aggregated_data, & + output_array=output_variable%temp_buffer_real64_3d, & + landpt=landpt, & + patch=patch) + call output_file%write_darray( & + var_name=output_variable%name, & + values=output_variable%temp_buffer_real64_3d, & + decomp=output_variable%decomp, & + fill_value=FILL_VALUE_REAL64, & + frame=time_index) + class default + call cable_abort("Unexpected aggregator type", __FILE__, __LINE__) + end select + + end subroutine write_variable_grid_cell_average + +end module diff --git a/src/offline/cable_restart.F90 b/src/offline/cable_restart.F90 new file mode 100644 index 000000000..6182ebfef --- /dev/null +++ b/src/offline/cable_restart.F90 @@ -0,0 +1,501 @@ +module cable_restart_mod + use iso_fortran_env, only: int32, real32, real64 + + use cable_abort_module, only: cable_abort + + use cable_def_types_mod, only: mp, mp_global + use cable_def_types_mod, only: mland_global + use cable_def_types_mod, only: ms, msn, nrb, ncp, ncs + + use cable_io_vars_module, only: patch_decomp_start + use cable_io_vars_module, only: patch + use cable_io_vars_module, only: timeunits, calendar, time_coord + + 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_IOTYPE_CLASSIC + 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_decomp_util_mod, only: dim_spec_t + use cable_netcdf_decomp_util_mod, only: io_decomp_patch_to_patch + implicit none + private + + public :: cable_restart_mod_init + public :: cable_restart_mod_end + public :: cable_restart_variable_write + public :: cable_restart_variable_write_darray + public :: cable_restart_write_time + + ! TODO(Sean): is an interface overkill here? It does make things more intuitive for distributed I/O + + interface cable_restart_variable_write_darray + module procedure cable_restart_variable_write_darray_int32_1d + module procedure cable_restart_variable_write_darray_int32_2d + module procedure cable_restart_variable_write_darray_int32_3d + module procedure cable_restart_variable_write_darray_real32_1d + module procedure cable_restart_variable_write_darray_real32_2d + module procedure cable_restart_variable_write_darray_real32_3d + module procedure cable_restart_variable_write_darray_real64_1d + module procedure cable_restart_variable_write_darray_real64_2d + module procedure cable_restart_variable_write_darray_real64_3d + end interface cable_restart_variable_write_darray + + interface cable_restart_variable_write + module procedure cable_restart_variable_write_int32_1d + module procedure cable_restart_variable_write_int32_2d + module procedure cable_restart_variable_write_int32_3d + module procedure cable_restart_variable_write_real32_1d + module procedure cable_restart_variable_write_real32_2d + module procedure cable_restart_variable_write_real32_3d + module procedure cable_restart_variable_write_real64_1d + module procedure cable_restart_variable_write_real64_2d + module procedure cable_restart_variable_write_real64_3d + end interface cable_restart_variable_write + + class(cable_netcdf_file_t), allocatable :: restart_output_file + + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_int32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_real32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_real64 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_soil_int32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_soil_real32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_soil_real64 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_snow_int32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_snow_real32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_snow_real64 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_rad_int32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_rad_real32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_rad_real64 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_plantcarbon_int32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_plantcarbon_real32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_plantcarbon_real64 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_soilcarbon_int32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_soilcarbon_real32 + class(cable_netcdf_decomp_t), allocatable, target :: decomp_patch_soilcarbon_real64 + +contains + + subroutine cable_restart_mod_init() + + type(dim_spec_t), allocatable :: mem_shape_patch(:) + type(dim_spec_t), allocatable :: mem_shape_patch_soil(:) + type(dim_spec_t), allocatable :: mem_shape_patch_snow(:) + type(dim_spec_t), allocatable :: mem_shape_patch_rad(:) + type(dim_spec_t), allocatable :: mem_shape_patch_plantcarbon(:) + type(dim_spec_t), allocatable :: mem_shape_patch_soilcarbon(:) + type(dim_spec_t), allocatable :: var_shape_patch(:) + type(dim_spec_t), allocatable :: var_shape_patch_soil(:) + type(dim_spec_t), allocatable :: var_shape_patch_snow(:) + type(dim_spec_t), allocatable :: var_shape_patch_rad(:) + type(dim_spec_t), allocatable :: var_shape_patch_plantcarbon(:) + type(dim_spec_t), allocatable :: var_shape_patch_soilcarbon(:) + + mem_shape_patch = [dim_spec_t('patch', mp)] + mem_shape_patch_soil = [dim_spec_t('patch', mp), dim_spec_t('soil', ms)] + mem_shape_patch_snow = [dim_spec_t('patch', mp), dim_spec_t('snow', msn)] + mem_shape_patch_rad = [dim_spec_t('patch', mp), dim_spec_t('rad', nrb)] + mem_shape_patch_plantcarbon = [dim_spec_t('patch', mp), dim_spec_t('plantcarbon', ncp)] + mem_shape_patch_soilcarbon = [dim_spec_t('patch', mp), dim_spec_t('soilcarbon', ncs)] + + var_shape_patch = [dim_spec_t('patch', mp_global)] + var_shape_patch_soil = [dim_spec_t('patch', mp_global), dim_spec_t('soil', ms)] + var_shape_patch_snow = [dim_spec_t('patch', mp_global), dim_spec_t('snow', msn)] + var_shape_patch_rad = [dim_spec_t('patch', mp_global), dim_spec_t('rad', nrb)] + var_shape_patch_plantcarbon = [dim_spec_t('patch', mp_global), dim_spec_t('plantcarbon', ncp)] + var_shape_patch_soilcarbon = [dim_spec_t('patch', mp_global), dim_spec_t('soilcarbon', ncs)] + + decomp_patch_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch, var_shape_patch, CABLE_NETCDF_INT) + decomp_patch_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch, var_shape_patch, CABLE_NETCDF_FLOAT) + decomp_patch_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch, var_shape_patch, CABLE_NETCDF_DOUBLE) + decomp_patch_soil_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soil, var_shape_patch_soil, CABLE_NETCDF_INT) + decomp_patch_soil_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soil, var_shape_patch_soil, CABLE_NETCDF_FLOAT) + decomp_patch_soil_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soil, var_shape_patch_soil, CABLE_NETCDF_DOUBLE) + decomp_patch_snow_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_snow, var_shape_patch_snow, CABLE_NETCDF_INT) + decomp_patch_snow_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_snow, var_shape_patch_snow, CABLE_NETCDF_FLOAT) + decomp_patch_snow_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_snow, var_shape_patch_snow, CABLE_NETCDF_DOUBLE) + decomp_patch_rad_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_rad, var_shape_patch_rad, CABLE_NETCDF_INT) + decomp_patch_rad_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_rad, var_shape_patch_rad, CABLE_NETCDF_FLOAT) + decomp_patch_rad_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_rad, var_shape_patch_rad, CABLE_NETCDF_DOUBLE) + decomp_patch_plantcarbon_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_plantcarbon, var_shape_patch_plantcarbon, CABLE_NETCDF_INT) + decomp_patch_plantcarbon_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_plantcarbon, var_shape_patch_plantcarbon, CABLE_NETCDF_FLOAT) + decomp_patch_plantcarbon_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_plantcarbon, var_shape_patch_plantcarbon, CABLE_NETCDF_DOUBLE) + decomp_patch_soilcarbon_int32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soilcarbon, var_shape_patch_soilcarbon, CABLE_NETCDF_INT) + decomp_patch_soilcarbon_real32 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soilcarbon, var_shape_patch_soilcarbon, CABLE_NETCDF_FLOAT) + decomp_patch_soilcarbon_real64 = io_decomp_patch_to_patch(patch_decomp_start, mem_shape_patch_soilcarbon, var_shape_patch_soilcarbon, CABLE_NETCDF_DOUBLE) + + restart_output_file = cable_netcdf_create_file("test_restart.nc", iotype=CABLE_NETCDF_IOTYPE_CLASSIC) ! TODO(Sean): use filename from namelist + + call restart_output_file%def_dims(["mland"], [mland_global]) + call restart_output_file%def_dims(["mp"], [mp_global]) + call restart_output_file%def_dims(["soil"], [ms]) + call restart_output_file%def_dims(["snow"], [msn]) + call restart_output_file%def_dims(["rad"], [nrb]) + call restart_output_file%def_dims(["soil_carbon_pools"], [ncs]) + call restart_output_file%def_dims(["plant_carbon_pools"], [ncp]) + call restart_output_file%def_dims(["time"], [1]) + + call restart_output_file%end_def() + + end subroutine cable_restart_mod_init + + subroutine cable_restart_mod_end() + + if (allocated(restart_output_file)) call restart_output_file%close() + + end subroutine cable_restart_mod_end + + subroutine define_variable(output_file, var_name, var_dims, var_type, long_name, units) + class(cable_netcdf_file_t), intent(inout) :: output_file + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call output_file%redef() + call output_file%def_var(var_name, var_dims, var_type) + call output_file%put_att(var_name, "long_name", long_name) + call output_file%put_att(var_name, "units", units) + call output_file%end_def() + + end subroutine define_variable + + subroutine associate_decomp_int32(var_name, decomp, data_shape) + character(len=*), intent(in) :: var_name + class(cable_netcdf_decomp_t), pointer, intent(inout) :: decomp + integer, dimension(:), intent(in) :: data_shape + + if (all(data_shape == [mp])) then + decomp => decomp_patch_int32 + else if (all(data_shape == [mp, ms])) then + decomp => decomp_patch_soil_int32 + else if (all(data_shape == [mp, msn])) then + decomp => decomp_patch_snow_int32 + else if (all(data_shape == [mp, nrb])) then + decomp => decomp_patch_rad_int32 + else if (all(data_shape == [mp, ncp])) then + decomp => decomp_patch_plantcarbon_int32 + else if (all(data_shape == [mp, ncs])) then + decomp => decomp_patch_soilcarbon_int32 + else + call cable_abort("Unexpected data shape for variable " // var_name, __FILE__, __LINE__) + end if + + end subroutine associate_decomp_int32 + + subroutine associate_decomp_real32(var_name, decomp, data_shape) + character(len=*), intent(in) :: var_name + class(cable_netcdf_decomp_t), pointer, intent(inout) :: decomp + integer, dimension(:), intent(in) :: data_shape + + if (all(data_shape == [mp])) then + decomp => decomp_patch_real32 + else if (all(data_shape == [mp, ms])) then + decomp => decomp_patch_soil_real32 + else if (all(data_shape == [mp, msn])) then + decomp => decomp_patch_snow_real32 + else if (all(data_shape == [mp, nrb])) then + decomp => decomp_patch_rad_real32 + else if (all(data_shape == [mp, ncp])) then + decomp => decomp_patch_plantcarbon_real32 + else if (all(data_shape == [mp, ncs])) then + decomp => decomp_patch_soilcarbon_real32 + else + call cable_abort("Unexpected data shape for variable " // var_name, __FILE__, __LINE__) + end if + + end subroutine associate_decomp_real32 + + subroutine associate_decomp_real64(var_name, decomp, data_shape) + character(len=*), intent(in) :: var_name + class(cable_netcdf_decomp_t), pointer, intent(inout) :: decomp + integer, dimension(:), intent(in) :: data_shape + + if (all(data_shape == [mp])) then + decomp => decomp_patch_real64 + else if (all(data_shape == [mp, ms])) then + decomp => decomp_patch_soil_real64 + else if (all(data_shape == [mp, msn])) then + decomp => decomp_patch_snow_real64 + else if (all(data_shape == [mp, nrb])) then + decomp => decomp_patch_rad_real64 + else if (all(data_shape == [mp, ncp])) then + decomp => decomp_patch_plantcarbon_real64 + else if (all(data_shape == [mp, ncs])) then + decomp => decomp_patch_soilcarbon_real64 + else + call cable_abort("Unexpected data shape for variable " // var_name, __FILE__, __LINE__) + end if + + end subroutine associate_decomp_real64 + + subroutine cable_restart_variable_write_darray_int32_1d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + integer(kind=int32), intent(in) :: data(:) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_int32(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_int32_1d + + subroutine cable_restart_variable_write_darray_int32_2d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + integer(kind=int32), intent(in) :: data(:, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_int32(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_int32_2d + + subroutine cable_restart_variable_write_darray_int32_3d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + integer(kind=int32), intent(in) :: data(:, :, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_int32(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_int32_3d + + subroutine cable_restart_variable_write_darray_real32_1d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real32), intent(in) :: data(:) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_real32(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_real32_1d + + subroutine cable_restart_variable_write_darray_real32_2d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real32), intent(in) :: data(:, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_real32(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_real32_2d + + subroutine cable_restart_variable_write_darray_real32_3d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real32), intent(in) :: data(:, :, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_real32(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_real32_3d + + subroutine cable_restart_variable_write_darray_real64_1d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real64), intent(in) :: data(:) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_real32(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_real64_1d + + subroutine cable_restart_variable_write_darray_real64_2d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real64), intent(in) :: data(:, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_real64(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_real64_2d + + subroutine cable_restart_variable_write_darray_real64_3d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real64), intent(in) :: data(:, :, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + class(cable_netcdf_decomp_t), pointer :: decomp + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call associate_decomp_real64(var_name, decomp, shape(data)) + call restart_output_file%write_darray(var_name, data, decomp) + + end subroutine cable_restart_variable_write_darray_real64_3d + + subroutine cable_restart_variable_write_int32_1d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + integer(kind=int32), intent(in) :: data(:) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_int32_1d + + subroutine cable_restart_variable_write_int32_2d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + integer(kind=int32), intent(in) :: data(:, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_int32_2d + + subroutine cable_restart_variable_write_int32_3d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + integer(kind=int32), intent(in) :: data(:, :, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_int32_3d + + subroutine cable_restart_variable_write_real32_1d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real32), intent(in) :: data(:) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_real32_1d + + subroutine cable_restart_variable_write_real32_2d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real32), intent(in) :: data(:, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_real32_2d + + subroutine cable_restart_variable_write_real32_3d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real32), intent(in) :: data(:, :, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_real32_3d + + subroutine cable_restart_variable_write_real64_1d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real64), intent(in) :: data(:) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_real64_1d + + subroutine cable_restart_variable_write_real64_2d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real64), intent(in) :: data(:, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_real64_2d + + subroutine cable_restart_variable_write_real64_3d(var_name, var_dims, data, var_type, long_name, units) + character(len=*), intent(in) :: var_name + character(len=*), intent(in), optional :: var_dims(:) + real(kind=real64), intent(in) :: data(:, :, :) + integer, intent(in) :: var_type + character(len=*), intent(in) :: long_name + character(len=*), intent(in) :: units + + call define_variable(restart_output_file, var_name, var_dims, var_type, long_name, units) + call restart_output_file%put_var(var_name, data) + + end subroutine cable_restart_variable_write_real64_3d + + subroutine cable_restart_write_time(time_value) + real, intent(in) :: time_value + + call restart_output_file%redef() + call restart_output_file%def_var("time", ["time"], CABLE_NETCDF_DOUBLE) + call restart_output_file%put_att("time", "units", timeunits) + call restart_output_file%put_att("time", "coordinate", time_coord) + call restart_output_file%put_att("time", "calendar", calendar) + call restart_output_file%end_def() + call restart_output_file%put_var("time", [time_value]) + + end subroutine cable_restart_write_time + +end module diff --git a/src/offline/cable_restart_write.F90 b/src/offline/cable_restart_write.F90 new file mode 100644 index 000000000..0a7410771 --- /dev/null +++ b/src/offline/cable_restart_write.F90 @@ -0,0 +1,568 @@ +module cable_restart_write_mod + use cable_restart_mod, only: cable_restart_write_time + use cable_restart_mod, only: cable_restart_variable_write + use cable_restart_mod, only: cable_restart_variable_write_darray + + use cable_common_module, only: cable_user + + use cable_def_types_mod, only: met_type + use cable_def_types_mod, only: soil_parameter_type + use cable_def_types_mod, only: veg_parameter_type + use cable_def_types_mod, only: soil_snow_type + use cable_def_types_mod, only: bgc_pool_type + use cable_def_types_mod, only: canopy_type + use cable_def_types_mod, only: roughness_type + use cable_def_types_mod, only: radiation_type + use cable_def_types_mod, only: balances_type + use cable_def_types_mod, only: mvtype, mstype + + use cable_io_vars_module, only: latitude, longitude + use cable_io_vars_module, only: landpt_global + use cable_io_vars_module, only: patch + + use cable_netcdf_mod, only: CABLE_NETCDF_INT + use cable_netcdf_mod, only: CABLE_NETCDF_FLOAT + + implicit none + private + + public :: cable_restart_write + +contains + + subroutine cable_restart_write(current_time, soil, veg, ssnow, canopy, rough, rad, bgc, bal, met) + real, intent(in) :: current_time !! Current simulation time + type(met_type), intent(in) :: met !! Meteorological data + type(soil_parameter_type), intent(in) :: soil !! Soil parameters + type(veg_parameter_type), intent(in) :: veg !! Vegetation parameters + type(soil_snow_type), intent(in) :: ssnow !! Soil and snow variables + type(bgc_pool_type), intent(in) :: bgc !! Carbon pool variables + type(canopy_type), intent(in) :: canopy !! Vegetation variables + type(roughness_type), intent(in) :: rough !! Roughness variables + type(radiation_type), intent(in) :: rad !! Radiation variables + type(balances_type), intent(in) :: bal !! Energy and water balance variables + + call cable_restart_write_time(current_time) + + call cable_restart_variable_write( & + var_name="longitude", & + var_dims=["mland"], & + data=longitude, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="", & + units="degrees_east" & + ) + + call cable_restart_variable_write( & + var_name="latitude", & + var_dims=["mland"], & + data=latitude, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="", & + units="degrees_north" & + ) + + call cable_restart_variable_write( & + var_name="nap", & + var_dims=["mland"], & + data=landpt_global(:)%nap, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Number of active patches", & + units="" & + ) + + call cable_restart_variable_write_darray( & + var_name="patchfrac", & + var_dims=["mp"], & + data=patch(:)%frac, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Fraction of vegetated grid cell area occupied by a vegetation/soil patch", & + units="" & + ) + + call cable_restart_variable_write( & + var_name="mvtype", & + data=[mvtype], & + var_type=CABLE_NETCDF_INT, & + long_name="Number of vegetation types", & + units="" & + ) + + call cable_restart_variable_write( & + var_name="mstype", & + data=[mstype], & + var_type=CABLE_NETCDF_INT, & + long_name="Number of soil types", & + units="" & + ) + + call cable_restart_variable_write_darray( & + var_name="tgg", & + var_dims=["mp","soil"], & + data=ssnow%tgg, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Average layer soil temperature", & + units="K" & + ) + + call cable_restart_variable_write_darray( & + var_name="wb", & + var_dims=["mp", "soil"], & + data=ssnow%wb, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Average layer volumetric soil moisture", & + units="vol/vol" & + ) + + call cable_restart_variable_write_darray( & + var_name="wbice", & + var_dims=["mp","soil"], & + data=ssnow%wbice, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Average layer volumetric soil ice", & + units="vol/vol" & + ) + + call cable_restart_variable_write_darray( & + var_name="tss", & + var_dims=["mp"], & + data=ssnow%tss, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Combined soil/snow temperature", & + units="K" & + ) + + call cable_restart_variable_write_darray( & + var_name="albsoilsn", & + var_dims=["mp", "rad"], & + data=ssnow%albsoilsn, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Combined soil/snow albedo", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="rtsoil", & + var_dims=["mp"], & + data=ssnow%rtsoil, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Turbulent resistance for soil", & + units="??" & + ) + + call cable_restart_variable_write_darray( & + var_name="gammzz", & + var_dims=["mp","soil"], & + data=ssnow%gammzz, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Heat capacity for each soil layer", & + units="J/kg/C" & + ) + + call cable_restart_variable_write_darray( & + var_name="runoff", & + var_dims=["mp"], & + data=ssnow%runoff, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Total runoff", & + units="mm/timestep" & + ) + + call cable_restart_variable_write_darray( & + var_name="rnof1", & + var_dims=["mp"], & + data=ssnow%rnof1, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Surface runoff", & + units="mm/timestep" & + ) + + call cable_restart_variable_write_darray( & + var_name="rnof2", & + var_dims=["mp"], & + data=ssnow%rnof2, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Subsurface runoff", & + units="mm/timestep" & + ) + + call cable_restart_variable_write_darray( & + var_name="tggsn", & + var_dims=["mp","snow"], & + data=ssnow%tggsn, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Average layer snow temperature", & + units="K" & + ) + + call cable_restart_variable_write_darray( & + var_name="ssdnn", & + var_dims=["mp"], & + data=ssnow%ssdnn, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Average snow density", & + units="kg/m^3" & + ) + + call cable_restart_variable_write_darray( & + var_name="ssdn", & + var_dims=["mp","snow"], & + data=ssnow%ssdn, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Average layer snow density", & + units="kg/m^3" & + ) + + call cable_restart_variable_write_darray( & + var_name="snowd", & + var_dims=["mp"], & + data=ssnow%snowd, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Liquid water equivalent snow depth", & + units="mm" & + ) + + call cable_restart_variable_write_darray( & + var_name="snage", & + var_dims=["mp"], & + data=ssnow%snage, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Snow age", & + units="??" & + ) + + call cable_restart_variable_write_darray( & + var_name="smass", & + var_dims=["mp","snow"], & + data=ssnow%smass, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Average layer snow mass", & + units="kg/m^2" & + ) + + call cable_restart_variable_write_darray( & + var_name="sdepth", & + var_dims=["mp", "snow"], & + data=ssnow%sdepth, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Snow layer depth", & + units="m" & + ) + + call cable_restart_variable_write_darray( & + var_name="osnowd", & + var_dims=["mp"], & + data=ssnow%osnowd, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Previous time step snow depth in water equivalent", & + units="mm" & + ) + + call cable_restart_variable_write_darray( & + var_name="isflag", & + var_dims=["mp"], & + data=ssnow%isflag, & + var_type=CABLE_NETCDF_INT, & + long_name="Snow layer scheme flag", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="cansto", & + var_dims=["mp"], & + data=canopy%cansto, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Canopy surface water storage", & + units="mm" & + ) + + call cable_restart_variable_write_darray( & + var_name="ghflux", & + var_dims=["mp"], & + data=canopy%ghflux, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="????", & + units="W/m^2?" & + ) + + call cable_restart_variable_write_darray( & + var_name="sghflux", & + var_dims=["mp"], & + data=canopy%sghflux, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="????", & + units="W/m^2?" & + ) + + call cable_restart_variable_write_darray( & + var_name="ga", & + var_dims=["mp"], & + data=canopy%ga, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Ground heat flux", & + units="W/m^2" & + ) + + call cable_restart_variable_write_darray( & + var_name="dgdtg", & + var_dims=["mp"], & + data=canopy%dgdtg, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Derivative of ground heat flux wrt soil temperature", & + units="W/m^2/K" & + ) + + call cable_restart_variable_write_darray( & + var_name="fev", & + var_dims=["mp"], & + data=canopy%fev, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Latent heat flux from vegetation", & + units="W/m^2" & + ) + + call cable_restart_variable_write_darray( & + var_name="fes", & + var_dims=["mp"], & + data=canopy%fes, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Latent heat flux from soil", & + units="W/m^2" & + ) + + call cable_restart_variable_write_darray( & + var_name="fhs", & + var_dims=["mp"], & + data=canopy%fhs, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Sensible heat flux from soil", & + units="W/m^2" & + ) + + call cable_restart_variable_write_darray( & + var_name="cplant", & + var_dims=["mp", "plant_carbon_pools"], & + data=bgc%cplant, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Plant carbon stores", & + units="gC/m^2" & + ) + + call cable_restart_variable_write_darray( & + var_name="csoil", & + var_dims=["mp", "soil_carbon_pools"], & + data=bgc%csoil, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Soil carbon stores", & + units="gC/m^2" & + ) + + call cable_restart_variable_write_darray( & + var_name="wbtot0", & + var_dims=["mp"], & + data=bal%wbtot0, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Initial time step soil water total", & + units="mm" & + ) + + call cable_restart_variable_write_darray( & + var_name="osnowd0", & + var_dims=["mp"], & + data=bal%osnowd0, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Initial time step snow water total", & + units="mm" & + ) + + call cable_restart_variable_write_darray( & + var_name="albedo", & + var_dims=["mp","rad"], & + data=rad%albedo, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Albedo for shortwave and NIR radiation", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="trad", & + var_dims=["mp"], & + data=rad%trad, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Surface radiative temperature (soil/snow/veg inclusive)", & + units="K" & + ) + + call cable_restart_variable_write_darray( & + var_name="iveg", & + var_dims=["mp"], & + data=veg%iveg, & + var_type=CABLE_NETCDF_INT, & + long_name="Vegetation type", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="isoil", & + var_dims=["mp"], & + data=soil%isoilm, & + var_type=CABLE_NETCDF_INT, & + long_name="Soil type", & + units="-" & + ) + + call cable_restart_variable_write( & + var_name="zse", & + var_dims=["soil"], & + data=soil%zse, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Depth of each soil layer", & + units="m" & + ) + + call cable_restart_variable_write_darray( & + var_name="albsoil", & + var_dims=["mp","rad"], & + data=soil%albsoil, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Soil reflectance", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="GWwb", & + var_dims=["mp"], & + data=ssnow%GWwb, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Groundwater water content", & + units="mm3/mm3" & + ) + + if (cable_user%soil_struc == 'sli' .or. cable_user%fwsoil_switch == 'Haverd2013') then + + call cable_restart_variable_write_darray( & + var_name="gamma", & + var_dims=["mp"], & + data=veg%gamma, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Parameter in root efficiency function (Lai and Katul 2000)", & + units="-" & + ) + + end if + + if (cable_user%soil_struc == 'sli') then + + call cable_restart_variable_write_darray( & + var_name="nhorizons", & + var_dims=["mp"], & + data=soil%nhorizons, & + var_type=CABLE_NETCDF_INT, & + long_name="Number of soil horizons", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="ishorizon", & + var_dims=["mp"], & + data=soil%ishorizon, & + var_type=CABLE_NETCDF_INT, & + long_name="Horizon number", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="clitt", & + var_dims=["mp"], & + data=veg%clitt, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Litter layer carbon content", & + units="tC/ha" & + ) + + call cable_restart_variable_write_darray( & + var_name="ZR", & + var_dims=["mp"], & + data=veg%ZR, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Maximum rooting depth", & + units="cm" & + ) + + call cable_restart_variable_write_darray( & + var_name="F10", & + var_dims=["mp"], & + data=veg%F10, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Fraction of roots in top 10 cm", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="S", & + var_dims=["mp"], & + data=ssnow%S, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Fractional soil moisture content relative to saturated value", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="Tsoil", & + var_dims=["mp"], & + data=ssnow%Tsoil, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Soil temperature", & + units="degC" & + ) + + call cable_restart_variable_write_darray( & + var_name="snowliq", & + var_dims=["mp"], & + data=ssnow%snowliq, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Liquid water content of snowpack", & + units="mm" & + ) + + call cable_restart_variable_write_darray( & + var_name="sconds", & + var_dims=["mp"], & + data=ssnow%sconds, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Thermal conductivity of snowpack", & + units="W/m/K" & + ) + + call cable_restart_variable_write_darray( & + var_name="h0", & + var_dims=["mp"], & + data=ssnow%h0, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Pond height above soil", & + units="m" & + ) + + call cable_restart_variable_write_darray( & + var_name="nsnow", & + var_dims=["mp"], & + data=ssnow%nsnow, & + var_type=CABLE_NETCDF_INT, & + long_name="Number of snow layers", & + units="-" & + ) + + call cable_restart_variable_write_darray( & + var_name="Tsurface", & + var_dims=["mp"], & + data=ssnow%Tsurface, & + var_type=CABLE_NETCDF_FLOAT, & + long_name="Soil or snow surface temperature", & + units="degC" & + ) + + end if + + end subroutine cable_restart_write + +end module diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 14aa29449..15ca60e40 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,13 +113,22 @@ 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_commit + use cable_output_prototype_v2_mod, only: cable_output_update + use cable_output_prototype_v2_mod, only: cable_output_write_parameters + use cable_output_definitions_mod, only: cable_output_definitions_set + use cable_netcdf_mod, only: cable_netcdf_mod_init, cable_netcdf_mod_end + use cable_restart_mod, only: cable_restart_mod_init, cable_restart_mod_end + use cable_restart_write_mod, only: cable_restart_write USE cable_checks_module, ONLY: constant_check_range USE cable_write_module, ONLY: nullify_write USE cable_IO_vars_module, ONLY: timeunits,calendar USE cable_cbm_module, ONLY : cbm !mpidiff USE cable_climate_mod - + ! modules related to CASA-CNP USE casadimension, ONLY: icycle USE casavariable, ONLY: casafile, casa_biome, casa_pool, casa_flux, & @@ -273,10 +283,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 +475,14 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi call cable_io_decomp_init(io_decomp) + if (output%restart) call cable_restart_mod_init() + + if (.not. casaonly) then + call cable_output_mod_init() + call cable_output_definitions_set(io_decomp, canopy, soil) + call cable_output_commit() + end if + ENDIF ! CALL 1 ! globally (WRT code) accessible kend through USE cable_common_module @@ -572,6 +593,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 +613,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 +634,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 +760,16 @@ 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( & + time_index=ktau, & + dels=dels, & + leaps=leaps, & + start_year=start_year, & + patch=patch, & + landpt=landpt, & + met=met & + ) ENDIF @@ -926,6 +974,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi IF ( SpinConv .AND. .NOT. CASAONLY ) THEN ! Close output file and deallocate main variables: + call cable_output_mod_end() CALL close_output_file( bal, air, bgc, canopy, met, & rad, rough, soil, ssnow, & sum_flux, veg ) @@ -1004,9 +1053,12 @@ 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_restart_write(dels * ktau, soil, veg, ssnow, canopy, rough, rad, bgc, bal, met) + CALL cable_restart_mod_end() + END IF !mpidiff IF (cable_user%CALL_climate) & CALL WRITE_CLIMATE_RESTART_NC ( climate, ktauday ) @@ -1014,6 +1066,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site, mpi !--- LN ------------------------------------------[ ENDIF + call cable_netcdf_mod_end() IF ( TRIM(cable_user%MetType) .NE. "gswp" .AND. & diff --git a/src/util/aggregator.F90 b/src/util/aggregator.F90 new file mode 100644 index 000000000..f52e9da1a --- /dev/null +++ b/src/util/aggregator.F90 @@ -0,0 +1,174 @@ +module aggregator_mod + + use iso_fortran_env, only: int32, real32, real64 + use cable_abort_module, only: cable_abort + + use aggregator_types_mod + + implicit none + private + + public :: aggregator_t + public :: aggregator_handle_t + public :: aggregator_int32_1d_t + public :: aggregator_int32_2d_t + public :: aggregator_int32_3d_t + public :: aggregator_real32_1d_t + public :: aggregator_real32_2d_t + public :: aggregator_real32_3d_t + public :: aggregator_real64_1d_t + public :: aggregator_real64_2d_t + public :: aggregator_real64_3d_t + public :: aggregator_mod_init + public :: aggregator_mod_end + public :: new_aggregator + public :: store_aggregator + + type aggregator_store_t + class(aggregator_t), allocatable :: aggregator + end type aggregator_store_t + + interface new_aggregator + 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_1d + module procedure new_aggregator_real32_2d + module procedure new_aggregator_real32_3d + module procedure new_aggregator_real64_1d + module procedure new_aggregator_real64_2d + module procedure new_aggregator_real64_3d + end interface + + integer, parameter :: DEFAULT_MAX_AGGREGATORS = 1000 + integer :: num_aggregators = 0 + type(aggregator_store_t), allocatable, target :: aggregators(:) + + +contains + + subroutine aggregator_mod_init() + integer :: max_aggregators + max_aggregators = DEFAULT_MAX_AGGREGATORS ! TODO(Sean): Make this configurable + allocate(aggregators(max_aggregators)) + num_aggregators = 0 + end subroutine + + subroutine aggregator_mod_end() + if (allocated(aggregators)) deallocate(aggregators) + end subroutine + + function store_aggregator(aggregator) result(aggregator_handle) + class(aggregator_t), intent(in) :: aggregator + type(aggregator_handle_t) :: aggregator_handle + integer :: index + + num_aggregators = num_aggregators + 1 + if (num_aggregators > size(aggregators)) then + ! Note: we cannot resize the aggregators array as this would + ! invalidate any pointers to its elements elsewhere. + ! TODO(Sean): provide a recommendation to increase the max number of aggregators + call cable_abort("Exceeded maximum number of aggregators.", __FILE__, __LINE__) + end if + + index = num_aggregators + + ! Copy the aggregator into the store + aggregators(index)%aggregator = aggregator + + ! Create a handle pointing to the stored aggregator + aggregator_handle%aggregator => aggregators(index)%aggregator + + end function store_aggregator + + function new_aggregator_int32_1d_t(source_data, method) result(agg) + integer(kind=int32), dimension(:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_int32_1d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_int32_1d_t + + function new_aggregator_int32_2d_t(source_data, method) result(agg) + integer(kind=int32), dimension(:,:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_int32_2d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_int32_2d_t + + function new_aggregator_int32_3d_t(source_data, method) result(agg) + integer(kind=int32), dimension(:,:,:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_int32_3d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_int32_3d_t + + function new_aggregator_real32_1d(source_data, method) result(agg) + real(kind=real32), dimension(:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_real32_1d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_real32_1d + + function new_aggregator_real32_2d(source_data, method) result(agg) + real(kind=real32), dimension(:,:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_real32_2d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_real32_2d + + function new_aggregator_real32_3d(source_data, method) result(agg) + real(kind=real32), dimension(:,:,:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_real32_3d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_real32_3d + + function new_aggregator_real64_1d(source_data, method) result(agg) + real(kind=real64), dimension(:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_real64_1d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_real64_1d + + function new_aggregator_real64_2d(source_data, method) result(agg) + real(kind=real64), dimension(:,:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_real64_2d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_real64_2d + + function new_aggregator_real64_3d(source_data, method) result(agg) + real(kind=real64), dimension(:,:,:), intent(inout), target :: source_data + character(len=*), intent(in) :: method + type(aggregator_real64_3d_t) :: agg + + agg%source_data => source_data + call agg%set_method(method) + + end function new_aggregator_real64_3d + +end module diff --git a/src/util/aggregator_types.F90 b/src/util/aggregator_types.F90 new file mode 100644 index 000000000..009174e21 --- /dev/null +++ b/src/util/aggregator_types.F90 @@ -0,0 +1,390 @@ +module aggregator_types_mod + use iso_fortran_env, only: int32, real32, real64 + use cable_abort_module, only: cable_abort + implicit none + private + + public :: aggregator_t + public :: aggregator_handle_t + public :: aggregator_int32_1d_t + public :: aggregator_int32_2d_t + public :: aggregator_int32_3d_t + public :: aggregator_real32_1d_t + public :: aggregator_real32_2d_t + public :: aggregator_real32_3d_t + public :: aggregator_real64_1d_t + public :: aggregator_real64_2d_t + public :: aggregator_real64_3d_t + + 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 + 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 aggregator_handle_t + class(aggregator_t), pointer :: aggregator => null() + contains + procedure :: init => aggregator_handle_init + procedure :: accumulate => aggregator_handle_accumulate + procedure :: reset => aggregator_handle_reset + end type aggregator_handle_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_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_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 + +contains + + subroutine aggregator_handle_init(this) + class(aggregator_handle_t), intent(inout) :: this + + call this%aggregator%init() + + end subroutine aggregator_handle_init + + subroutine aggregator_handle_accumulate(this) + class(aggregator_handle_t), intent(inout) :: this + + call this%aggregator%accumulate() + + end subroutine aggregator_handle_accumulate + + subroutine aggregator_handle_reset(this) + class(aggregator_handle_t), intent(inout) :: this + + call this%aggregator%reset() + + end subroutine aggregator_handle_reset + + subroutine aggregator_init(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + 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_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_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%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 + + subroutine mean_accumulate(this) + class(aggregator_t), intent(inout) :: this + + select type (this) + 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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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 + +end module diff --git a/src/util/cable_grid_reductions.F90 b/src/util/cable_grid_reductions.F90 new file mode 100644 index 000000000..7796cd7ee --- /dev/null +++ b/src/util/cable_grid_reductions.F90 @@ -0,0 +1,143 @@ +module cable_grid_reductions_mod + + use iso_fortran_env, only: real32, real64 + + use cable_io_vars_module, only: patch_type, land_type + + implicit none + private + + public :: grid_cell_average + + 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 + +end module diff --git a/src/util/cable_timing_utils.F90 b/src/util/cable_timing_utils.F90 new file mode 100644 index 000000000..a3c0353cd --- /dev/null +++ b/src/util/cable_timing_utils.F90 @@ -0,0 +1,60 @@ +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_day + integer :: last_day_of_month_in_accumulated_days(months_in_year) ! TODO(Sean): better variable name? + + select case (frequency) + ! TODO(Sean): implement case for custom hourly frequencies + case ('all') + match = .true. + case ('daily') + time_steps_per_day = seconds_per_hour * hours_per_day / int(dels) + match = mod(ktau, time_steps_per_day) == 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/netcdf/cable_netcdf.F90 b/src/util/netcdf/cable_netcdf.F90 index 086c9c307..7246f23e3 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 @@ -150,6 +173,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 +190,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) @@ -607,16 +635,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 +666,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..f83a38049 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 @@ -100,16 +101,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 +135,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 +151,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 diff --git a/src/util/netcdf/nf90/cable_netcdf_nf90.F90 b/src/util/netcdf/nf90/cable_netcdf_nf90.F90 index dfa8a08d4..1bf6bfcf0 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 @@ -136,6 +143,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 +186,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 +226,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 +252,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) diff --git a/src/util/netcdf/pio/cable_netcdf_pio.F90 b/src/util/netcdf/pio/cable_netcdf_pio.F90 index 3a2cacb20..3c8c807b0 100644 --- a/src/util/netcdf/pio/cable_netcdf_pio.F90 +++ b/src/util/netcdf/pio/cable_netcdf_pio.F90 @@ -25,6 +25,7 @@ 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_inq_varid @@ -35,8 +36,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 +76,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 @@ -151,6 +158,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 +253,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 +301,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 +327,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)