From 7d53910c8aa146667018d523eed2628c20527ebd Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 12 Dec 2025 11:39:41 +1100 Subject: [PATCH 01/28] Remove unused file --- .../explicit/cable_explicit_main.F90.cnp | 401 ------------------ 1 file changed, 401 deletions(-) delete mode 100644 src/coupled/AM3/control/cable/interface/explicit/cable_explicit_main.F90.cnp diff --git a/src/coupled/AM3/control/cable/interface/explicit/cable_explicit_main.F90.cnp b/src/coupled/AM3/control/cable/interface/explicit/cable_explicit_main.F90.cnp deleted file mode 100644 index f5e4263f6..000000000 --- a/src/coupled/AM3/control/cable/interface/explicit/cable_explicit_main.F90.cnp +++ /dev/null @@ -1,401 +0,0 @@ -MODULE cable_explicit_main_mod - -CONTAINS - -SUBROUTINE cable_explicit_main( & - ! IN: UM/JULES model/grid parameters, fields, mappings - mype, timestep_len, timestep_number, row_length, rows, land_pts, & - nsurft, npft, sm_levels, dzsoil, land_index, surft_pts, surft_index, & - cosine_zenith_angle, latitude, longitude, Fland, tile_frac, & - - ! IN: soil parameters !1 is only allowable index in UM - bexp, hcon, satcon, sathh, smvcst, smvcwt, smvccl, albsoil, & - - ! IN: Met forcing: - lw_down, sw_surft, ls_rain, ls_snow, & - tl_1, qw_1, vshr_land, pstar, z1_tq, z1_uv, canopy_tile, & - ! This an outlier IN here. INOUT @ implicit. (was)OUT at extras - ! I think we are dealing with it OK now but confusion could be removed - snow_tile, & - - ! IN: canopy height, LAI seasonally presecribed, potentially prognostic - ! IN: CO2 mass mixing ratio - canht_pft, lai_pft, CO2_MMR, & - - ! TYPEs passed from top_level to maintain scope, access to UM STASH - ! ------------------------------------------------------------------ - ! IN: tiled soil/snow prognostics - IN here. INOUT @ implicit - ! INOUT: Carries fields needed by CABLE b/n pathways (rad, explicit etc) - ! Currently carrying CABLE TYPEs (canopy%, rad% etc). - ! IN: pars carries vegin/soilin - potentially redundant - progs, work, pars, & - - ! IN: CASA-CNP prognostics - IN here. INOUT @ implicit - progs_cnp, & - - ! OUT: UM fields UNPACKed from CABLE (@ explicit) - ftl_tile, fqw_tile, tstar_tile, dtstar_surft, & - u_s, u_s_std_tile, cd_tile, ch_tile, & - radnet_tile, fraca, resfs, resft, z0h_tile, z0m_tile, & - recip_l_mo_tile, epot_tile, npp_pft_acc, resp_w_pft_acc ) - -! subrs -USE cable_explicit_driv_mod, ONLY: cable_explicit_driver -USE cable_expl_unpack_mod, ONLY: cable_expl_unpack -USE init_active_tile_mask_mod, ONLY: init_active_tile_mask_cbl - -! data: TYPE definitions of passed asarguments -USE progs_cbl_vars_mod, ONLY: progs_cbl_vars_type ! CABLE requires extra progs -USE work_vars_mod_cbl, ONLY: work_vars_type ! and some kept thru timestep -USE params_io_mod_cbl, ONLY: params_io_data_type -USE params_io_mod_cbl, ONLY: params_io_type -USE progs_cnp_vars_mod, ONLY: progs_cnp_vars_type ! CASA requires progs - -! data: Scalars -USE grid_constants_mod_cbl, ONLY: nrb, nrs, mp -USE cable_common_module, ONLY: knode_gl, ktau_gl, kwidth_gl, cable_runtime, & - cable_user, redistrb, satuparam, wiltparam, & - l_casacnp -USE casadimension, ONLY: icycle ! 0=No CASA- [1=C,2=CN,3=CNP] -!Leave for reference -!! USE atm_fields_real_mod, ONLY : -!! NPP_PFT_ACC, RSP_W_PFT_ACC, & -!! aquifer_moist_cable,aquifer_thickness_cable, & -!! slope_avg_cable,slope_std_cable,& -!! visc_sublayer_depth,aquifer_perm_cable,& -!! aquifer_draindens_cable - -IMPLICIT NONE - -INTEGER, INTENT(IN) :: mype ! # processor -REAL, INTENT(IN) :: timestep_len ! # seconds (cucurrently 1200) -INTEGER, INTENT(IN) :: timestep_number ! # timestep (cucurrently 3 per hr) -INTEGER, INTENT(IN) :: row_length ! # columns in spatial grid -INTEGER, INTENT(IN) :: rows ! # rows in spatial grid -INTEGER, INTENT(IN) :: land_pts ! # land points being processed -INTEGER, INTENT(IN) :: nsurft ! # tiles -INTEGER, INTENT(IN) :: npft ! # plant functional types -INTEGER, INTENT(IN) :: sm_levels ! # soil layers -REAL, INTENT(IN) :: dzsoil(sm_levels) ! soil layer thicknesses - -INTEGER, INTENT(IN) :: land_index(land_pts) ! land point indices - ! recipe back to (i,j) cell -INTEGER, INTENT(IN) :: surft_pts(nsurft) ! # land points per tile -INTEGER, INTENT(IN) :: surft_index(land_pts, nsurft) ! tile points indices - ! recipe back to land_index - -REAL, INTENT(IN) :: canht_pft(land_pts, npft) ! canopy height (seasonal) -REAL, INTENT(IN) :: lai_pft(land_pts, npft) ! LAI (seasonal) -REAL, INTENT(IN) :: fland(land_pts) ! land fraction (<1 for coastal) -REAL, INTENT(IN) :: co2_mmr ! prescribed MMR -REAL, INTENT(IN) :: tile_frac(land_pts, nsurft) ! tile fraction - -REAL, INTENT(IN) :: cosine_zenith_angle(row_length,rows) -REAL, INTENT(IN) :: latitude (row_length,rows) -REAL, INTENT(IN) :: longitude(row_length,rows) - -! soil parameters -REAL, INTENT(IN) :: bexp (land_pts, sm_levels) ! parameter b in Campbell eqn -REAL, INTENT(IN) :: sathh(land_pts, sm_levels) ! -REAL, INTENT(IN) :: smvcst(land_pts, sm_levels) ! -REAL, INTENT(IN) :: smvcwt(land_pts, sm_levels) ! -REAL, INTENT(IN) :: smvccl(land_pts, sm_levels) ! -REAL, INTENT(IN) :: hcon(land_pts) ! Soil thermal conductivity (W/m/K). -REAL, INTENT(IN) :: albsoil(land_pts) ! bare soil albedo -REAL, INTENT(IN) :: satcon(land_pts, sm_levels) ! hydraulic conductivity - ! @ saturation [mm/s] - -! "forcing" -REAL, INTENT(IN) :: lw_down(row_length,rows) -REAL, INTENT(IN) :: ls_rain(row_length,rows) -REAL, INTENT(IN) :: ls_snow(row_length,rows) -REAL, INTENT(IN) :: sw_surft(land_pts, nsurft) -REAL, INTENT(IN) :: tl_1(row_length,rows) -REAL, INTENT(IN) :: qw_1(row_length,rows) -REAL, INTENT(IN) :: vshr_land(row_length,rows) -REAL, INTENT(IN) :: pstar(row_length,rows) -REAL, INTENT(IN) :: z1_tq(row_length,rows) -REAL, INTENT(IN) :: z1_uv(row_length,rows) - -! prognostics -REAL, INTENT(IN) :: canopy_tile(land_pts, nsurft) -REAL, INTENT(IN) :: snow_tile(land_pts, nsurft) ! Lying snow [kg/m2] - -! TYPEs passed from top_level to maintain scope, access to UM STASH -! IN: tiled soil/snow prognostics - IN here. INOUT @ implicit -! INOUT: Carries fields needed by CABLE b/n pathways (rad, explicit etc) -! Currently carrying CABLE TYPEs (canopy%, rad% etc). -! IN: pars carries vegin/soilin - potentially redundant -TYPE(progs_cbl_vars_type), INTENT(IN) :: progs -TYPE(params_io_data_type), INTENT(IN) :: pars -TYPE(work_vars_type), INTENT(INOUT) :: work -! IN: CASA-CNP prognostics - IN here. INOUT @ implicit -TYPE(progs_cnp_vars_type), INTENT(IN) :: progs_cnp - -! OUT: UM fields UNPACKed from CABLE (@ explicit) -REAL, INTENT(OUT) :: ftl_tile(land_pts,nsurft) ! surface FTL for land tiles -REAL, INTENT(OUT) :: fqw_tile(land_pts,nsurft) ! surface FQW for land tiles -REAL, INTENT(OUT) :: tstar_tile(land_pts,nsurft) ! radiative surf. temperature -REAL, INTENT(OUT) :: dtstar_surft(land_pts,nsurft) ! -REAL, INTENT(OUT) :: u_s(row_length,rows) ! friction velocity (m/s) -REAL, INTENT(OUT) :: u_s_std_tile(land_pts,nsurft)! -REAL, INTENT(OUT) :: cd_tile(land_pts,nsurft) ! Drag coefficient -REAL, INTENT(OUT) :: ch_tile(land_pts,nsurft) ! Transfer coefficient -REAL, INTENT(OUT) :: radnet_tile(land_pts,nsurft) ! Surface net radiation -REAL, INTENT(OUT) :: z0h_tile(land_pts,nsurft) ! roughness -REAL, INTENT(OUT) :: z0m_tile(land_pts,nsurft) ! roughness -REAL, INTENT(OUT) :: epot_tile(land_pts,nsurft) ! -REAL, INTENT(OUT) :: recip_l_mo_tile(land_pts,nsurft) ! Reciprocal:Monin-Obukhov - ! length for tiles (m^-1) -REAL, INTENT(OUT) :: fraca(land_pts,nsurft) ! Fraction - surface moisture -REAL, INTENT(OUT) :: RESFS(land_pts,nsurft) - ! Combined soil, stomatal & aerodynamic resistance - ! factor for fraction (1-FRACA) of snow-free land tiles -REAL, INTENT(OUT) :: RESFT(land_pts,nsurft) - ! Total resistance factor. - ! FRACA+(1-FRACA)*RESFS for snow-free l_tile_pts, - ! 1 for snow. -REAL, INTENT(IN) :: npp_pft_acc(land_pts,npft) -REAL, INTENT(IN) :: resp_w_pft_acc (land_pts,npft) - -!___ local vars, may be passed as args downstream -LOGICAL :: cbl_standalone = .FALSE. !needs to be set from namelist -LOGICAL :: jls_standalone = .FALSE. !needs to be set from namelist -LOGICAL :: jls_radiation = .FALSE. !needs to be set from n amelist - -INTEGER :: isnow_flg_cable(land_pts, nsurft) -REAL :: radians_degrees -REAL :: latitude_deg(row_length,rows) -REAL :: longitude_deg(row_length,rows) -REAL :: sw_down_ij(row_length,rows,nrs) -REAL :: sw_down_TOT(row_length,rows) -REAL :: sw_down_DIR(row_length,rows) -REAL :: sw_down_VIS(row_length,rows) -REAL :: sw_down_NIR(row_length,rows) -REAL :: beamFrac_VIS(row_length,rows) -REAL :: beamFrac_NIR(row_length,rows) -REAL :: beamFrac_TOT(row_length,rows) - -LOGICAL, ALLOCATABLE :: l_tile_pts(:,:) - -INTEGER :: i,j,k,l,n -LOGICAL, SAVE :: zero_points_warning = .TRUE. - -CHARACTER(LEN=*), PARAMETER :: subr_name = "cable_explicit_main" - -IF( land_pts == 0 ) THEN - IF( zero_points_warning ) THEN - WRITE(6,*) "Reached CABLE ", subr_name, & - " even though zero land_points on processor ", mype - END IF - zero_points_warning = .FALSE. - RETURN -END IF - -!--- Set up some cable-globals ------------------------------------------------- -cable_runtime%um = .TRUE. -cable_runtime%um_explicit = .TRUE. - -! this done every call (maybe we hould pass this through work%) -!------------------------------------------------------------------------------ -! Determine the number of active tiles -mp = SUM(surft_pts) - -IF( .NOT. ALLOCATED(l_tile_pts) ) ALLOCATE( l_tile_pts(land_pts, nsurft) ) - -! Define mapping mask. i.e. l_tile_pts =TRUE (active) , where tile_frac > 0 -CALL init_active_tile_mask_cbl( l_tile_pts, land_pts, nsurft, tile_frac ) -!------------------------------------------------------------------------------- - -!!extracted from ap/um/rose-app.conf au-aa809@2729 -![namelist:cable] -cable_user%diag_soil_resp='ON' -cable_user%fwsoil_switch='Haverd2013' -cable_user%gs_switch='medlyn' -cable_user%gw_model=.false. -cable_user%l_rev_corr=.true. -cable_user%l_revised_coupling=.true. -cable_user%or_evap=.false. -!cable_user%soil_thermal_fix=.true. -cable_user%soil_thermal_fix=.false.! fudge - worked to dt=4 -cable_user%ssnow_potev='HDM' -icycle=3 !previously blocked?? -l_casacnp=.true. !previously blocked?? -redistrb=.false. -satuparam=0.8 -wiltparam=0.5 - -! initialize processor number, timestep len -knode_gl = mype -kwidth_gl = INT(timestep_len) -ktau_gl = timestep_number - -!--- Convert lat/long to degrees -radians_degrees = 180.0 / ( 4.0*atan(1.0) ) ! 180 / PI -latitude_deg = latitude * radians_degrees -longitude_deg = longitude * radians_degrees - -isnow_flg_cable = INT(progs%ThreeLayerSnowFlag_CABLE) - -!--- Fix SW for CABLE ---------------------------------------------------------------------------- -sw_down_ij(:,:,:) = 0.0 -sw_down_TOT(:,:) = 0.0 -sw_down_DIR(:,:) = 0.0 -sw_down_VIS(:,:) = 0.0 -sw_down_NIR(:,:) = 0.0 -beamFrac_VIS(:,:) = 0.0 -beamFrac_NIR(:,:) = 0.0 -beamFrac_TOT(:,:) = 0.0 - -IF(jls_standalone) THEN - - DO n = 1, nsurft - ! loop over number of points per tile - DO k = 1, surft_pts(n) - l = surft_index(k, n) - j = (land_index(l) - 1) / row_length + 1 - i = land_index(l) - (j-1) * row_length - sw_down_VIS(i, j) = sw_surft(l,n) / 2.0 - END DO - END DO - sw_down_NIR(:,:) = sw_down_VIS(:,:) - -ELSE - - ! in all cases zenith angle needs to be applied - DO n = 1, nrs - sw_down_ij(:,:,n) = work%sw_down_ij(:,:,n) * cosine_zenith_angle(:,:) - END DO - - ! SUM over ALL components of sw_down_ij - sw_down_TOT(:,:) = sw_down_ij(:,:,1) + sw_down_ij(:,:,2) + & - sw_down_ij(:,:,3) + sw_down_ij(:,:,4) - - ! SUM DIRect components of sw_down_ij(in VIS & NIR ) - sw_down_DIR(:,:) = sw_down_ij(:,:,1) + sw_down_ij(:,:,3) - - ! SUM VIS components of sw_down_ij(incl DIR & DIF ) - sw_down_VIS(:,:) = sw_down_ij(:,:,1) + sw_down_ij(:,:,2) - - ! SUM NIR components of sw_down_ij(incl DIR & DIF ) - sw_down_NIR(:,:) = sw_down_ij(:,:,3) + sw_down_ij(:,:,4) - - ! beam(DIR) fraction in VISible spectrum - beamFrac_VIS(:,:) = sw_down_ij(:,:,1) / MAX( 0.1, sw_down_VIS(:,:) ) - - ! beam(DIR) fraction in NIR spectrum - beamFrac_NIR(:,:) = sw_down_ij(:,:,3) / MAX( 0.1, sw_down_NIR(:,:) ) - - ! beam(DIR) fraction for all solar - beamFrac_TOT(:,:) = sw_down_DIR(:,:) / MAX( 0.1, sw_down_TOT(:,:) ) - -ENDIF - -!---------------------------------------------------------------------------- -!--- CALL _driver to run specific and necessary components of CABLE with IN - -!--- args PACKED to force CABLE -!------------------------------------------------------------------------------- -CALL cable_explicit_driver( & - ! IN: UM/JULES/CABLE model/grid parameters, fields, mappings - mype, row_length, rows, land_pts, nsurft, npft, sm_levels, dzsoil, & - timestep_len, timestep_number, mp, nrb, land_index, surft_pts, surft_index, & - l_tile_pts, latitude_deg, longitude_deg, cosine_zenith_angle, Fland, & - tile_frac, L_casacnp, & - - ! IN: soil parameters !1 is only allowable index in UM - bexp, hcon, satcon, sathh, smvcst, smvcwt, smvccl, albsoil, & - - ! IN: SW forcing: manipulated for CABLE - sw_down_VIS, sw_down_NIR, beamFrac_VIS, beamFrac_NIR, beamFrac_TOT, & - - ! IN: Met forcing: - lw_down, ls_rain, ls_snow, & - tl_1, qw_1, vshr_land, pstar, z1_tq, z1_uv, canopy_tile, & - ! This an outlier IN here. INOUT @ implicit. (was)OUT at extras - ! I think we are dealing with it OK now but confusion could be removed - snow_tile, & - - ! IN: canopy height, LAI seasonally presecribed, potentially prognostic - ! IN: CO2 mass mixing ratio - canht_pft, lai_pft, CO2_MMR, & - - ! IN: carries vegin/soilin - potentially redundant - pars, & - - ! IN: tiled soil/snow prognostics - IN here. INOUT @ implicit - progs%soiltemp_CABLE, progs%soilmoisture_CABLE, progs%FrozenSoilFrac_CABLE, & - isnow_flg_cable, progs%SnowDepth_CABLE, progs%SnowMass_CABLE, & - progs%SnowTemp_CABLE, progs%SnowDensity_CABLE, progs%snowage_CABLE, & - progs%snowosurft, progs%OneLyrSnowDensity_CABLE, & - - ! IN: casa-CNP prognostics - IN here. INOUT @ implicit - progs_cnp% C_pool_casa, progs_cnp% N_pool_casa, progs_cnp% P_pool_casa, & - progs_cnp% soil_order_casa, & - progs_cnp% N_dep_casa, progs_cnp% N_fix_casa, & - progs_cnp% P_dust_casa, progs_cnp% P_weath_casa, & - progs_cnp% LAI_casa, progs_cnp% phenphase_casa, & - progs_cnp% wood_flux_C, progs_cnp% wood_flux_N, progs_cnp% wood_flux_P, & - progs_cnp% wood_hvest_C, progs_cnp% wood_hvest_N, progs_cnp% wood_hvest_P, & - progs_cnp% thinning, & - - ! INOUT: CABLE TYPEs roughly grouped fields per module - work%rad, work%met, work%rough, work%canopy, work%veg, work%soil, & - work%ssnow, work%bal, work%air, work%bgc, work%sum_flux, & - - ! INOUT: CASA TYPEs roughly grouped fields per module - work%casapool, work%casaflux, & - work%sum_casapool, work%sum_casaflux, work%casabiome, & - work%casamet, work%casabal, work%phen, & - - !IN: persistent veg%iveg, soil%isoilm are initialized on first rad/alb call - work%veg%iveg, work%soil%isoilm, & - !OUT: currently being passed back to UM in veg%hc, veg%vlai - work%veg%hc, work%veg%vlai, & - - !IN: currently being passed from prev radiation call through work% - ! jhan:quirky, snow (in turn reduced LAI due to snow) can evolve through a - ! constant rad dt. However reducedLAIdue2snow used ubiquitously as trigger - ! Further, snow does NOT evolve in explicit AND reducedLAIdue2snow absent - ! in implicit - work%reducedLAIdue2snow, & - - !GW - !visc_sublayer_depth, smgw_tile, slope_avg, slope_std, - !dz_gw, perm_gw, drain_gw, - !casa progs - !CPOOL_TILE, NPOOL_TILE, PPOOL_TILE, SOIL_ORDER, NIDEP, - !NIFIX, PWEA, PDUST, GLAI, PHENPHASE, - - !IN: if not passed a dangling argument would ensue - npp_pft_acc, resp_w_pft_acc ) - -!---------------------------------------------------------------------------- -!--- CALL _unpack to unpack variables from CABLE back to UM format to return -!---------------------------------------------------------------------------- -call cable_expl_unpack( & - ! IN: UM/JULES/CABLE model/grid parameters, fields, mappings - row_length, rows, land_pts, nsurft, npft, mp, land_index, surft_pts, & - surft_index, l_tile_pts, fland, tile_frac, latitude, longitude, & - - !OUT: UM fields to be updated - ftl_tile, fqw_tile, tstar_tile, dtstar_surft , u_s, u_s_std_tile, cd_tile, & - ch_tile, radnet_tile, fraca, resfs, resft, z0h_tile, z0m_tile, & - recip_l_mo_tile, epot_tile, & - - !IN: UM fields to be updated FROM these CABLE fields - work%canopy%fh, work%canopy%fes, work%canopy%fev, work%canopy%us, & - work%canopy%cdtq,work%canopy%fwet, work%canopy%wetfac_cs, & - work%canopy%rnet, work%canopy%zetar, work%canopy%epot, work%rad%trad, & - work%rad%otrad, work%rad%transd, work%rough%z0m, work%rough%zref_tq, & - - !IN: UM fields used in derivation of fields to be updated - work%ssnow%snowd, work%ssnow%cls, work%air%rlam, work%air%rho, work%met%ua ) - -cable_runtime%um_explicit = .FALSE. - -RETURN - -END SUBROUTINE cable_explicit_main - -END MODULE cable_explicit_main_mod - From b042689066d4b2a5629a153e435911240d85166c Mon Sep 17 00:00:00 2001 From: Whyborn Date: Tue, 16 Dec 2025 15:24:08 +1100 Subject: [PATCH 02/28] Start substitution of datetime with initialisation --- src/offline/cable_serial.F90 | 67 +++++++++++++++++++++++------------- 1 file changed, 43 insertions(+), 24 deletions(-) diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 366aa58de..ab421bb6d 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -268,6 +268,33 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) integer, dimension(:), allocatable, save :: cstart,cend,nap real(r_2), dimension(:,:,:), allocatable, save :: patchfrac_new + type(datetime) :: ts_start, ts_end + type(timedelta) :: dt + + dt = timedelta(seconds=dels) + ktauday = nint(d2s / dels) + + ! Set the calendar- some forcing cases have special rules + select case (trim(cable_user%mettype)) + case ("plum") + if (plume%leapyears) then + call setcalendar("gregorian") + else + call setcalendar("noleaps") + endif + + case ("cru") + call setcalendar("noleaps") + + case default + if (leaps) then + call setcalendar("gregorian") + else + call setcalendar("noleaps") + endif + + end select + ! END header ! INISTUFF @@ -280,19 +307,25 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) NREP: DO RRRR = 1, NRRRR YEAR: DO YYYY= CABLE_USER%YearStart, CABLE_USER%YearEnd + + ts_start = datetime(year=YYYY, month=1, day=1) + ts_end = ts_start + dt + + LOY = daysInYear(YYYY) - CurYear = YYYY + ! Set kend for all cases bar princeton and gswp + kend = ktauday * LOY !ccc Set "calendar" for netcdf time attribute and ! number of days in the year - calendar = "noleap" - LOY = 365 - IF ( leaps ) THEN - calendar = "standard" - END IF - IF ( leaps .AND. IS_LEAPYEAR( YYYY ) ) THEN - LOY = 366 - END IF + !calendar = "noleap" + !LOY = 365 + !IF ( leaps ) THEN + !calendar = "standard" + !END IF + !IF ( leaps .AND. IS_LEAPYEAR( YYYY ) ) THEN + !LOY = 366 + !END IF ! Check for gswp run SELECT CASE (TRIM(cable_user%MetType)) @@ -307,7 +340,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) IF ( RRRR .EQ. 1 ) THEN CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) - IF (leaps.AND.is_leapyear(YYYY).AND.kend.EQ.2920) THEN + IF (is_leapyear(YYYY).AND.kend.EQ.2920) THEN STOP 'LEAP YEAR INCOMPATIBILITY WITH INPUT MET !' ENDIF IF ( NRRRR .GT. 1 ) THEN @@ -332,23 +365,12 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) kend = ktauday * LOY ENDIF - CASE ('plum') - ! PLUME experiment setup using WATCH - IF ( .NOT. PLUME%LeapYears ) LOY = 365 - kend = NINT(24.0*3600.0/dels) * LOY - - CASE ('cru') - ! TRENDY experiment using CRU-NCEP - LOY = 365 - kend = NINT(24.0*3600.0/dels) * LOY - CASE ('gswp3') ncciy = CurYear CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) CASE ('site') ! site experiment eg AmazonFace (spinup or transient run type) - kend = NINT(24.0*3600.0/dels) * LOY ! get koffset to add to time-step of sitemet SELECT CASE (TRIM(site%RunType)) CASE ('historical') @@ -373,9 +395,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) END SELECT - ! somethings (e.g. CASA-CNP) only need to be done once per day - ktauday=INT(24.0*3600.0/dels) - ! Checks where parameters and initialisations should be loaded from. ! If they can be found in either the met file or restart file, they will ! load from there, with the met file taking precedence. Otherwise, they'll From 0a75dc14a24659a7a5fa4b4a12a8447228885575 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Tue, 16 Dec 2025 15:32:39 +1100 Subject: [PATCH 03/28] Add datetime module --- CMakeLists.txt | 1 + src/util/datetime_module.f90 | 1449 ++++++++++++++++++++++++++++++++++ 2 files changed, 1450 insertions(+) create mode 100644 src/util/datetime_module.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 354051dd0..95a320b0d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -243,6 +243,7 @@ else() src/params/cable_maths_constants_mod.F90 src/util/cable_runtime_opts_mod.F90 src/util/cable_common.F90 + src/util/datetime_module.f90 src/shared/casa_offline_inout.F90 src/shared/casa_ncdf.F90 src/offline/cable_iovars.F90 diff --git a/src/util/datetime_module.f90 b/src/util/datetime_module.f90 new file mode 100644 index 000000000..297e8c269 --- /dev/null +++ b/src/util/datetime_module.f90 @@ -0,0 +1,1449 @@ +module datetime_module + + use, intrinsic :: iso_fortran_env, only: int64, real32, real64, & + stderr => error_unit + use, intrinsic :: iso_c_binding, only: c_char, c_int, c_int64_t, c_null_char, c_associated, C_PTR + + implicit none + + private + + public :: datetime, timedelta, clock + public :: date2num + public :: datetimeRange + public :: daysInMonth + public :: daysInYear + public :: isLeapYear + public :: num2date + public :: machinetimezone + public :: strptime + public :: epochdatetime + public :: localtime + public :: gmtime + public :: tm2date + public :: tm_struct + public :: c_strftime + public :: c_strptime + public :: set_calendar + + real(real64), parameter :: zero = 0._real64, one = 1._real64 + + ! Constant multipliers to transform a number of some time unit to another + real(real64), parameter :: d2h = 24._real64 ! day -> hour + real(real64), parameter :: h2d = one / d2h ! hour -> day + real(real64), parameter :: d2m = d2h * 60._real64 ! day -> minute + real(real64), parameter :: m2d = one / d2m ! minute -> day + real(real64), parameter :: m2h = one / 60 ! minute -> hour + real(real64), parameter :: s2d = m2d / 60 ! second -> day + real(real64), parameter :: d2s = 86400._real64 ! day -> second + real(real64), parameter :: h2s = 3600._real64 ! hour -> second + real(real64), parameter :: h2m = 60._real64 ! hour -> minute + real(real64), parameter :: s2h = one / h2s ! second -> hour + real(real64), parameter :: m2s = 60._real64 ! minute -> second + real(real64), parameter :: s2m = one / m2s ! second -> minute + + integer, parameter :: MAXSTRLEN = 99 ! maximum string length for strftime + + private :: calendarType, gregorian, julian, noLeaps, three60day + enum, bind(c) + enumerator :: calendarType = 0 + enumerator :: gregorian = 1 + enumerator :: julian = 2 + enumerator :: noLeaps = 3 + enumerator :: three60day = 4 + end enum + + integer(kind(calendarType)) :: calendar = gregorian + + type :: datetime + + private + + integer :: year = 1 ! year [1-HUGE(year)] + integer :: month = 1 ! month in year [1-12] + integer :: day = 1 ! day in month [1-31] + integer :: hour = 0 ! hour in day [0-23] + integer :: minute = 0 ! minute in hour [0-59] + integer :: second = 0 ! second in minute [0-59] + integer :: millisecond = 0 ! milliseconds in second [0-999] + real(real64) :: tz = 0 ! timezone offset from UTC [hours] + + contains + + ! getter functions + procedure, pass(self), public :: getYear + procedure, pass(self), public :: getMonth + procedure, pass(self), public :: getDay + procedure, pass(self), public :: getHour + procedure, pass(self), public :: getMinute + procedure, pass(self), public :: getSecond + procedure, pass(self), public :: getMillisecond + procedure, pass(self), public :: getTz + + ! public methods + procedure, pass(self), public :: isocalendar + procedure, pass(self), public :: isoformat + procedure, pass(self), public :: isValid + procedure, nopass, public :: now + procedure, pass(self), public :: secondsSinceEpoch + procedure, pass(self), public :: strftime + procedure, pass(self), public :: tm + procedure, pass(self), public :: tzOffset + procedure, pass(self), public :: isoweekday + procedure, pass(self), public :: isoweekdayLong + procedure, pass(self), public :: isoweekdayShort + procedure, pass(self), public :: utc + procedure, pass(self), public :: weekday + procedure, pass(self), public :: weekdayLong + procedure, pass(self), public :: weekdayShort + procedure, pass(self), public :: yearday + + ! private methods + procedure, pass(self), private :: addMilliseconds + procedure, pass(self), private :: addSeconds + procedure, pass(self), private :: addMinutes + procedure, pass(self), private :: addHours + procedure, pass(self), private :: addDays + + ! operator overloading procedures + procedure, pass(d0), private :: datetime_plus_timedelta + procedure, pass(d0), private :: timedelta_plus_datetime + procedure, pass(d0), private :: datetime_minus_datetime + procedure, pass(d0), private :: datetime_minus_timedelta + procedure, pass(d0), private :: datetime_eq + procedure, pass(d0), private :: datetime_neq + procedure, pass(d0), private :: datetime_gt + procedure, pass(d0), private :: datetime_ge + procedure, pass(d0), private :: datetime_lt + procedure, pass(d0), private :: datetime_le + + generic :: operator(+) => datetime_plus_timedelta, & + timedelta_plus_datetime + generic :: operator(-) => datetime_minus_datetime, & + datetime_minus_timedelta + generic :: operator(==) => datetime_eq + generic :: operator(/=) => datetime_neq + generic :: operator(>) => datetime_gt + generic :: operator(>=) => datetime_ge + generic :: operator(<) => datetime_lt + generic :: operator(<=) => datetime_le + + end type datetime + + interface datetime + module procedure :: datetime_constructor + endinterface datetime + + type :: timedelta + private + + integer :: days = 0 + integer :: hours = 0 + integer :: minutes = 0 + integer :: seconds = 0 + integer :: milliseconds = 0 + + contains + + procedure, pass(self), public :: getDays + procedure, pass(self), public :: getHours + procedure, pass(self), public :: getMinutes + procedure, pass(self), public :: getSeconds + procedure, pass(self), public :: getMilliseconds + + procedure, public :: total_seconds + + procedure, private :: timedelta_plus_timedelta + procedure, private :: timedelta_minus_timedelta + procedure, private :: unary_minus_timedelta + procedure, private :: timedelta_eq + procedure, private :: timedelta_neq + procedure, private :: timedelta_gt + procedure, private :: timedelta_ge + procedure, private :: timedelta_lt + procedure, private :: timedelta_le + + generic :: operator(+) => timedelta_plus_timedelta + generic :: operator(-) => timedelta_minus_timedelta, unary_minus_timedelta + generic :: operator(==) => timedelta_eq + generic :: operator(/=) => timedelta_neq + generic :: operator(>) => timedelta_gt + generic :: operator(>=) => timedelta_ge + generic :: operator(<) => timedelta_lt + generic :: operator(<=) => timedelta_le + + end type timedelta + + interface timedelta + module procedure :: timedelta_constructor + endinterface timedelta + + type,bind(c) :: tm_struct + ! Derived type for compatibility with C and C++ struct tm. + ! Enables calling strftime and strptime using iso_c_binding. + ! See http://www.cplusplus.com/reference/ctime/tm for reference. + integer(c_int) :: tm_sec = 0 ! Seconds [0-60] (1 leap second) + integer(c_int) :: tm_min = 0 ! Minutes [0-59] + integer(c_int) :: tm_hour = 0 ! Hours [0-23] + integer(c_int) :: tm_mday = 0 ! Day [1-31] + integer(c_int) :: tm_mon = 0 ! Month [0-11] + integer(c_int) :: tm_year = 0 ! Year - 1900 + integer(c_int) :: tm_wday = 0 ! Day of week [0-6] + integer(c_int) :: tm_yday = 0 ! Days in year [0-365] + integer(c_int) :: tm_isdst = 0 ! DST [-1/0/1] + end type tm_struct + + interface + + type(c_ptr) function c_strftime(str, slen, format, tm) bind(c, name='strftime') + ! Returns a formatted time string, given input time struct and format. + ! See https://www.cplusplus.com/reference/ctime/strftime for reference. + import :: c_char, c_int, tm_struct, C_PTR + character(kind=c_char), intent(out) :: str(*) ! result string + integer(c_int), value, intent(in) :: slen ! string length + character(kind=c_char), intent(in) :: format(*) ! time format + type(tm_struct), intent(in) :: tm ! tm_struct instance + end function c_strftime + + integer(c_int) function c_strptime(str,format,tm) bind(c,name='strptime') + ! Interface to POSIX strptime. + ! Returns a time struct object based on the input time string str and format. + ! See http://man7.org/linux/man-pages/man3/strptime.3.html for reference. + import :: c_char, c_int, tm_struct + character(kind=c_char), intent(in) :: str(*) ! input string + character(kind=c_char), intent(in) :: format(*) ! time format + type(tm_struct), intent(out) :: tm ! result tm_struct + end function c_strptime + + end interface + + + type :: clock + type(datetime) :: startTime + type(datetime) :: stopTime + type(datetime) :: currentTime + type(timedelta) :: tickInterval + logical :: alarm = .false. + logical :: started = .false. + logical :: stopped = .false. + contains + procedure :: reset + procedure :: tick + end type clock + +contains + + + pure elemental subroutine reset(self) + ! Resets the clock to its start time. + class(clock), intent(in out) :: self + self % currentTime = self % startTime + self % started = .false. + self % stopped = .false. + end subroutine reset + + + pure elemental subroutine tick(self) + ! Increments the currentTime of the clock instance by one tickInterval. + class(clock), intent(in out) :: self + if (self % stopped) return + if (.not. self % started) then + self % started = .true. + self % currentTime = self % startTime + end if + self % currentTime = self % currentTime + self % tickInterval + if (self % currentTime >= self % stopTime) self % stopped = .true. + end subroutine tick + + subroutine set_calendar(calendarString) + ! Set the calendar for the module + character(len=*), intent(in) :: calendarString + + if (trim(calendarString) == "gregorian") then + calendar = gregorian + elseif (trim(calendarString) == "julian") then + calendar = julian + elseif (trim(calendarString) == "noleaps") then + calendar = noLeaps + elseif (trim(calendarString) == "360day") then + calendar = three60day + else + write(stderr,*) calendarString//" is not a valid calendar. "//& + "Valid calendars are gregorian, julian, noleaps or 360day." + end if + + end subroutine set_calendar + + pure elemental type(datetime) function datetime_constructor( & + year, month, day, hour, minute, second, millisecond, tz) + ! Constructor function for the `datetime` class. + integer, intent(in), optional :: year, month, day, hour, minute, second, millisecond + real(real64), intent(in), optional :: tz ! timezone offset in hours + + datetime_constructor % year = 1 + if (present(year)) datetime_constructor % year = year + + datetime_constructor % month = 1 + if (present(month)) datetime_constructor % month = month + + datetime_constructor % day = 1 + if (present(day)) datetime_constructor % day = day + + datetime_constructor % hour = 0 + if (present(hour)) datetime_constructor % hour = hour + + datetime_constructor % minute = 0 + if (present(minute)) datetime_constructor % minute = minute + + datetime_constructor % second = 0 + if (present(second)) datetime_constructor % second = second + + datetime_constructor % millisecond = 0 + if (present(millisecond)) datetime_constructor % millisecond = millisecond + + datetime_constructor % tz = 0 + if (present(tz)) datetime_constructor % tz = tz + + end function datetime_constructor + + + pure elemental integer function getYear(self) + ! Returns the year component + class(datetime), intent(in) :: self + getYear = self % year + end function getYear + + + pure elemental integer function getMonth(self) + ! Returns the year component + class(datetime), intent(in) :: self + getMonth = self % month + end function getMonth + + + pure elemental integer function getDay(self) + ! Returns the year component + class(datetime), intent(in) :: self + getDay = self % day + end function getDay + + + pure elemental integer function getHour(self) + ! Returns the year component + class(datetime), intent(in) :: self + getHour = self % hour + end function getHour + + + pure elemental integer function getMinute(self) + ! Returns the year component + class(datetime), intent(in) :: self + getMinute = self % minute + end function getMinute + + + pure elemental integer function getSecond(self) + ! Returns the year component + class(datetime), intent(in) :: self + getSecond = self % second + end function getSecond + + + pure elemental integer function getMillisecond(self) + ! Returns the year component + class(datetime), intent(in) :: self + getMillisecond = self % millisecond + end function getMillisecond + + + pure elemental real(real64) function getTz(self) + ! Returns the timezone offset component + class(datetime), intent(in) :: self + getTz = self % tz + end function getTz + + + pure elemental subroutine addMilliseconds(self, ms) + ! Adds an integer number of milliseconds to self. Called by `datetime` + ! addition (`+`) and subtraction (`-`) operators. + class(datetime), intent(in out) :: self + integer, intent(in) :: ms + self % millisecond = self % millisecond + ms + do + if (self % millisecond >= 1000) then + call self % addSeconds(self % millisecond / 1000) + self % millisecond = mod(self % millisecond, 1000) + else if (self % millisecond < 0) then + call self % addSeconds(self % millisecond / 1000 - 1) + self % millisecond = mod(self % millisecond, 1000) + 1000 + else + exit + end if + end do + end subroutine addMilliseconds + + + pure elemental subroutine addSeconds(self, s) + ! Adds an integer number of seconds to self. Called by `datetime` + ! addition (`+`) and subtraction (`-`) operators. + class(datetime), intent(in out) :: self + integer, intent(in) :: s + self % second = self % second + s + do + if (self % second >= 60) then + call self % addMinutes(self % second / 60) + self % second = mod(self % second, 60) + else if (self % second < 0) then + call self % addMinutes(self % second / 60 - 1) + self % second = mod(self % second, 60) + 60 + else + exit + end if + end do + end subroutine addSeconds + + + pure elemental subroutine addMinutes(self,m) + ! Adds an integer number of minutes to self. Called by `datetime` + ! addition (`+`) and subtraction (`-`) operators. + class(datetime), intent(in out) :: self + integer, intent(in) :: m + self % minute = self % minute + m + do + if (self % minute >= 60) then + call self % addHours(self % minute / 60) + self % minute = mod(self % minute, 60) + else if (self % minute < 0) then + call self % addHours(self % minute / 60 - 1) + self % minute = mod(self % minute, 60) + 60 + else + exit + end if + end do + end subroutine addMinutes + + + pure elemental subroutine addHours(self,h) + ! Adds an integer number of hours to self. Called by `datetime` + ! addition (`+`) and subtraction (`-`) operators. + class(datetime), intent(in out) :: self + integer, intent(in) :: h + self % hour = self % hour + h + do + if (self % hour >= 24) then + call self % addDays(self % hour / 24) + self % hour = mod(self % hour, 24) + else if (self % hour < 0) then + call self % addDays(self % hour / 24 - 1) + self % hour = mod(self % hour, 24) + 24 + else + exit + end if + end do + end subroutine addHours + + + pure elemental subroutine addDays(self, d) + ! Adds an integer number of dayss to self. Called by `datetime` + ! addition (`+`) and subtraction (`-`) operators. + class(datetime), intent(in out) :: self + integer, intent(in) :: d + integer :: daysInCurrentMonth + self % day = self % day + d + do + daysInCurrentMonth = daysInMonth(self % month, self % year) + if (self % day > daysInCurrentMonth) then + self % day = self % day - daysInCurrentMonth + self % month = self % month+1 + if (self % month > 12) then + self % year = self % year + self % month/12 + self % month = mod(self % month, 12) + end if + else if (self % day < 1) then + self % month = self % month-1 + if (self % month < 1) then + self % year = self % year + self % month / 12 - 1 + self % month = 12 + mod(self % month, 12) + end if + self % day = self % day + daysInMonth(self % month, self % year) + else + exit + end if + end do + end subroutine addDays + + + pure elemental character(23) function isoformat(self,sep) + ! Returns character string with time in ISO 8601 format. + class(datetime), intent(in) :: self + character, intent(in), optional :: sep + character :: separator + + separator = 'T' + if (present(sep)) separator = sep + + ! TODO below is a bit cumbersome and was implemented + ! at a time before the interface to strftime. Now we + ! could do something like: + ! + ! isoformat = self % strftime('%Y-%m-%d'//separator//'%H:%M:%S') + ! + isoformat = int2str(self % year, 4)//'-'// & + int2str(self % month, 2)//'-'// & + int2str(self % day, 2)//separator//& + int2str(self % hour, 2)//':'// & + int2str(self % minute, 2)//':'// & + int2str(self % second, 2)//'.'// & + int2str(self % millisecond,3) + + end function isoformat + + + pure elemental logical function isValid(self) + ! Checks whether the `datetime` instance has valid component values. + ! Returns `.true.` if the `datetime` instance is valid, and `.false.` + ! otherwise. + class(datetime), intent(in) :: self + + ! assume valid + isValid = .true. + + if (self % year < 1) then + isValid = .false. + return + end if + + if (self % month < 1 .or. self % month > 12) then + isValid = .false. + return + end if + + if (self % day < 1 .or. & + self % day > daysInMonth(self % month,self % year)) then + isValid = .false. + return + end if + + if (self % hour < 0 .or. self % hour > 23) then + isValid = .false. + return + end if + + if (self % minute < 0 .or. self % minute > 59) then + isValid = .false. + return + end if + + if (self % second < 0 .or. self % second > 59) then + isValid = .false. + return + end if + + if (self % millisecond < 0 .or. self % millisecond > 999) then + isValid = .false. + return + end if + + end function isValid + + + type(datetime) function now() + ! Returns a `datetime` instance with current time. + character(5) :: zone + integer :: values(8) + integer :: hour, minute + + ! Obtain local machine time zone information + call date_and_time(zone=zone, values=values) + + read(zone(1:3), '(i3)') hour + read(zone(4:5), '(i2)') minute + + now = datetime(year = values(1), month = values(2), day = values(3), & + hour = values(5), minute = values(6), second = values(7), & + millisecond = values(8)) + + now % tz = hour + minute * m2h + + end function now + + + pure elemental integer function weekday(self) + ! Returns the day of the week calculated using Zeller's congruence. + ! Returned value is an integer scalar in the range [0-6], such that: + ! + ! 0: Sunday + ! 1: Monday + ! 2: Tuesday + ! 3: Wednesday + ! 4: Thursday + ! 5: Friday + ! 6: Saturday + class(datetime), intent(in) :: self + integer :: year, month, j, k + + year = self % year + month = self % month + + if (month <= 2) then + month = month + 12 + year = year - 1 + end if + + j = year / 100 + k = mod(year, 100) + + ! Assume other calendars return nonsense + if (calendar == gregorian) then + weekday = mod(self % day + ((month + 1) * 26) / 10 + k + k / 4 + j / 4 + 5 * j, 7) -1 + elseif (calendar == julian) then + weekday = mod(self % day + ((month + 1) * 26) / 10 + k + k / 4 + 5 - j, 7) - 1 + end if + + if (weekday < 0) weekday = 6 + + end function weekday + + + pure elemental integer function isoweekday(self) + ! Returns the day of the week per ISO 8601 returned from weekday(). + ! Returned value is an integer scalar in the range [1-7]. + class(datetime), intent(in) :: self + isoweekday = self % weekday() + if (isoweekday == 0) isoweekday = 7 + end function isoweekday + + + pure elemental character(9) function weekdayLong(self) + ! Returns the full name of the day of the week. + class(datetime), intent(in) :: self + character(9), parameter :: & + days(*) = ['Sunday ', 'Monday ', 'Tuesday ','Wednesday', & + 'Thursday ', 'Friday ', 'Saturday '] + weekdayLong = days(self % weekday() + 1) + end function weekdayLong + + + pure elemental character(9) function isoweekdayLong(self) + ! Returns the full name of the day of the week for ISO 8601 + ! ordered weekdays. + class(datetime), intent(in) :: self + character(9), parameter :: & + days(7) = ['Monday ','Tuesday ','Wednesday','Thursday ', & + 'Friday ','Saturday ','Sunday '] + isoweekdayLong = days(self % isoweekday()) + end function isoweekdayLong + + + pure elemental character(3) function weekdayShort(self) + ! Returns the short (3-letter) name of the day of the week. + class(datetime), intent(in) :: self + character(3), parameter :: days(7) = ['Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat'] + weekdayShort = days(self % weekday() + 1) + end function weekdayShort + + + pure elemental character(3) function isoweekdayShort(self) + ! Returns the short (3-letter) name of the day of the week + ! based on ISO 8601 ordering. + class(datetime), intent(in) :: self + character(3), parameter :: days(7) = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun'] + isoweekdayShort = days(self % isoweekday()) + end function isoweekdayShort + + + function isocalendar(self) + ! Returns an array of 3 integers, year, week number, and week day, + ! as defined by ISO 8601 week date. Essentially a wrapper around C + ! `strftime` function. + class(datetime), intent(in) :: self + integer :: isocalendar(3) + integer :: year, week, wday + type(C_PTR) :: rc + character(20) :: string + + rc = c_strftime(string, len(string), '%G %V %u' // c_null_char, self % tm()) + if (.not. c_associated(rc)) then + write(stderr,*) "ERROR:datetime:strftime: format: %G %V %u" + return + endif + + read(string(1:4), '(i4)') year + read(string(6:7), '(i2)') week + read(string(9:9), '(i1)') wday + + isocalendar = [year, week, wday] + + end function isocalendar + + + integer(int64) function secondsSinceEpoch(self) + ! Returns an integer number of seconds since the UNIX Epoch (1 Jan 1970). + ! Since Windows does not have strftime('%s'), we implement this using + ! datetime itself. + class(datetime), intent(in) :: self + type(timedelta) :: delta + type(datetime) :: this_time, unix_time + integer :: sign, hours, minutes, tzsec + + this_time = datetime(self%year, self%month, self%day, & + self%hour, self%minute, self%second) + unix_time = datetime(1970, 1, 1, 0, 0, 0) + delta = this_time - unix_time + secondsSinceEpoch = delta%total_seconds() + + if(self % tz == 0_real64) return + + ! affect timezone information + if(self % tz < 0_real64) then + sign = -1 + else + sign = 1 + end if + + hours = int(abs(self % tz)) + minutes = nint((abs(self % tz) - hours) * 60) + + if (minutes == 60) then + minutes = 0 + hours = hours + 1 + end if + + tzsec = sign * (hours * h2s + minutes) + secondsSinceEpoch = secondsSinceEpoch - tzsec + + end function secondsSinceEpoch + + + function strftime(self, format) + ! Wrapper around C and C++ `strftime` function. + class(datetime), intent(in) :: self + character(*), intent(in) :: format + character(:), allocatable :: strftime + type(C_PTR) :: rc + character(MAXSTRLEN) :: resultString + resultString = "" + rc = c_strftime(resultString, len(resultString), trim(format) // c_null_char, & + self % tm()) + if (.not. c_associated(rc)) write(stderr, '(a)') "ERROR:datetime:strftime: format: " // trim(format) + strftime = resultString(1:len_trim(resultString)-1) !< strip null + end function strftime + + + pure elemental type(tm_struct) function tm(self) + ! Returns a `tm_struct` instance of the current `datetime`. + class(datetime), intent(in) :: self + tm % tm_sec = self % second + tm % tm_min = self % minute + tm % tm_hour = self % hour + tm % tm_mday = self % day + tm % tm_mon = self % month - 1 + tm % tm_year = self % year - 1900 + tm % tm_wday = self % weekday() + tm % tm_yday = self % yearday() - 1 + tm % tm_isdst = -1 + end function tm + + + pure elemental character(5) function tzOffset(self) + ! Returns a character string with timezone offset in hours from UTC, + ! in format +/-[hh][mm]. + class(datetime), intent(in) :: self + integer :: hours,minutes + + if (self % tz < 0) then + tzOffset(1:1) = '-' + else + tzOffset(1:1) = '+' + end if + + hours = int(abs(self % tz)) + minutes = nint((abs(self % tz) - hours) * 60) + + if (minutes == 60) then + minutes = 0 + hours = hours + 1 + end if + + write(tzOffset(2:5), '(2i2.2)') hours, minutes + + end function tzOffset + + + pure elemental type(datetime) function utc(self) + ! Returns the `datetime` instance at Coordinated Universal Time (UTC). + class(datetime), intent(in) :: self + integer :: hours, minutes, sgn + hours = int(abs(self % tz)) + minutes = nint((abs(self % tz) - hours) * 60) + sgn = int(sign(one, self % tz)) + utc = self - timedelta(hours=sgn * hours, minutes=sgn * minutes) + utc % tz = 0 + end function utc + + + pure elemental integer function yearday(self) + ! Returns the integer day of the year (ordinal date). + class(datetime), intent(in) :: self + integer :: month + yearday = 0 + do month = 1, self % month-1 + yearday = yearday + daysInMonth(month, self % year) + end do + yearday = yearday + self % day + end function yearday + + + pure elemental function datetime_plus_timedelta(d0,t) result(d) + ! Adds a `timedelta` instance to a `datetime` instance, and returns a + ! new `datetime` instance. Overloads the operator `+`. + class(datetime), intent(in) :: d0 + class(timedelta), intent(in) :: t + type(datetime) :: d + + integer :: milliseconds, seconds, minutes, hours, days + + d = datetime(year = d0 % getYear(), & + month = d0 % getMonth(), & + day = d0 % getDay(), & + hour = d0 % getHour(), & + minute = d0 % getMinute(), & + second = d0 % getSecond(), & + millisecond = d0 % getMillisecond(), & + tz = d0 % getTz()) + + milliseconds = t % getMilliseconds() + seconds = t % getSeconds() + minutes = t % getMinutes() + hours = t % getHours() + days = t % getDays() + + if (milliseconds /= 0) call d % addMilliseconds(milliseconds) + if (seconds /= 0) call d % addSeconds(seconds) + if (minutes /= 0) call d % addMinutes(minutes) + if (hours /= 0) call d % addHours(hours) + if (days /= 0) call d % addDays(days) + + end function datetime_plus_timedelta + + + pure elemental function timedelta_plus_datetime(t,d0) result(d) + ! Adds a `timedelta` instance to a `datetime` instance, and returns a + ! new `datetime` instance. Overloads the operator `+`. + class(timedelta), intent(in) :: t + class(datetime), intent(in) :: d0 + type(datetime) :: d + d = d0 + t + end function timedelta_plus_datetime + + + pure elemental function datetime_minus_timedelta(d0,t) result(d) + ! Subtracts a `timedelta` instance from a `datetime` instance and + ! returns a new `datetime` instance. Overloads the operator `-`. + class(datetime), intent(in) :: d0 + class(timedelta), intent(in) :: t + type(datetime) :: d + d = d0 + (-t) + end function datetime_minus_timedelta + + + pure elemental function datetime_minus_datetime(d0,d1) result(t) + ! Subtracts a `datetime` instance from another `datetime` instance, + ! and returns a `timedelta` instance. Overloads the operator `-`. + class(datetime), intent(in) :: d0, d1 + type(timedelta) :: t + real(real64) :: daysDiff + integer :: days,hours,minutes,seconds,milliseconds + integer :: sign_ + + daysDiff = date2num(d0)-date2num(d1) + + if (daysDiff < 0) then + sign_ = -1 + daysDiff = ABS(daysDiff) + else + sign_ = 1 + end if + + days = int(daysDiff) + hours = int((daysDiff-days)*d2h) + minutes = int((daysDiff-days-hours*h2d)*d2m) + seconds = int((daysDiff-days-hours*h2d-minutes*m2d)*d2s) + milliseconds = nint((daysDiff-days-hours*h2d-minutes*m2d& + -seconds*s2d)*d2s*1e3_real64) + + t = timedelta(sign_*days,sign_*hours,sign_*minutes,sign_*seconds,& + sign_*milliseconds) + + end function datetime_minus_datetime + + + pure elemental logical function datetime_gt(d0,d1) result(res) + ! `datetime` comparison operator that returns `.true.` if `d0` is + ! greater than `d1` and `.false.` otherwise. Overloads the + ! operator `>`. + class(datetime), intent(in) :: d0, d1 + type(datetime) :: d0_utc, d1_utc + + ! Convert to UTC before making comparison + d0_utc = d0 % utc() + d1_utc = d1 % utc() + + ! Compare years + if (d0_utc % year > d1_utc % year) then + res = .true. + else if (d0_utc % year < d1_utc % year) then + res = .false. + else + + ! Compare months + if (d0_utc % month > d1_utc % month) then + res = .true. + else if (d0_utc % month < d1_utc % month) then + res = .false. + else + + ! Compare days + if (d0_utc % day > d1_utc % day) then + res = .true. + else if (d0_utc % day < d1_utc % day) then + res = .false. + else + + ! Compare hours + if (d0_utc % hour > d1_utc % hour) then + res = .true. + else if (d0_utc % hour < d1_utc % hour) then + res = .false. + else + + ! Compare minutes + if (d0_utc % minute > d1_utc % minute) then + res = .true. + else if (d0_utc % minute < d1_utc % minute) then + res = .false. + else + + ! Compare seconds + if (d0_utc % second > d1_utc % second) then + res = .true. + else if (d0_utc % second < d1_utc % second) then + res = .false. + else + + ! Compare milliseconds + if (d0_utc % millisecond > d1_utc % millisecond) then + res = .true. + else + res = .false. + end if + + end if + end if + end if + end if + end if + end if + + end function datetime_gt + + + pure elemental logical function datetime_lt(d0,d1) result(res) + ! `datetime` comparison operator that returns `.true.` if `d0` is + ! less than `d1` and `.false.` otherwise. Overloads the operator `<`. + class(datetime), intent(in) :: d0, d1 + res = d1 > d0 + end function datetime_lt + + + pure elemental logical function datetime_eq(d0,d1) result(res) + ! `datetime` comparison operator that returns `.true.` if `d0` is + ! equal to `d1` and `.false.` otherwise. Overloads the operator `==`. + class(datetime), intent(in) :: d0, d1 + type(datetime) :: d0_utc, d1_utc + + ! Convert to UTC before making comparison + d0_utc = d0 % utc() + d1_utc = d1 % utc() + + res = d0_utc % year == d1_utc % year .and. & + d0_utc % month == d1_utc % month .and. & + d0_utc % day == d1_utc % day .and. & + d0_utc % hour == d1_utc % hour .and. & + d0_utc % minute == d1_utc % minute .and. & + d0_utc % second == d1_utc % second .and. & + d0_utc % millisecond == d1_utc % millisecond + + end function datetime_eq + + + pure elemental logical function datetime_neq(d0,d1) result(res) + ! `datetime` comparison operator that eturns `.true.` if `d0` is + ! not equal to `d1` and `.false.` otherwise. Overloads the operator `/=`. + class(datetime), intent(in) :: d0, d1 + res = .not. d0 == d1 + end function datetime_neq + + + pure elemental logical function datetime_ge(d0,d1) result(res) + ! `datetime` comparison operator. Returns `.true.` if `d0` is greater + ! than or equal to `d1` and `.false.` otherwise. Overloads the + ! operator `>=`. + class(datetime), intent(in) :: d0, d1 + res = d0 > d1 .or. d0 == d1 + end function datetime_ge + + + pure elemental logical function datetime_le(d0,d1) result(res) + ! `datetime` comparison operator. Returns `.true.` if `d0` is less + ! than or equal to `d1`, and `.false.` otherwise. Overloads the + ! operator `<=`. + class(datetime), intent(in) :: d0, d1 + res = d1 > d0 .or. d0 == d1 + end function datetime_le + + + pure elemental logical function isLeapYear(year) + ! Returns `.true.` if year is leap year and `.false.` otherwise. + integer, intent(in) :: year + + if (calendar == gregorian) then + isLeapYear = (mod(year,4) == 0 .and. .not. mod(year,100) == 0)& + .or. (mod(year,400) == 0) + elseif (calendar == julian) then + isLeapYear = mod(year,4) == 0 + elseif (calendar == noLeaps .or. calendar == three60day) then + isLeapYear = .false. + end if + + end function isLeapYear + + + pure function datetimeRange(d0, d1, t) + ! Given start and end `datetime` instances `d0` and `d1` and time + ! increment as `timedelta` instance `t`, returns an array of + ! `datetime` instances. The number of elements is the number of whole + ! time increments contained between datetimes `d0` and `d1`. + type(datetime), intent(in) :: d0, d1 + type(timedelta), intent(in) :: t + real(real64) :: datenum0, datenum1, eps, increment + type(datetime), allocatable :: datetimeRange(:) + integer :: n, nm + eps = 1e-10_real64 + datenum0 = date2num(d0) + datenum1 = date2num(d1) + increment = t % total_seconds() * s2d + nm = floor((datenum1 - datenum0 + eps) / increment) + 1 + allocate(datetimeRange(nm)) + do n = 1, nm + datetimeRange(n) = num2date(datenum0 + (n - 1) * increment) + end do + end function datetimeRange + + + pure elemental integer function daysInMonth(month,year) + ! Given integer month and year, returns an integer number + ! of days in that particular month. + integer, intent(in) :: month, year + + integer :: days(12) + + if (calendar == three60day) then + days = [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30] + else + days = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + end if + + if (month < 1 .or. month > 12) then + ! Should raise an error and abort here, however we want to keep + ! the pure and elemental attributes. Make sure this function is + ! called with the month argument in range. + daysInMonth = 0 + return + end if + + if (month == 2 .and. isLeapYear(year)) then + daysInMonth = 29 + else + daysInMonth = days(month) + end if + + end function daysInMonth + + + pure elemental integer function daysInYear(year) + ! Returns the number of days in year. + integer, intent(in) :: year + if (calendar == three60day) then + daysInYear = 360 + else + if (isLeapYear(year)) then + daysInYear = 366 + else + daysInYear = 365 + end if + end if + end function daysInYear + + + pure elemental real(real64) function date2num(d) + ! Given a datetime instance d, returns number of days since + ! `0001-01-01 00:00:00`, taking into account the timezone offset. + type(datetime), intent(in) :: d + type(datetime) :: d_utc + integer :: year + + ! Convert to UTC first + d_utc = d % utc() + + ! d_utc % year must be positive: + if (d_utc % year < 1) then + date2num = 0 + return + end if + + date2num = 0 + do year = 1,d_utc % year-1 + date2num = date2num + daysInYear(year) + end do + + date2num = date2num & + + d_utc % yearday() & + + d_utc % hour*h2d & + + d_utc % minute*m2d& + + (d_utc % second+1e-3_real64*d_utc % millisecond)*s2d + + end function date2num + + + pure elemental type(datetime) function num2date(num) + ! Given number of days since `0001-01-01 00:00:00`, returns a + ! correspoding `datetime` instance. + real(real64), intent(in) :: num + integer :: year, month, day, hour, minute, second, millisecond + real(real64) :: days, totseconds + + ! num must be positive + if (num < 0) then + num2date = datetime(1) + return + end if + + days = num + + year = 1 + do + if (int(days) <= daysInYear(year))exit + days = days-daysInYear(year) + year = year+1 + end do + + month = 1 + do + if (inT(days) <= daysInMonth(month,year))exit + days = days-daysInMonth(month,year) + month = month+1 + end do + + day = int(days) + totseconds = (days-day)*d2s + hour = int(totseconds*s2h) + minute = int((totseconds-hour*h2s)*s2m) + second = int(totseconds-hour*h2s-minute*m2s) + millisecond = nint((totseconds-int(totseconds))*1e3_real64) + + num2date = datetime(year,month,day,hour,minute,second,millisecond,tz=zero) + + ! Handle a special case caused by floating-point arithmethic: + if (num2date % millisecond == 1000) then + num2date % millisecond = 0 + call num2date % addSeconds(1) + end if + + if (num2date % second == 60) then + num2date % second = 0 + call num2date % addMinutes(1) + end if + if (num2date % minute == 60) then + num2date % minute = 0 + call num2date % addHours(1) + end if + if (num2date % hour == 24) then + num2date % hour = 0 + call num2date % addDays(1) + end if + + end function num2date + + + real(real64) function machinetimezone() + ! Return a real value instance of local machine's timezone. + character(len=5) :: zone + integer :: values(8) + integer :: hour, minute + + ! Obtain local machine time zone information + call date_and_time(zone=zone, values=values) + read(zone(1:3), '(i3)') hour + read(zone(4:5), '(i2)') minute + + if(hour<0)then + machinetimezone = real(hour, kind=real64) - real(minute, kind=real64) * m2h + else + machinetimezone = real(hour, kind=real64) + real(minute, kind=real64) * m2h + end if + end function machinetimezone + + + type(datetime) function strptime(str,format,tz) + ! A wrapper function around C/C++ strptime function. + ! Returns a `datetime` instance. + character(*), intent(in) :: str, format + real(real64), intent(in), optional :: tz + integer :: rc + type(tm_struct) :: tm + rc = c_strptime(trim(str) // c_null_char, trim(format) // c_null_char, tm) + if (rc == 0) then + write(stderr, *) "ERROR:datetime:strptime: failed to parse string: ", str + return + endif + strptime = tm2date(tm,tz) + end function strptime + + + pure elemental type(datetime) function epochdatetime() + epochdatetime = datetime(1970,1,1,0,0,0,0,tz=zero) + end function epochdatetime + + + pure elemental type(datetime) function localtime(epoch, tz) + ! Returns a `datetime` instance from epoch. + ! tz can be obtained from `machinetimezone` + integer(int64),intent(in) :: epoch + real(real64),intent(in) :: tz !! local machine time zone information + type(datetime) :: datetime_from_epoch + type(timedelta) :: td + integer :: day, sec + integer(int64) :: localseconds + + datetime_from_epoch = epochdatetime() + + localseconds = nint(tz * h2s) + epoch + !suppress overflow + day = floor(localseconds/d2s, kind=real32) + sec = localseconds - day * d2s + td = timedelta(days=day, seconds=sec) + datetime_from_epoch % tz = tz + localtime = datetime_from_epoch + td + end function localtime + + + pure elemental type(datetime) function gmtime(epoch) + ! Returns a `datetime` instance from epoch. + integer(int64),intent(in) :: epoch + type(datetime) :: datetime_from_epoch + type(timedelta) :: td + integer :: day, sec + datetime_from_epoch = epochdatetime() + !suppress overflow + day = floor(epoch/d2s, kind=real32) + sec = epoch - day * d2s + td = timedelta(days=day, seconds=sec) + gmtime = datetime_from_epoch + td + end function gmtime + + + pure elemental type(datetime) function tm2date(ctime, tz) + ! Given a `tm_struct` instance, returns a corresponding `datetime` + ! instance. + type(tm_struct), intent(in) :: ctime + real(real64), intent(in), optional :: tz ! time zone + + tm2date % millisecond = 0 + tm2date % second = ctime % tm_sec + tm2date % minute = ctime % tm_min + tm2date % hour = ctime % tm_hour + tm2date % day = ctime % tm_mday + tm2date % month = ctime % tm_mon+1 + tm2date % year = ctime % tm_year+1900 + + ! tm_struct have no information of timze zone. + ! but if you run this library with C language's time.h, + ! localtime function deals system's timezone. + ! So, if you want to similar way, you can set tz value with + ! this library's `machinetimezone` function. + if(present(tz))then + tm2date % tz = tz + else + tm2date % tz = 0.0_real64 + end if + + end function tm2date + + + pure function int2str(i, length) + ! Converts an integer `i` into a character string of requested length, + ! pre-pending zeros if necessary. + integer, intent(in) :: i, length + character(length) :: int2str + character(2) :: string + write(string, '(i2)') length + write(int2str, '(i' // string // '.' // string //')') i + end function int2str + + + pure elemental type(timedelta) function timedelta_constructor( & + days, hours, minutes, seconds, milliseconds) + ! Constructor function for the `timedelta` class. + integer, intent(in), optional :: days, hours, minutes, seconds, milliseconds + + timedelta_constructor % days = 0 + if (present(days)) timedelta_constructor % days = days + + timedelta_constructor % hours = 0 + if (present(hours)) timedelta_constructor % hours = hours + + timedelta_constructor % minutes = 0 + if (present(minutes)) timedelta_constructor % minutes = minutes + + timedelta_constructor % seconds = 0 + if (present(seconds)) timedelta_constructor % seconds = seconds + + timedelta_constructor % milliseconds = 0 + if (present(milliseconds)) timedelta_constructor % milliseconds = milliseconds + + end function timedelta_constructor + + + pure elemental integer function getDays(self) + ! Returns the number of days. + class(timedelta), intent(in) :: self + getDays = self % days + end function getDays + + + pure elemental integer function getHours(self) + ! Returns the number of hours. + class(timedelta), intent(in) :: self + getHours = self % hours + end function getHours + + + pure elemental integer function getMinutes(self) + ! Returns the number of minutes. + class(timedelta), intent(in) :: self + getMinutes = self % minutes + end function getMinutes + + + pure elemental integer function getSeconds(self) + ! Returns the number of seconds. + class(timedelta), intent(in) :: self + getSeconds = self % seconds + end function getSeconds + + + pure elemental integer function getMilliseconds(self) + ! Returns the number of milliseconds. + class(timedelta), intent(in) :: self + getMilliseconds = self % milliseconds + end function getMilliseconds + + + pure elemental real(real64) function total_seconds(self) + ! Returns a total number of seconds contained in a `timedelta` + ! instance. + class(timedelta), intent(in) :: self + total_seconds = self % days*86400._real64 & + + self % hours*3600._real64 & + + self % minutes*60._real64 & + + self % seconds & + + self % milliseconds*1e-3_real64 + end function total_seconds + + + pure elemental type(timedelta) function timedelta_plus_timedelta(t0,t1) result(t) + ! Adds two `timedelta` instances together and returns a `timedelta` + ! instance. Overloads the operator `+`. + class(timedelta), intent(in) :: t0, t1 + t = timedelta(days = t0 % days + t1 % days, & + hours = t0 % hours + t1 % hours, & + minutes = t0 % minutes + t1 % minutes, & + seconds = t0 % seconds + t1 % seconds, & + milliseconds = t0 % milliseconds + t1 % milliseconds) + end function timedelta_plus_timedelta + + + pure elemental type(timedelta) function timedelta_minus_timedelta(t0,t1) result(t) + ! Subtracts a `timedelta` instance from another. Returns a + ! `timedelta` instance. Overloads the operator `-`. + class(timedelta), intent(in) :: t0, t1 + t = t0 + (-t1) + end function timedelta_minus_timedelta + + + pure elemental type(timedelta) function unary_minus_timedelta(t0) result(t) + ! Takes a negative of a `timedelta` instance. Overloads the operator `-`. + class(timedelta), intent(in) :: t0 + t % days = -t0 % days + t % hours = -t0 % hours + t % minutes = -t0 % minutes + t % seconds = -t0 % seconds + t % milliseconds = -t0 % milliseconds + end function unary_minus_timedelta + + + pure elemental logical function timedelta_eq(td0,td1) result(res) + ! `timedelta` object comparison operator. Returns `.true.` if `td0` + ! is equal to `td1` and `.false.` otherwise. Overloads the operator + ! `==`. + class(timedelta), intent(in) :: td0, td1 + res = td0 % total_seconds() == td1 % total_seconds() + end function timedelta_eq + + + pure elemental logical function timedelta_neq(td0,td1) result(res) + ! `timedelta` object comparison operator. Returns `.true.` if `td0` + ! is not equal to `td1` and `.false.` otherwise. Overloads the + ! operator `/=`. + class(timedelta), intent(in) :: td0, td1 + res = td0 % total_seconds() /= td1 % total_seconds() + end function timedelta_neq + + + pure elemental logical function timedelta_gt(td0,td1) result(res) + ! `timedelta` object comparison operator. Returns `.true.` if + ! `td0` is greater than `td1` and `.false.` otherwise. Overloads the + ! operator `>`. + class(timedelta), intent(in) :: td0, td1 + res = td0 % total_seconds() > td1 % total_seconds() + end function timedelta_gt + + + pure elemental logical function timedelta_ge(td0,td1) result(res) + ! `timedelta` object comparison operator. Returns `.true.` if `td0` + ! is greater than or equal to `td1` and `.false.` otherwise. + ! Overloads the operator >=. + class(timedelta), intent(in) :: td0, td1 + res = td0 % total_seconds() >= td1 % total_seconds() + end function timedelta_ge + + + pure elemental logical function timedelta_lt(td0,td1) result(res) + ! `timedelta` object comparison operator. Returns `.true.` if `td0` + ! is less than `td1` and `.false.` otherwise. Overloads the operator + ! `<`. + class(timedelta), intent(in) :: td0, td1 + res = td0 % total_seconds() < td1 % total_seconds() + end function timedelta_lt + + + pure elemental logical function timedelta_le(td0,td1) result(res) + ! `timedelta` object comparison operator. Returns `.true.` if `td0` + ! is less than or equal to `td1` and `.false.` otherwise. Overloads + ! the operator `<=`. + class(timedelta), intent(in) :: td0, td1 + res = td0 % total_seconds() <= td1 % total_seconds() + end function timedelta_le + +end module datetime_module From fe75799de59548dbe2f52db01ba28b70f7ef3288 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Tue, 16 Dec 2025 16:02:11 +1100 Subject: [PATCH 04/28] Minor changes to make it compile --- src/offline/cable_serial.F90 | 4 +++- src/util/datetime_module.f90 | 10 +++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index ab421bb6d..3bf952867 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -150,6 +150,8 @@ MODULE cable_serial USE landuse_variable USE bgcdriver_mod, ONLY : bgcdriver USE casa_offline_inout_module, ONLY : WRITE_CASA_RESTART_NC, WRITE_CASA_OUTPUT_NC + +use datetime_module IMPLICIT NONE PRIVATE @@ -271,7 +273,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) type(datetime) :: ts_start, ts_end type(timedelta) :: dt - dt = timedelta(seconds=dels) + dt = timedelta(seconds=int(dels)) ktauday = nint(d2s / dels) ! Set the calendar- some forcing cases have special rules diff --git a/src/util/datetime_module.f90 b/src/util/datetime_module.f90 index 297e8c269..8c46ca87b 100644 --- a/src/util/datetime_module.f90 +++ b/src/util/datetime_module.f90 @@ -6,7 +6,7 @@ module datetime_module implicit none - private + !private public :: datetime, timedelta, clock public :: date2num @@ -24,7 +24,7 @@ module datetime_module public :: tm_struct public :: c_strftime public :: c_strptime - public :: set_calendar + public :: setcalendar real(real64), parameter :: zero = 0._real64, one = 1._real64 @@ -53,7 +53,7 @@ module datetime_module enumerator :: three60day = 4 end enum - integer(kind(calendarType)) :: calendar = gregorian + integer(kind(calendarType)), private :: calendar = gregorian type :: datetime @@ -255,7 +255,7 @@ pure elemental subroutine tick(self) if (self % currentTime >= self % stopTime) self % stopped = .true. end subroutine tick - subroutine set_calendar(calendarString) + subroutine setcalendar(calendarString) ! Set the calendar for the module character(len=*), intent(in) :: calendarString @@ -272,7 +272,7 @@ subroutine set_calendar(calendarString) "Valid calendars are gregorian, julian, noleaps or 360day." end if - end subroutine set_calendar + end subroutine setcalendar pure elemental type(datetime) function datetime_constructor( & year, month, day, hour, minute, second, millisecond, tz) From 4500c9efd0a16ef814f2257109d241234d4c6a87 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Tue, 16 Dec 2025 16:49:06 +1100 Subject: [PATCH 05/28] Add end of period checks and apply them to is_casa_time --- src/offline/cable_serial.F90 | 3 ++- src/shared/casa_ncdf.F90 | 44 ++++++--------------------------- src/util/datetime_module.f90 | 47 ++++++++++++++++++++++++++++++++++++ 3 files changed, 57 insertions(+), 37 deletions(-) diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 3bf952867..a89b8382c 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -311,7 +311,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) YEAR: DO YYYY= CABLE_USER%YearStart, CABLE_USER%YearEnd ts_start = datetime(year=YYYY, month=1, day=1) - ts_end = ts_start + dt LOY = daysInYear(YYYY) @@ -499,6 +498,8 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ! time step loop over ktau DO ktau=kstart, kend + ts_end = ts_start + dt + WRITE(logn,*) 'Current time: ', ts_start WRITE(logn,*) 'Progress -',REAL(ktau)/REAL(kend)*100.0 ! increment total timstep counter diff --git a/src/shared/casa_ncdf.F90 b/src/shared/casa_ncdf.F90 index 593f922b9..0bc47b857 100644 --- a/src/shared/casa_ncdf.F90 +++ b/src/shared/casa_ncdf.F90 @@ -36,6 +36,8 @@ !#define UM_BUILD YES MODULE casa_ncdf_module + use datetime_module, only: datetime, isnewday, isnewmonth, isnewyear + IMPLICIT NONE #ifndef UM_BUILD @@ -402,7 +404,7 @@ SUBROUTINE DOYSOD2YMDHMS( YYYY,DOY,SOD,MM,DD,HOUR,MINUTE,SECOND ) END SUBROUTINE DOYSOD2YMDHMS - FUNCTION IS_CASA_TIME(iotype, yyyy, ktau, kstart, koffset, kend, ktauday, logn) + FUNCTION IS_CASA_TIME(currtime, iotype, logn) USE cable_common_module, ONLY: CABLE_USER ! Correctly determine if it is time to dump-read or standard-write @@ -413,33 +415,12 @@ FUNCTION IS_CASA_TIME(iotype, yyyy, ktau, kstart, koffset, kend, ktauday, logn) !applications. iovars is an offline module and so not appropriate to include !here. Suggested FIX is to move decs of vars needed (e.g. leaps) to here, and !then use common in iovars -#ifdef Vanessas_common - USE cable_IO_vars_module, ONLY: leaps -#endif IMPLICIT NONE LOGICAL :: IS_CASA_TIME - INTEGER ,INTENT(IN) :: yyyy, ktau, kstart, koffset, kend, ktauday, logn + type(datetime), intent(in) :: currtime + INTEGER ,INTENT(IN) :: logn CHARACTER,INTENT(IN) :: iotype*5 - LOGICAL :: is_eod, is_eom, is_eoy - INTEGER :: doy, m - INTEGER, DIMENSION(12) :: MONTH - - is_eom = .FALSE. - is_eoy = .FALSE. - IS_CASA_TIME = .FALSE. - - MONTH = (/ 31,28,31,30,31,30,31,31,30,31,30,31 /) - is_eod = ( MOD((ktau-kstart+1+koffset),ktauday).EQ.0 ) - IF ( .NOT. is_eod ) RETURN ! NO if it is not end of day - -#ifdef Vanessas_common - IF ( IS_LEAPYEAR( YYYY ) .AND. leaps ) THEN - MONTH(2) = 29 - ELSE - MONTH(2) = 28 - ENDIF -#endif ! Check for reading from dump-file (hard-wired to daily casa-timestep) IF ( iotype .EQ. "dread" ) THEN @@ -450,19 +431,10 @@ FUNCTION IS_CASA_TIME(iotype, yyyy, ktau, kstart, koffset, kend, ktauday, logn) ! Check for writing of casa standard output ELSE IF ( iotype .EQ. "write" ) THEN - doy = NINT(REAL(ktau-kstart+1+koffset)/REAL(ktauday)) - DO m = 1, 12 - IF ( doy .EQ. SUM(MONTH(1:m)) ) THEN - is_eom = .TRUE. - IF ( m .EQ. 12 ) is_eoy = .TRUE. - EXIT - ENDIF - END DO - SELECT CASE ( TRIM(CABLE_USER%CASA_OUT_FREQ) ) - CASE ("daily" ) ; IS_CASA_TIME = .TRUE. - CASE ("monthly" ) ; IF ( is_eom ) IS_CASA_TIME = .TRUE. - CASE ("annually") ; IF ( is_eoy ) IS_CASA_TIME = .TRUE. + CASE ("daily" ) ; IS_CASA_TIME = isnewday(currtime) + CASE ("monthly" ) ; IS_CASA_TIME = isnewmonth(currtime) + CASE ("annually") ; IS_CASA_TIME = isnewyear(currtime) END SELECT ELSE WRITE(logn,*)"Wrong statement 'iotype'", iotype, "in call to IS_CASA_TIME" diff --git a/src/util/datetime_module.f90 b/src/util/datetime_module.f90 index 8c46ca87b..6fb3165bd 100644 --- a/src/util/datetime_module.f90 +++ b/src/util/datetime_module.f90 @@ -1019,6 +1019,23 @@ pure elemental logical function isLeapYear(year) end function isLeapYear + pure function nDeltas(d0, d1, t) + ! Given start and end `datetime` instances `d0` and `d1` and time + ! increment as `timedelta` instance `t`, return the number of `timedelta` + ! instances `t` between d0 and d1. + type(datetime), intent(in) :: d0, d1 + type(timedelta), intent(in) :: t + real(real64) :: datenum0, datenum1, eps, increment + integer :: nDeltas + + eps = 1e-10_real64 + datenum0 = date2num(d0) + datenum1 = date2num(d1) + increment = t%total_seconds() * s2d + nDeltas = floor((datenum1 - datenum0 + eps) / increment) + 1 + end function nDeltas + + pure function datetimeRange(d0, d1, t) ! Given start and end `datetime` instances `d0` and `d1` and time ! increment as `timedelta` instance `t`, returns an array of @@ -1085,6 +1102,36 @@ pure elemental integer function daysInYear(year) end if end function daysInYear + + pure elemental logical function isNewDay(d) + ! Determines whether the given `datetime` `d` is a the start of a month. + type(datetime), intent(in) :: d + + isNewDay = (d%getHour() == 0 .and. d%getMinute() == 0 .and.& + d%getSecond() == 0 .and. d%getMillisecond() == 0) + + end function isNewDay + + + pure elemental logical function isNewMonth(d) + ! Determines whether the given `datetime` `d` is the start of a month. + type(datetime), intent(in) :: d + + isNewYear = (d%getDay() == 1 .and. d%getHour() == 0 .and.& + d%getMinute() == 0 .and. d%getSecond() == 0 .and.& + d%getMillisecond() == 0) + end function isNewYear + + + pure elemental logical function isNewYear(d) + ! Determines whether the given `datetime` `d` is the start of a year. + type(datetime), intent(in) :: d + + isNewYear = (d%getMonth() == 1 .and. d%getDay() == 1 .and.& + d%getHour() == 0 .and. d%getMinute() == 0 .and.& + d%getSecond() == 0 .and. d%getMillisecond() == 0) + end function isNewYear + pure elemental real(real64) function date2num(d) ! Given a datetime instance d, returns number of days since From 823798639450df6335ef161eb9ecb281afe395e0 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Mon, 22 Dec 2025 14:45:23 +1100 Subject: [PATCH 06/28] More instances of specific timings being updated --- src/offline/cable_input.F90 | 6 ++-- src/offline/cable_serial.F90 | 56 +++++++++++++++++++----------------- 2 files changed, 34 insertions(+), 28 deletions(-) diff --git a/src/offline/cable_input.F90 b/src/offline/cable_input.F90 index 84c3b6a66..90b0315bf 100644 --- a/src/offline/cable_input.F90 +++ b/src/offline/cable_input.F90 @@ -301,7 +301,7 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ) IMPLICIT NONE ! Input arguments - REAL, INTENT(OUT) :: dels ! time step size + type(timedelta), intent(out) :: dt REAL, INTENT(IN) :: TFRZ INTEGER, INTENT(INOUT) :: koffset ! offset between met file and desired period INTEGER, INTENT(OUT) :: kend ! number of time steps in simulation @@ -353,6 +353,7 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ) precipTot, & ! used for spinup adj avPrecipInMet ! used for spinup adj CHARACTER(LEN=10) :: todaydate, nowtime ! used to timestamp log file + REAL :: dels ! time step size REAL(4),DIMENSION(1) :: data1 ! temp variable for netcdf reading REAL(4),DIMENSION(1,1) :: data2 ! temp variable for netcdf reading REAL(4), DIMENSION(:), ALLOCATABLE :: temparray1 ! temp read in variable @@ -1597,6 +1598,7 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ) ' some synthesised (as above).' END IF + dt = timedelta(seconds=dels) !=================^^ End met variables search^^======================= END SUBROUTINE open_met_file !============================================================================== @@ -1678,7 +1680,7 @@ SUBROUTINE get_met_data(spinup,spinConv,met,soil,rad, & END SELECT ELSE ! increment hour-of-day by time step size: - met%hod(landpt(i)%cstart) = met%hod(landpt(i)%cstart) + dels/3600.0 + met%hod(landpt(i)%cstart) = met%hod(landpt(i)%cstart) + dt%total_seconds()/3600.0 END IF ! IF(met%hod(landpt(i)%cstart)<0.0) THEN ! may be -ve since longitude diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index a89b8382c..e10fe31cc 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -340,7 +340,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) END IF IF ( RRRR .EQ. 1 ) THEN - CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) + CALL open_met_file( dt, koffset, kend, spinup, CTFRZ ) IF (is_leapyear(YYYY).AND.kend.EQ.2920) THEN STOP 'LEAP YEAR INCOMPATIBILITY WITH INPUT MET !' ENDIF @@ -368,7 +368,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) CASE ('gswp3') ncciy = CurYear - CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) + CALL open_met_file( dt, koffset, kend, spinup, CTFRZ ) CASE ('site') ! site experiment eg AmazonFace (spinup or transient run type) @@ -499,7 +499,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) DO ktau=kstart, kend ts_end = ts_start + dt - WRITE(logn,*) 'Current time: ', ts_start + WRITE(logn,*) 'Current time: ', isoformat(ts_start) WRITE(logn,*) 'Progress -',REAL(ktau)/REAL(kend)*100.0 ! increment total timstep counter @@ -508,7 +508,8 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ! globally (WRT code) accessible kend through USE cable_common_module ktau_gl = ktau_tot - idoy =INT( MOD(REAL(CEILING(REAL((ktau+koffset)/ktauday))),REAL(LOY))) + !idoy =INT( MOD(REAL(CEILING(REAL((ktau+koffset)/ktauday))),REAL(LOY))) + idoy = ts_start%getday() IF ( idoy .EQ. 0 ) idoy = LOY ! needed for CASA-CNP @@ -569,14 +570,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) END SELECT - IF (TRIM(cable_user%MetType).EQ.'' ) THEN - CurYear = met%year(1) - IF ( leaps .AND. IS_LEAPYEAR( CurYear ) ) THEN - LOY = 366 - ELSE - LOY = 365 - ENDIF - ENDIF met%ofsd = met%fsd(:,1) + met%fsd(:,2) canopy%oldcansto=canopy%cansto ! Zero out lai where there is no vegetation acc. to veg. index @@ -584,7 +577,8 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ! At first time step of year, set tile area according to updated LU areas ! and zero casa fluxes - IF (ktau == 1) THEN + if (startofyear(ts_start)) then + !IF (ktau == 1) THEN IF (icycle>1) CALL casa_cnpflux(casaflux,casapool,casabal,.TRUE.) IF ( CABLE_USER%POPLUC) CALL POPLUC_set_patchfrac(POPLUC,LUC_EXPT) ENDIF @@ -616,8 +610,9 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) - ELSE IF ( IS_CASA_TIME("dread", yyyy, ktau, kstart, & - koffset, kend, ktauday, logn) ) THEN ! CLN READ FROM FILE INSTEAD ! + !ELSE IF ( IS_CASA_TIME("dread", yyyy, ktau, kstart, & + !koffset, kend, ktauday, logn) ) THEN ! CLN READ FROM FILE INSTEAD ! + else if (is_casa_time(ts_Start, "dread", logn)) then WRITE(CYEAR,FMT="(I4)")CurYear + INT((ktau-kstart+koffset)/(LOY*ktauday)) ncfile = TRIM(casafile%c2cdumppath)//'c2c_'//CYEAR//'_dump.nc' casa_it = NINT( REAL(ktau / ktauday) ) @@ -637,7 +632,8 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) CABLE_USER%CASA_DUMP_READ, CABLE_USER%CASA_DUMP_WRITE, & LALLOC ) - IF(MOD((ktau-kstart+1),ktauday)==0) THEN + !IF(MOD((ktau-kstart+1),ktauday)==0) THEN + if (isnewday(ts_end)) then !mpidiff ! update time-aggregates of casa pools and fluxes @@ -647,8 +643,9 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ENDIF - IF( ((MOD((ktau-kstart+1),ktauday)==0) .AND. & - MOD((ktau-kstart+1)/ktauday,LOY)==0) )THEN ! end of year + !IF( ((MOD((ktau-kstart+1),`ktauday)==0) .AND. & + !MOD((ktau-kstart+1)/ktauday,LOY)==0) )THEN ! end of year + if (isnewyear(ts_end)) then IF (CABLE_USER%POPLUC) THEN ! Dynamic LUC CALL LUCdriver( casabiome,casapool,casaflux,POP, & @@ -672,23 +669,29 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) IF(icycle >0) THEN - IF ( IS_CASA_TIME("write", yyyy, ktau, kstart, & - koffset, kend, ktauday, logn) ) THEN + !IF ( IS_CASA_TIME("write", yyyy, ktau, kstart, & + !koffset, kend, ktauday, logn) ) THEN + if (is_casa_time(ts_start, "write", logn)) then ctime = ctime +1 !mpidiff CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, & .FALSE. , .TRUE. , count_sum_casa) - CALL WRITE_CASA_OUTPUT_NC (veg, casamet, sum_casapool, casabal, sum_casaflux, & - CASAONLY, ctime, ( ktau.EQ.kend .AND. YYYY .EQ. & - cable_user%YearEnd.AND. RRRR .EQ.NRRRR ) ) + !CALL WRITE_CASA_OUTPUT_NC (veg, casamet, sum_casapool, casabal, sum_casaflux, & + !CASAONLY, ctime, ( ktau.EQ.kend .AND. YYYY .EQ. & + !cable_user%YearEnd.AND. RRRR .EQ.NRRRR ) ) + call write_casa_output_nc(veg, casamet, sum_casapool, casabal,& + sum_casaflux, CASAONLY, ctime, isnewday(ts_start) .and.& + ts_start%getyear() == cable_user%YearEnd .and. RRRR == NRRRR) !mpidiff count_sum_casa = 0 CALL zero_sum_casa(sum_casapool, sum_casaflux) ENDIF - IF (((.NOT.spinup).OR.(spinup.AND.spinConv)).AND. & - MOD((ktau-kstart+1),ktauday)==0) THEN + !IF (((.NOT.spinup).OR.(spinup.AND.spinConv)).AND. & + !MOD((ktau-kstart+1),ktauday)==0) THEN + if (((.not. spinup) .or. (spinup .and. spinConv)) .and.& + isnewday(ts_end)) then IF ( CABLE_USER%CASA_DUMP_WRITE ) THEN SELECT CASE (TRIM(cable_user%MetType)) @@ -789,7 +792,8 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) !$ !$ enddo - IF( ktau == kend ) THEN + !IF( ktau == kend ) THEN + if (isnewyear(ts_end)) then nkend = nkend+1 PRINT *, "NB. Offline-serial runs spinup cycles:", nkend CALL compare_consistency_check_values(new_sumbal) From fb2f4622dba2ae8f9d59032e04807715ed930083 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 10:41:41 +1100 Subject: [PATCH 07/28] Take dt out of cable_input --- src/offline/cable_input.F90 | 6 ++---- src/offline/cable_serial.F90 | 10 ++++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/offline/cable_input.F90 b/src/offline/cable_input.F90 index 90b0315bf..84c3b6a66 100644 --- a/src/offline/cable_input.F90 +++ b/src/offline/cable_input.F90 @@ -301,7 +301,7 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ) IMPLICIT NONE ! Input arguments - type(timedelta), intent(out) :: dt + REAL, INTENT(OUT) :: dels ! time step size REAL, INTENT(IN) :: TFRZ INTEGER, INTENT(INOUT) :: koffset ! offset between met file and desired period INTEGER, INTENT(OUT) :: kend ! number of time steps in simulation @@ -353,7 +353,6 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ) precipTot, & ! used for spinup adj avPrecipInMet ! used for spinup adj CHARACTER(LEN=10) :: todaydate, nowtime ! used to timestamp log file - REAL :: dels ! time step size REAL(4),DIMENSION(1) :: data1 ! temp variable for netcdf reading REAL(4),DIMENSION(1,1) :: data2 ! temp variable for netcdf reading REAL(4), DIMENSION(:), ALLOCATABLE :: temparray1 ! temp read in variable @@ -1598,7 +1597,6 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ) ' some synthesised (as above).' END IF - dt = timedelta(seconds=dels) !=================^^ End met variables search^^======================= END SUBROUTINE open_met_file !============================================================================== @@ -1680,7 +1678,7 @@ SUBROUTINE get_met_data(spinup,spinConv,met,soil,rad, & END SELECT ELSE ! increment hour-of-day by time step size: - met%hod(landpt(i)%cstart) = met%hod(landpt(i)%cstart) + dt%total_seconds()/3600.0 + met%hod(landpt(i)%cstart) = met%hod(landpt(i)%cstart) + dels/3600.0 END IF ! IF(met%hod(landpt(i)%cstart)<0.0) THEN ! may be -ve since longitude diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index e10fe31cc..f3f0f2abd 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -340,8 +340,9 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) END IF IF ( RRRR .EQ. 1 ) THEN - CALL open_met_file( dt, koffset, kend, spinup, CTFRZ ) - IF (is_leapyear(YYYY).AND.kend.EQ.2920) THEN + CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) + dt = timedelta(seconds=dels) + IF (isleapyear(YYYY).AND.kend.EQ.2920) THEN STOP 'LEAP YEAR INCOMPATIBILITY WITH INPUT MET !' ENDIF IF ( NRRRR .GT. 1 ) THEN @@ -368,7 +369,8 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) CASE ('gswp3') ncciy = CurYear - CALL open_met_file( dt, koffset, kend, spinup, CTFRZ ) + CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) + dt = timedelta(seconds=dels) CASE ('site') ! site experiment eg AmazonFace (spinup or transient run type) @@ -577,7 +579,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ! At first time step of year, set tile area according to updated LU areas ! and zero casa fluxes - if (startofyear(ts_start)) then + if (isnewyear(ts_start)) then !IF (ktau == 1) THEN IF (icycle>1) CALL casa_cnpflux(casaflux,casapool,casabal,.TRUE.) IF ( CABLE_USER%POPLUC) CALL POPLUC_set_patchfrac(POPLUC,LUC_EXPT) From 2a06a3d4eb4b7cc1282841ef1db155b53dfcacf4 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 10:44:31 +1100 Subject: [PATCH 08/28] Fix some typos --- src/util/datetime_module.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util/datetime_module.f90 b/src/util/datetime_module.f90 index 6fb3165bd..8990ea5d1 100644 --- a/src/util/datetime_module.f90 +++ b/src/util/datetime_module.f90 @@ -1117,10 +1117,10 @@ pure elemental logical function isNewMonth(d) ! Determines whether the given `datetime` `d` is the start of a month. type(datetime), intent(in) :: d - isNewYear = (d%getDay() == 1 .and. d%getHour() == 0 .and.& + isNewMonth = (d%getDay() == 1 .and. d%getHour() == 0 .and.& d%getMinute() == 0 .and. d%getSecond() == 0 .and.& d%getMillisecond() == 0) - end function isNewYear + end function isNewMonth pure elemental logical function isNewYear(d) From 6c73c32e189e5c6aac7aa570c237b877cce2bac0 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 10:56:55 +1100 Subject: [PATCH 09/28] try removing seemingly unused nday and ndays --- src/offline/CASAONLY_LUC.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index 62103fba6..fcab5cfd4 100644 --- a/src/offline/CASAONLY_LUC.F90 +++ b/src/offline/CASAONLY_LUC.F90 @@ -66,8 +66,7 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & INTEGER :: myearspin,nyear, yyyy, nyear_dump CHARACTER(LEN=99) :: ncfile CHARACTER(LEN=4) :: cyear - INTEGER :: ktau,ktauday,nday,idoy,ktaux,ktauy,nloop - INTEGER, SAVE :: ndays + INTEGER :: ktau,ktauday,idoy,ktaux,ktauy,nloop REAL, DIMENSION(mp) :: cleaf2met, cleaf2str, croot2met, croot2str, cwood2cwd REAL, DIMENSION(mp) :: nleaf2met, nleaf2str, nroot2met, nroot2str, nwood2cwd REAL, DIMENSION(mp) :: pleaf2met, pleaf2str, proot2met, proot2str, pwood2cwd @@ -98,7 +97,6 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & ENDIF ktauday=INT(24.0*3600.0/dels) - nday=(kend-kstart+1)/ktauday ctime = 0 CALL zero_sum_casa(sum_casapool, sum_casaflux) count_sum_casa = 0 From dd14d057438801fc0974f89114d10212f5aee108 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 11:00:38 +1100 Subject: [PATCH 10/28] Temporarily reintroduce old IS_CASA_TIME function --- src/shared/casa_ncdf.F90 | 76 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/src/shared/casa_ncdf.F90 b/src/shared/casa_ncdf.F90 index 0bc47b857..d4f1b242e 100644 --- a/src/shared/casa_ncdf.F90 +++ b/src/shared/casa_ncdf.F90 @@ -444,7 +444,83 @@ FUNCTION IS_CASA_TIME(currtime, iotype, logn) END FUNCTION IS_CASA_TIME + FUNCTION IS_CASA_TIME(iotype, yyyy, ktau, kstart, koffset, kend, ktauday, logn) + USE cable_common_module, ONLY: CABLE_USER + ! Correctly determine if it is time to dump-read or standard-write + ! casa output from cable_serial. + ! Writing casa-dump data is handled in casa_cable and therefore not \ + ! captured here + !cable_common module was intended to be unequivocally common to all + !applications. iovars is an offline module and so not appropriate to include + !here. Suggested FIX is to move decs of vars needed (e.g. leaps) to here, and + !then use common in iovars +#ifdef Vanessas_common + USE cable_IO_vars_module, ONLY: leaps +#endif + IMPLICIT NONE + + + LOGICAL :: IS_CASA_TIME + INTEGER ,INTENT(IN) :: yyyy, ktau, kstart, koffset, kend, ktauday, logn + CHARACTER,INTENT(IN) :: iotype*5 + LOGICAL :: is_eod, is_eom, is_eoy + INTEGER :: doy, m + INTEGER, DIMENSION(12) :: MONTH + + + is_eom = .FALSE. + is_eoy = .FALSE. + IS_CASA_TIME = .FALSE. + + + MONTH = (/ 31,28,31,30,31,30,31,31,30,31,30,31 /) + is_eod = ( MOD((ktau-kstart+1+koffset),ktauday).EQ.0 ) + IF ( .NOT. is_eod ) RETURN ! NO if it is not end of day + + +#ifdef Vanessas_common + IF ( IS_LEAPYEAR( YYYY ) .AND. leaps ) THEN + MONTH(2) = 29 + ELSE + MONTH(2) = 28 + ENDIF +#endif + + + ! Check for reading from dump-file (hard-wired to daily casa-timestep) + IF ( iotype .EQ. "dread" ) THEN + IF ( CABLE_USER%CASA_DUMP_READ ) IS_CASA_TIME = .TRUE. + ! Check for writing of casa dump output + ELSE IF ( iotype .EQ. "dwrit" ) THEN + IF ( CABLE_USER%CASA_DUMP_WRITE ) IS_CASA_TIME = .TRUE. + ! Check for writing of casa standard output + ELSE IF ( iotype .EQ. "write" ) THEN + + + doy = NINT(REAL(ktau-kstart+1+koffset)/REAL(ktauday)) + DO m = 1, 12 + IF ( doy .EQ. SUM(MONTH(1:m)) ) THEN + is_eom = .TRUE. + IF ( m .EQ. 12 ) is_eoy = .TRUE. + EXIT + ENDIF + END DO + + + SELECT CASE ( TRIM(CABLE_USER%CASA_OUT_FREQ) ) + CASE ("daily" ) ; IS_CASA_TIME = .TRUE. + CASE ("monthly" ) ; IF ( is_eom ) IS_CASA_TIME = .TRUE. + CASE ("annually") ; IF ( is_eoy ) IS_CASA_TIME = .TRUE. + END SELECT + ELSE + WRITE(logn,*)"Wrong statement 'iotype'", iotype, "in call to IS_CASA_TIME" + WRITE(* ,*)"Wrong statement 'iotype'", iotype, "in call to IS_CASA_TIME" + STOP -1 + ENDIF + + + END FUNCTION IS_CASA_TIME END MODULE casa_ncdf_module From 7677535d128253e4567179f0c7f025a0be28f9fa Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 11:17:00 +1100 Subject: [PATCH 11/28] Change CASAONLY_LUC version of IS_CASA_TIME --- src/offline/CASAONLY_LUC.F90 | 17 ++++---- src/offline/cable_serial.F90 | 2 +- src/shared/casa_ncdf.F90 | 78 ------------------------------------ 3 files changed, 8 insertions(+), 89 deletions(-) diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index fcab5cfd4..b7f800932 100644 --- a/src/offline/CASAONLY_LUC.F90 +++ b/src/offline/CASAONLY_LUC.F90 @@ -1,4 +1,4 @@ -SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & +SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, & casaflux,casamet,casabal,phen,POP,climate,LALLOC,LUC_EXPT, POPLUC, & sum_casapool, sum_casaflux ) @@ -24,14 +24,13 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & USE casa_cable USE casa_inout_module USE biogeochem_mod, ONLY : biogeochem -USE casa_offline_inout_module, ONLY : WRITE_CASA_OUTPUT_NC + USE casa_offline_inout_module, ONLY : WRITE_CASA_OUTPUT_NC + use datetime_module, only: datetime IMPLICIT NONE !!CLN CHARACTER(LEN=99), INTENT(IN) :: fcnpspin - REAL, INTENT(IN) :: dels - INTEGER, INTENT(IN) :: kstart - INTEGER, INTENT(IN) :: kend + REAL, INTENT(IN) :: dels, ktauday INTEGER, INTENT(IN) :: LALLOC TYPE (veg_parameter_type), INTENT(INOUT) :: veg ! vegetation parameters TYPE (soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters @@ -48,6 +47,7 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & TYPE (casa_pool) , INTENT(INOUT) :: sum_casapool TYPE (casa_flux) , INTENT(INOUT) :: sum_casaflux + type(datetime), intent(in) :: curr_time TYPE (casa_met) :: casaspin @@ -65,8 +65,7 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & ! local variables INTEGER :: myearspin,nyear, yyyy, nyear_dump CHARACTER(LEN=99) :: ncfile - CHARACTER(LEN=4) :: cyear - INTEGER :: ktau,ktauday,idoy,ktaux,ktauy,nloop + INTEGER :: ktau,ktauday,idoy REAL, DIMENSION(mp) :: cleaf2met, cleaf2str, croot2met, croot2str, cwood2cwd REAL, DIMENSION(mp) :: nleaf2met, nleaf2str, nroot2met, nroot2str, nwood2cwd REAL, DIMENSION(mp) :: pleaf2met, pleaf2str, proot2met, proot2str, pwood2cwd @@ -105,7 +104,6 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & myearspin = CABLE_USER%YEAREND - CABLE_USER%YEARSTART + 1 yyyy = CABLE_USER%YEARSTART - 1 - DO nyear=1,myearspin yyyy = yyyy+1 WRITE(*,*) 'casaonly_LUC', YYYY @@ -268,8 +266,7 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & ENDIF ! end of year - IF ( IS_CASA_TIME("write", yyyy, ktau, kstart, & - 0, kend, ktauday, logn) ) THEN + IF ( IS_CASA_TIME(ts_start, "write", logn) ) THEN ctime = ctime +1 CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, & diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index f3f0f2abd..9dbe8a589 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -468,7 +468,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ELSEIF ( casaonly .AND. (.NOT. spincasa) ) THEN !.AND. cable_user%popluc) THEN - CALL CASAONLY_LUC(dels,kstart,kend,veg,soil,casabiome,casapool, & + CALL CASAONLY_LUC(dels, ts_start, ktauday ,veg,soil,casabiome,casapool, & casaflux,casamet,casabal,phen,POP,climate,LALLOC, LUC_EXPT, POPLUC, & sum_casapool, sum_casaflux) SPINon = .FALSE. diff --git a/src/shared/casa_ncdf.F90 b/src/shared/casa_ncdf.F90 index d4f1b242e..357e7d4aa 100644 --- a/src/shared/casa_ncdf.F90 +++ b/src/shared/casa_ncdf.F90 @@ -444,84 +444,6 @@ FUNCTION IS_CASA_TIME(currtime, iotype, logn) END FUNCTION IS_CASA_TIME - FUNCTION IS_CASA_TIME(iotype, yyyy, ktau, kstart, koffset, kend, ktauday, logn) - - USE cable_common_module, ONLY: CABLE_USER - ! Correctly determine if it is time to dump-read or standard-write - ! casa output from cable_serial. - ! Writing casa-dump data is handled in casa_cable and therefore not \ - ! captured here - !cable_common module was intended to be unequivocally common to all - !applications. iovars is an offline module and so not appropriate to include - !here. Suggested FIX is to move decs of vars needed (e.g. leaps) to here, and - !then use common in iovars -#ifdef Vanessas_common - USE cable_IO_vars_module, ONLY: leaps -#endif - IMPLICIT NONE - - - LOGICAL :: IS_CASA_TIME - INTEGER ,INTENT(IN) :: yyyy, ktau, kstart, koffset, kend, ktauday, logn - CHARACTER,INTENT(IN) :: iotype*5 - LOGICAL :: is_eod, is_eom, is_eoy - INTEGER :: doy, m - INTEGER, DIMENSION(12) :: MONTH - - - is_eom = .FALSE. - is_eoy = .FALSE. - IS_CASA_TIME = .FALSE. - - - MONTH = (/ 31,28,31,30,31,30,31,31,30,31,30,31 /) - is_eod = ( MOD((ktau-kstart+1+koffset),ktauday).EQ.0 ) - IF ( .NOT. is_eod ) RETURN ! NO if it is not end of day - - -#ifdef Vanessas_common - IF ( IS_LEAPYEAR( YYYY ) .AND. leaps ) THEN - MONTH(2) = 29 - ELSE - MONTH(2) = 28 - ENDIF -#endif - - - ! Check for reading from dump-file (hard-wired to daily casa-timestep) - IF ( iotype .EQ. "dread" ) THEN - IF ( CABLE_USER%CASA_DUMP_READ ) IS_CASA_TIME = .TRUE. - ! Check for writing of casa dump output - ELSE IF ( iotype .EQ. "dwrit" ) THEN - IF ( CABLE_USER%CASA_DUMP_WRITE ) IS_CASA_TIME = .TRUE. - ! Check for writing of casa standard output - ELSE IF ( iotype .EQ. "write" ) THEN - - - doy = NINT(REAL(ktau-kstart+1+koffset)/REAL(ktauday)) - DO m = 1, 12 - IF ( doy .EQ. SUM(MONTH(1:m)) ) THEN - is_eom = .TRUE. - IF ( m .EQ. 12 ) is_eoy = .TRUE. - EXIT - ENDIF - END DO - - - SELECT CASE ( TRIM(CABLE_USER%CASA_OUT_FREQ) ) - CASE ("daily" ) ; IS_CASA_TIME = .TRUE. - CASE ("monthly" ) ; IF ( is_eom ) IS_CASA_TIME = .TRUE. - CASE ("annually") ; IF ( is_eoy ) IS_CASA_TIME = .TRUE. - END SELECT - ELSE - WRITE(logn,*)"Wrong statement 'iotype'", iotype, "in call to IS_CASA_TIME" - WRITE(* ,*)"Wrong statement 'iotype'", iotype, "in call to IS_CASA_TIME" - STOP -1 - ENDIF - - - END FUNCTION IS_CASA_TIME - END MODULE casa_ncdf_module From 6915dd3a48f07648cfec4f271b8cc34f202af870 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 11:29:49 +1100 Subject: [PATCH 12/28] Remove doubly defined variable and re-add required variable --- src/offline/CASAONLY_LUC.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index b7f800932..7facb7cfa 100644 --- a/src/offline/CASAONLY_LUC.F90 +++ b/src/offline/CASAONLY_LUC.F90 @@ -63,9 +63,9 @@ SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, REAL(r_2), DIMENSION(:), ALLOCATABLE, SAVE :: avg_xnplimit, avg_xkNlimiting,avg_xklitter, avg_xksoil ! local variables - INTEGER :: myearspin,nyear, yyyy, nyear_dump + INTEGER :: myearspin,nyear, yyyy, nyear_dump, cyear CHARACTER(LEN=99) :: ncfile - INTEGER :: ktau,ktauday,idoy + INTEGER :: ktau,idoy REAL, DIMENSION(mp) :: cleaf2met, cleaf2str, croot2met, croot2str, cwood2cwd REAL, DIMENSION(mp) :: nleaf2met, nleaf2str, nroot2met, nroot2str, nwood2cwd REAL, DIMENSION(mp) :: pleaf2met, pleaf2str, proot2met, proot2str, pwood2cwd From 2d6443b19e7a269a6480d89da133c8965f4f2ba8 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 11:37:26 +1100 Subject: [PATCH 13/28] Remove unused ktau variable from biogeochem --- src/offline/CASAONLY_LUC.F90 | 2 +- src/offline/cable_mpimaster.F90 | 2 +- src/offline/cable_mpiworker.F90 | 6 +++--- src/offline/spincasacnp.F90 | 4 ++-- src/science/casa-cnp/bgcdriver.F90 | 4 ++-- src/science/casa-cnp/biogeochem_casa.F90 | 2 +- 6 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index 7facb7cfa..1d00ce0aa 100644 --- a/src/offline/CASAONLY_LUC.F90 +++ b/src/offline/CASAONLY_LUC.F90 @@ -154,7 +154,7 @@ SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, climate%qtemp_max_last_year(:) = casamet%mtempspin(:,idoy) - CALL biogeochem(ktau,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter, & xksoil,xkleaf,xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 591caba01..a74fc484d 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -8265,7 +8265,7 @@ SUBROUTINE master_CASAONLY_LUC( dels,kstart,kend,veg,casabiome,casapool, & !$ ! zero annual sums !$ if (idoy==1) CALL casa_cnpflux(casaflux,casapool,casabal,.TRUE.) -!$ CALL biogeochem(ktau,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & +!$ CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & !$ casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter, & !$ xksoil,xkleaf,xkleafcold,xkleafdry,& !$ cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index 030f62af1..dd9567670 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -7042,7 +7042,7 @@ SUBROUTINE worker_spincasacnp( dels,kstart,kend,mloop,veg,soil,casabiome,casapoo DO idoy=1,mdyear ktau=(idoy-1)*ktauday +1 CALL MPI_Recv (MPI_BOTTOM, 1, casa_dump_t, 0, idoy, icomm, stat, ierr) - CALL biogeochem(ktau,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter, & xksoil,xkleaf,xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & @@ -7210,7 +7210,7 @@ SUBROUTINE worker_spincasacnp( dels,kstart,kend,mloop,veg,soil,casabiome,casapoo ktau=(idoy-1)*ktauday +1 CALL MPI_Recv (MPI_BOTTOM, 1, casa_dump_t, 0, idoy, icomm, stat, ierr) - CALL biogeochem(ktauy,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter,xksoil,xkleaf,& xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & @@ -7332,7 +7332,7 @@ SUBROUTINE worker_CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & CALL MPI_Recv (MPI_BOTTOM, 1, casa_dump_t, 0, idoy, icomm, stat, ierr) - CALL biogeochem(ktau,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter, & xksoil,xkleaf,xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & diff --git a/src/offline/spincasacnp.F90 b/src/offline/spincasacnp.F90 index 72f0d3bfc..374c24ba9 100755 --- a/src/offline/spincasacnp.F90 +++ b/src/offline/spincasacnp.F90 @@ -173,7 +173,7 @@ SUBROUTINE spincasacnp( dels,kstart,kend,mloop,veg,soil,casabiome,casapool, & ENDIF ! write(6699,*) casaflux%cgpp(1), climate%mtemp(1), casaflux%crmplant(1,1) - CALL biogeochem(ktau,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter, & xksoil,xkleaf,xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & @@ -376,7 +376,7 @@ SUBROUTINE spincasacnp( dels,kstart,kend,mloop,veg,soil,casabiome,casapool, & - CALL biogeochem(ktauy,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter,xksoil,xkleaf,& xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & diff --git a/src/science/casa-cnp/bgcdriver.F90 b/src/science/casa-cnp/bgcdriver.F90 index 34cc6ddae..103d01e7d 100644 --- a/src/science/casa-cnp/bgcdriver.F90 +++ b/src/science/casa-cnp/bgcdriver.F90 @@ -105,7 +105,7 @@ SUBROUTINE bgcdriver(ktau,kstart,kend,dels,met,ssnow,canopy,veg,soil, & ENDIF # endif - CALL biogeochem(ktau,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate, xnplimit,xkNlimiting,xklitter,xksoil, & xkleaf,xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & @@ -146,7 +146,7 @@ SUBROUTINE bgcdriver(ktau,kstart,kend,dels,met,ssnow,canopy,veg,soil, & IF( MOD((ktau-kstart+1),ktauday) == 0 ) THEN ! end of day - CALL biogeochem(ktau,dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & + CALL biogeochem(dels,idoy,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter,xksoil,xkleaf, & xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & diff --git a/src/science/casa-cnp/biogeochem_casa.F90 b/src/science/casa-cnp/biogeochem_casa.F90 index ef4f860b6..7bd9c2c63 100644 --- a/src/science/casa-cnp/biogeochem_casa.F90 +++ b/src/science/casa-cnp/biogeochem_casa.F90 @@ -4,7 +4,7 @@ MODULE biogeochem_mod CONTAINS - SUBROUTINE biogeochem(ktau,dels,idoY,LALLOC,veg,soil,casabiome,casapool,casaflux, & + SUBROUTINE biogeochem(dels,idoY,LALLOC,veg,soil,casabiome,casapool,casaflux, & casamet,casabal,phen,POP,climate,xnplimit,xkNlimiting,xklitter,xksoil,xkleaf,xkleafcold,xkleafdry,& cleaf2met,cleaf2str,croot2met,croot2str,cwood2cwd, & nleaf2met,nleaf2str,nroot2met,nroot2str,nwood2cwd, & From a10efb5aca41af6b227ad7bdccfa6f325d95b005 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Thu, 8 Jan 2026 11:40:12 +1100 Subject: [PATCH 14/28] Remove ktau from biogeochem body --- src/science/casa-cnp/biogeochem_casa.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/science/casa-cnp/biogeochem_casa.F90 b/src/science/casa-cnp/biogeochem_casa.F90 index 7bd9c2c63..b29a7642d 100644 --- a/src/science/casa-cnp/biogeochem_casa.F90 +++ b/src/science/casa-cnp/biogeochem_casa.F90 @@ -18,7 +18,6 @@ SUBROUTINE biogeochem(dels,idoY,LALLOC,veg,soil,casabiome,casapool,casaflux, & USE casa_rplant_module, ONLY: casa_rplant IMPLICIT NONE -INTEGER, INTENT(IN) :: ktau REAL, INTENT(IN) :: dels INTEGER, INTENT(IN) :: idoy INTEGER, INTENT(IN) :: LALLOC From 4e43da68f96e15fb4246ed90b770750e41169578 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 10:55:15 +1100 Subject: [PATCH 15/28] Reorder function args to be more logical --- src/offline/CASAONLY_LUC.F90 | 11 +++++------ src/offline/cable_serial.F90 | 11 +---------- src/shared/casa_ncdf.F90 | 10 +++++----- 3 files changed, 11 insertions(+), 21 deletions(-) diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index 1d00ce0aa..5e742cb27 100644 --- a/src/offline/CASAONLY_LUC.F90 +++ b/src/offline/CASAONLY_LUC.F90 @@ -1,4 +1,4 @@ -SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, & +SUBROUTINE CASAONLY_LUC(dels, curr_time, ktauday, veg,soil,casabiome,casapool, & casaflux,casamet,casabal,phen,POP,climate,LALLOC,LUC_EXPT, POPLUC, & sum_casapool, sum_casaflux ) @@ -63,9 +63,10 @@ SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, REAL(r_2), DIMENSION(:), ALLOCATABLE, SAVE :: avg_xnplimit, avg_xkNlimiting,avg_xklitter, avg_xksoil ! local variables - INTEGER :: myearspin,nyear, yyyy, nyear_dump, cyear + INTEGER :: myearspin,nyear, yyyy, nyear_dump CHARACTER(LEN=99) :: ncfile - INTEGER :: ktau,idoy + CHARACTER(LEN=4) :: cyear + INTEGER :: idoy REAL, DIMENSION(mp) :: cleaf2met, cleaf2str, croot2met, croot2str, cwood2cwd REAL, DIMENSION(mp) :: nleaf2met, nleaf2str, nroot2met, nroot2str, nwood2cwd REAL, DIMENSION(mp) :: pleaf2met, pleaf2str, proot2met, proot2str, pwood2cwd @@ -95,7 +96,6 @@ SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, Iw = POP%Iwood ENDIF - ktauday=INT(24.0*3600.0/dels) ctime = 0 CALL zero_sum_casa(sum_casapool, sum_casaflux) count_sum_casa = 0 @@ -127,7 +127,6 @@ SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, CALL read_casa_dump( ncfile,casamet, casaflux, phen,climate, 1,1,.TRUE. ) !!CLN901 format(A99) DO idoy=1,mdyear - ktau=(idoy-1)*ktauday +ktauday casamet%tairk(:) = casamet%Tairkspin(:,idoy) casamet%tsoil(:,1) = casamet%Tsoilspin_1(:,idoy) @@ -266,7 +265,7 @@ SUBROUTINE CASAONLY_LUC( dels, curr_time, ktauday, veg,soil,casabiome,casapool, ENDIF ! end of year - IF ( IS_CASA_TIME(ts_start, "write", logn) ) THEN + IF ( IS_CASA_TIME("write", curr_time, logn) ) THEN ctime = ctime +1 CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, & diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 9dbe8a589..0a75bcc85 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -608,12 +608,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ssnow%rnof2 = ssnow%rnof2*dels ssnow%runoff = ssnow%runoff*dels - - - - - !ELSE IF ( IS_CASA_TIME("dread", yyyy, ktau, kstart, & - !koffset, kend, ktauday, logn) ) THEN ! CLN READ FROM FILE INSTEAD ! else if (is_casa_time(ts_Start, "dread", logn)) then WRITE(CYEAR,FMT="(I4)")CurYear + INT((ktau-kstart+koffset)/(LOY*ktauday)) ncfile = TRIM(casafile%c2cdumppath)//'c2c_'//CYEAR//'_dump.nc' @@ -670,10 +664,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ! WRITE CASA OUTPUT IF(icycle >0) THEN - - !IF ( IS_CASA_TIME("write", yyyy, ktau, kstart, & - !koffset, kend, ktauday, logn) ) THEN - if (is_casa_time(ts_start, "write", logn)) then + if (is_casa_time("write", ts_start, logn)) then ctime = ctime +1 !mpidiff CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, & diff --git a/src/shared/casa_ncdf.F90 b/src/shared/casa_ncdf.F90 index 357e7d4aa..e9a80e674 100644 --- a/src/shared/casa_ncdf.F90 +++ b/src/shared/casa_ncdf.F90 @@ -404,7 +404,7 @@ SUBROUTINE DOYSOD2YMDHMS( YYYY,DOY,SOD,MM,DD,HOUR,MINUTE,SECOND ) END SUBROUTINE DOYSOD2YMDHMS - FUNCTION IS_CASA_TIME(currtime, iotype, logn) + FUNCTION IS_CASA_TIME(iotype, curr_time, logn) USE cable_common_module, ONLY: CABLE_USER ! Correctly determine if it is time to dump-read or standard-write @@ -418,7 +418,7 @@ FUNCTION IS_CASA_TIME(currtime, iotype, logn) IMPLICIT NONE LOGICAL :: IS_CASA_TIME - type(datetime), intent(in) :: currtime + type(datetime), intent(in) :: curr_time INTEGER ,INTENT(IN) :: logn CHARACTER,INTENT(IN) :: iotype*5 @@ -432,9 +432,9 @@ FUNCTION IS_CASA_TIME(currtime, iotype, logn) ELSE IF ( iotype .EQ. "write" ) THEN SELECT CASE ( TRIM(CABLE_USER%CASA_OUT_FREQ) ) - CASE ("daily" ) ; IS_CASA_TIME = isnewday(currtime) - CASE ("monthly" ) ; IS_CASA_TIME = isnewmonth(currtime) - CASE ("annually") ; IS_CASA_TIME = isnewyear(currtime) + CASE ("daily" ) ; IS_CASA_TIME = isnewday(curr_time) + CASE ("monthly" ) ; IS_CASA_TIME = isnewmonth(curr_time) + CASE ("annually") ; IS_CASA_TIME = isnewyear(curr_time) END SELECT ELSE WRITE(logn,*)"Wrong statement 'iotype'", iotype, "in call to IS_CASA_TIME" From ee020afb21c3194724d216f9ad672e3e3ffb974e Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 11:07:55 +1100 Subject: [PATCH 16/28] Fix data type of ktauday --- src/offline/CASAONLY_LUC.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index 5e742cb27..6145a17f9 100644 --- a/src/offline/CASAONLY_LUC.F90 +++ b/src/offline/CASAONLY_LUC.F90 @@ -30,8 +30,8 @@ SUBROUTINE CASAONLY_LUC(dels, curr_time, ktauday, veg,soil,casabiome,casapool, & IMPLICIT NONE !!CLN CHARACTER(LEN=99), INTENT(IN) :: fcnpspin - REAL, INTENT(IN) :: dels, ktauday - INTEGER, INTENT(IN) :: LALLOC + REAL, INTENT(IN) :: dels + INTEGER, INTENT(IN) :: LALLOC, ktauday TYPE (veg_parameter_type), INTENT(INOUT) :: veg ! vegetation parameters TYPE (soil_parameter_type), INTENT(INOUT) :: soil ! soil parameters TYPE (casa_biome), INTENT(INOUT) :: casabiome From 8470c00d29cbf790ee3c06bd6725d903f9431635 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 11:21:58 +1100 Subject: [PATCH 17/28] Re-add accidentally removed line --- src/offline/CASAONLY_LUC.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index 6145a17f9..704b69686 100644 --- a/src/offline/CASAONLY_LUC.F90 +++ b/src/offline/CASAONLY_LUC.F90 @@ -104,6 +104,7 @@ SUBROUTINE CASAONLY_LUC(dels, curr_time, ktauday, veg,soil,casabiome,casapool, & myearspin = CABLE_USER%YEAREND - CABLE_USER%YEARSTART + 1 yyyy = CABLE_USER%YEARSTART - 1 + DO nyear=1,myearspin yyyy = yyyy+1 WRITE(*,*) 'casaonly_LUC', YYYY From 6ef1d3f1b9baed1fc1d7d1c548ee5c5a82de60b1 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 11:27:42 +1100 Subject: [PATCH 18/28] Fix is_casa_time args and move dt definition --- src/offline/cable_serial.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 0a75bcc85..0cbd8ab16 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -273,7 +273,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) type(datetime) :: ts_start, ts_end type(timedelta) :: dt - dt = timedelta(seconds=int(dels)) ktauday = nint(d2s / dels) ! Set the calendar- some forcing cases have special rules @@ -341,7 +340,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) IF ( RRRR .EQ. 1 ) THEN CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) - dt = timedelta(seconds=dels) IF (isleapyear(YYYY).AND.kend.EQ.2920) THEN STOP 'LEAP YEAR INCOMPATIBILITY WITH INPUT MET !' ENDIF @@ -370,7 +368,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) CASE ('gswp3') ncciy = CurYear CALL open_met_file( dels, koffset, kend, spinup, CTFRZ ) - dt = timedelta(seconds=dels) CASE ('site') ! site experiment eg AmazonFace (spinup or transient run type) @@ -398,6 +395,9 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) END SELECT + ! Set the timestep here, after any met specific modifications of dels + dt = timedelta(second=int(dels)) + ! Checks where parameters and initialisations should be loaded from. ! If they can be found in either the met file or restart file, they will ! load from there, with the met file taking precedence. Otherwise, they'll @@ -608,7 +608,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ssnow%rnof2 = ssnow%rnof2*dels ssnow%runoff = ssnow%runoff*dels - else if (is_casa_time(ts_Start, "dread", logn)) then + else if (is_casa_time("dread", ts_start, logn)) then WRITE(CYEAR,FMT="(I4)")CurYear + INT((ktau-kstart+koffset)/(LOY*ktauday)) ncfile = TRIM(casafile%c2cdumppath)//'c2c_'//CYEAR//'_dump.nc' casa_it = NINT( REAL(ktau / ktauday) ) From 0d8a67d7ac863b326bd890c601840f7b38c1e785 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 12:06:43 +1100 Subject: [PATCH 19/28] missing s in keyword arg --- src/offline/cable_serial.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 0cbd8ab16..d6e73ef85 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -396,7 +396,7 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) END SELECT ! Set the timestep here, after any met specific modifications of dels - dt = timedelta(second=int(dels)) + dt = timedelta(seconds=int(dels)) ! Checks where parameters and initialisations should be loaded from. ! If they can be found in either the met file or restart file, they will From 4a4c19a48f8a458df13d8e14967ae6c71aa7ab6b Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 14:51:56 +1100 Subject: [PATCH 20/28] Try to get the MPI building as well, and fix ts_start not being updated --- src/offline/cable_mpimaster.F90 | 80 +++++++++++++++------------------ src/offline/cable_mpiworker.F90 | 71 +++++++++++++++-------------- src/offline/cable_serial.F90 | 3 ++ 3 files changed, 78 insertions(+), 76 deletions(-) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index a74fc484d..bc6b2c467 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -224,6 +224,8 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) USE landuse_variable USE bgcdriver_mod, ONLY : bgcdriver USE casa_offline_inout_module, ONLY : WRITE_CASA_RESTART_NC, WRITE_CASA_OUTPUT_NC + + use datetime_module IMPLICIT NONE ! MPI: @@ -329,6 +331,29 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) integer, dimension(:), allocatable, save :: cstart,cend,nap real(r_2), dimension(:,:,:), allocatable, save :: patchfrac_new + type(datetime) :: ts_start, ts_end + type(timedelta) :: dt + + ! Set the calendar- some forcing cases have special rules + select case (trim(cable_user%mettype)) + case ("plum") + if (plume%leapyears) then + call setcalendar("gregorian") + else + call setcalendar("noleaps") + endif + + case ("cru") + call setcalendar("noleaps") + + case default + if (leaps) then + call setcalendar("gregorian") + else + call setcalendar("noleaps") + endif + + end select ! END header @@ -340,16 +365,12 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) CurYear = YYYY - !ccc Set calendar attribute: dependant on the value of `leaps` - ! dependant on the MetType and set during initialisation. - calendar = "noleap" - LOY = 365 - IF ( leaps ) THEN - calendar = "standard" - ENDIF - IF ( leaps .AND. IS_LEAPYEAR( YYYY ) ) THEN - LOY = 366 - ENDIF + ts_start = datetime(year=YYYY, month=1, day=1) + + LOY = daysInYear(YYYY) + + ! Set kend for all cases bar princeton and gswp + kend = ktauday * LOY SELECT CASE (TRIM(cable_user%MetType)) CASE ('gswp', 'prin') @@ -458,6 +479,8 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) ! file themselves CALL MPI_Bcast (dels, 1, MPI_REAL, 0, comm, ierr) + dt = timedelta(seconds=int(dels)) + ! MPI: need to know extents before creating datatypes CALL find_extents @@ -673,14 +696,6 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) ENDIF - ! IF ( CASAONLY .AND. IS_CASA_TIME("dread", yyyy, iktau, kstart, koffset, & - ! kend, ktauday, logn) ) THEN - ! WRITE(CYEAR,FMT="(I4)")CurYear + INT((ktau-kstart+koffset)/(LOY*ktauday)) - ! ncfile = TRIM(casafile%c2cdumppath)//'c2c_'//CYEAR//'_dump.nc' - ! casa_it = NINT( REAL(iktau / ktauday) ) - ! CALL read_casa_dump( ncfile, casamet, casaflux,phen, casa_it, kend, .FALSE. ) - ! ENDIF - canopy%oldcansto=canopy%cansto ! Zero out lai where there is no vegetation acc. to veg. index @@ -710,6 +725,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) iktau = iktau + 1 oktau = oktau + 1 + ts_end = ts_start + dt WRITE(logn,*) 'Progress -',REAL(ktau)/REAL(kend)*100.0 met%year = imet%year @@ -748,15 +764,6 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) (TRIM(cable_user%MetType) .NE. 'gswp3') .AND. & (TRIM(cable_user%MetType) .NE. 'prin' )) CurYear = met%year(1) -!$ IF ( CASAONLY .AND. IS_CASA_TIME("dread", yyyy, iktau, kstart, koffset, & -!$ kend, ktauday, logn) ) THEN -!$ ! CLN READ FROM FILE INSTEAD ! -!$ WRITE(CYEAR,FMT="(I4)")CurYear + INT((ktau-kstart+koffset)/(LOY*ktauday)) -!$ ncfile = TRIM(casafile%c2cdumppath)//'c2c_'//CYEAR//'_dump.nc' -!$ casa_it = NINT( REAL(iktau / ktauday) ) -!$ CALL read_casa_dump( ncfile, casamet, casaflux, casa_it, kend, .FALSE. ) -!$ ENDIF - ! Zero out lai where there is no vegetation acc. to veg. index WHERE ( iveg%iveg(:) .GE. 14 ) iveg%vlai = 0. @@ -775,8 +782,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) CALL MPI_Waitall (wnp, recv_req, recv_stats, ierr) ! receive casa dump requirements from worker IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND. & - ( IS_CASA_TIME("dwrit", yyyy, oktau, kstart, & - koffset, kend, ktauday, logn) ) ) THEN + ( IS_CASA_TIME("dwrit", ts_start, logn) ) ) THEN CALL master_receive ( ocomm, oktau, casa_dump_ts ) ! CALL MPI_Waitall (wnp, recv_req, recv_stats, ierr) @@ -792,16 +798,6 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) CALL master_send_input (icomm, inp_ts, iktau) ! CALL MPI_Waitall (wnp, inp_req, inp_stats, ierr) -!$ IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND. & -!$ ( IS_CASA_TIME("dwrit", yyyy, oktau, kstart, & -!$ koffset, kend, ktauday, logn) ) ) THEN -!$ WRITE(CYEAR,FMT="(I4)") CurYear + INT((ktau-kstart)/(LOY*ktauday)) -!$ ncfile = TRIM(casafile%c2cdumppath)//'c2c_'//CYEAR//'_dump.nc' -!$ CALL write_casa_dump( ncfile, casamet , casaflux, idoy, & -!$ kend/ktauday ) -!$ -!$ ENDIF - IF (((.NOT.spinup).OR.(spinup.AND.spinConv)).AND. & MOD((ktau-kstart+1),ktauday)==0) THEN IF ( CABLE_USER%CASA_DUMP_WRITE ) THEN @@ -838,8 +834,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) IF ((.NOT.spinup).OR.(spinup.AND.spinConv)) THEN IF (icycle >0) THEN - IF ( IS_CASA_TIME("write", yyyy, oktau, kstart, & - koffset, kend, ktauday, logn) ) THEN + IF ( IS_CASA_TIME("write", ts_start, logn) ) THEN ctime = ctime +1 @@ -935,8 +930,7 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) CALL master_receive (ocomm, oktau, casa_ts) IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND. & - ( IS_CASA_TIME("dwrit", yyyy, oktau, kstart, & - koffset, kend, ktauday, logn) ) ) THEN + ( IS_CASA_TIME("dwrit", ts_start, logn) THEN CALL master_receive ( ocomm, oktau, casa_dump_ts ) ENDIF diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index dd9567670..79fd48c86 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -224,6 +224,10 @@ SUBROUTINE mpidrv_worker (comm) REAL,ALLOCATABLE, SAVE :: c1(:,:) REAL,ALLOCATABLE, SAVE :: rhoch(:,:) REAL,ALLOCATABLE, SAVE :: xk(:,:) + + type(datetime) :: ts_start, ts_end + type(timedelta) :: dt + ! END header ! Maciej: make sure the variable does not go out of scope @@ -264,16 +268,37 @@ SUBROUTINE mpidrv_worker (comm) ! bal, logn, vegparmnew, casabiome, casapool, & ! casaflux, casamet, casabal, phen, C%EMSOIL, & ! C%TFRZ ) + ! Set the calendar- some forcing cases have special rules + select case (trim(cable_user%mettype)) + case ("plum") + if (plume%leapyears) then + call setcalendar("gregorian") + else + call setcalendar("noleaps") + endif + + case ("cru") + call setcalendar("noleaps") + + case default + if (leaps) then + call setcalendar("gregorian") + else + call setcalendar("noleaps") + endif + + end select SPINLOOP:DO YEAR: DO YYYY= CABLE_USER%YearStart, CABLE_USER%YearEnd CurYear = YYYY - IF ( leaps .AND. IS_LEAPYEAR( YYYY ) ) THEN - LOY = 366 - ELSE - LOY = 365 - ENDIF + ts_start = datetime(year=YYYY, month=1, day=1) + + LOY = daysInYear(YYYY) + + ! Set kend for all cases bar princeton and gswp + kend = ktauday * LOY ! MPI: receive from master ending time fields CALL MPI_Bcast (kend, 1, MPI_INTEGER, 0, comm, ierr) @@ -285,6 +310,7 @@ SUBROUTINE mpidrv_worker (comm) ! MPI: bcast to workers so that they don't need to open the met ! file themselves CALL MPI_Bcast (dels, 1, MPI_REAL, 0, comm, ierr) + dt = timedelta(second=int(dels)) ! MPI: need to know extents before creating datatypes CALL find_extents @@ -442,6 +468,7 @@ SUBROUTINE mpidrv_worker (comm) KTAULOOP:DO ktau=kstart, kend ! increment total timstep counter + ts_end = ts_start + dt ktau_tot = ktau_tot + 1 WRITE(logn,*) 'ktau -',ktau_tot @@ -459,7 +486,7 @@ SUBROUTINE mpidrv_worker (comm) !$ nyear =INT((kend-kstart+1)/(365*ktauday)) ! some things (e.g. CASA-CNP) only need to be done once per day - idoy =INT( MOD((REAL(ktau+koffset)/REAL(ktauday)),REAL(LOY))) + idoy = ts_start%getday() IF ( idoy .EQ. 0 ) idoy = LOY ! needed for CASA-CNP @@ -483,8 +510,7 @@ SUBROUTINE mpidrv_worker (comm) CALL MPI_Recv (MPI_BOTTOM, 1, inp_t, 0, ktau_gl, icomm, stat, ierr) ! MPI: receive casa_dump_data for this step from the master - ELSEIF ( IS_CASA_TIME("dread", yyyy, ktau, kstart, koffset, & - kend, ktauday, logn) ) THEN + ELSEIF ( IS_CASA_TIME("dread", curr_time, logn) ) THEN CALL MPI_Recv (MPI_BOTTOM, 1, casa_dump_t, 0, ktau_gl, icomm, stat, ierr) END IF @@ -533,17 +559,9 @@ SUBROUTINE mpidrv_worker (comm) ! ENDIF - IF ( IS_CASA_TIME("write", yyyy, ktau, kstart, & - koffset, kend, ktauday, logn) ) THEN - ! write(logn,*) 'IN IS_CASA', casapool%cplant(:,1) - ! CALL MPI_Send (MPI_BOTTOM,1, casa_t,0,ktau_gl,ocomm,ierr) - ENDIF - - ! MPI: send the results back to the master IF( ((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND. & - IS_CASA_TIME("dwrit", yyyy, ktau, kstart, & - koffset, kend, ktauday, logn)) & + IS_CASA_TIME("dwrit", ts_start, logn)) & CALL MPI_Send (MPI_BOTTOM, 1, casa_dump_t, 0, ktau_gl, ocomm, ierr) ENDIF @@ -568,6 +586,9 @@ SUBROUTINE mpidrv_worker (comm) CALL1 = .FALSE. + ! Update the time + ts_start = ts_end + END DO KTAULOOP ! END Do loop over timestep ktau ! ELSE @@ -599,20 +620,6 @@ SUBROUTINE mpidrv_worker (comm) ENDIF - - - - - - - - - - - - - - IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)).AND. & CABLE_USER%CALL_POP) THEN @@ -620,10 +627,8 @@ SUBROUTINE mpidrv_worker (comm) ENDIF - END DO YEAR - IF (spincasa .OR. casaonly) THEN EXIT ENDIF diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 0cbd8ab16..f66c861dd 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -796,6 +796,9 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) CALL1 = .FALSE. + ! Set time of new timestep + ts_start = ts_end + END DO ! END Do loop over timestep ktau CALL1 = .FALSE. From 946a0ccc4e904159d648faf7d42eca650b5a76c6 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 14:56:53 +1100 Subject: [PATCH 21/28] add missing module --- src/offline/cable_mpiworker.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index 79fd48c86..84d928a7f 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -159,6 +159,8 @@ SUBROUTINE mpidrv_worker (comm) USE POP_Constants, ONLY: HEIGHT_BINS, NCOHORT_MAX USE cbl_soil_snow_init_special_module + + use datetime_module IMPLICIT NONE ! MPI: From 4ad52602cdbc802d7c935ea5e661f45bd3350081 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 15:39:41 +1100 Subject: [PATCH 22/28] Move calender selection to the offline driver --- src/offline/cable_mpimaster.F90 | 22 ---------------------- src/offline/cable_mpiworker.F90 | 20 -------------------- src/offline/cable_offline_driver.F90 | 8 ++++++++ src/offline/cable_serial.F90 | 21 --------------------- 4 files changed, 8 insertions(+), 63 deletions(-) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index bc6b2c467..ac8bb97a7 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -334,28 +334,6 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) type(datetime) :: ts_start, ts_end type(timedelta) :: dt - ! Set the calendar- some forcing cases have special rules - select case (trim(cable_user%mettype)) - case ("plum") - if (plume%leapyears) then - call setcalendar("gregorian") - else - call setcalendar("noleaps") - endif - - case ("cru") - call setcalendar("noleaps") - - case default - if (leaps) then - call setcalendar("gregorian") - else - call setcalendar("noleaps") - endif - - end select - - ! END header ! outer loop - spinup loop no. ktau_tot : diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index 84d928a7f..871b508e2 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -270,26 +270,6 @@ SUBROUTINE mpidrv_worker (comm) ! bal, logn, vegparmnew, casabiome, casapool, & ! casaflux, casamet, casabal, phen, C%EMSOIL, & ! C%TFRZ ) - ! Set the calendar- some forcing cases have special rules - select case (trim(cable_user%mettype)) - case ("plum") - if (plume%leapyears) then - call setcalendar("gregorian") - else - call setcalendar("noleaps") - endif - - case ("cru") - call setcalendar("noleaps") - - case default - if (leaps) then - call setcalendar("gregorian") - else - call setcalendar("noleaps") - endif - - end select SPINLOOP:DO YEAR: DO YYYY= CABLE_USER%YearStart, CABLE_USER%YearEnd diff --git a/src/offline/cable_offline_driver.F90 b/src/offline/cable_offline_driver.F90 index d63e4fe13..746f3ae6a 100644 --- a/src/offline/cable_offline_driver.F90 +++ b/src/offline/cable_offline_driver.F90 @@ -18,6 +18,7 @@ PROGRAM cable_offline_driver USE CABLE_PLUME_MIP, ONLY : PLUME_MIP_TYPE USE CABLE_CRU, ONLY : CRU_TYPE USE CABLE_site, ONLY : site_TYPE + use datetime_module, only: set_calendar IMPLICIT NONE @@ -58,6 +59,13 @@ PROGRAM cable_offline_driver STOP END SELECT + ! Set the calendar + if (leaps) then + call set_calendar("gregorian") + else + call set_calender("noleaps") + end + IF (mpi_grp%size == 1) THEN CALL serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ELSE diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 6412fa890..4a8a5dc44 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -275,27 +275,6 @@ SUBROUTINE serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) ktauday = nint(d2s / dels) - ! Set the calendar- some forcing cases have special rules - select case (trim(cable_user%mettype)) - case ("plum") - if (plume%leapyears) then - call setcalendar("gregorian") - else - call setcalendar("noleaps") - endif - - case ("cru") - call setcalendar("noleaps") - - case default - if (leaps) then - call setcalendar("gregorian") - else - call setcalendar("noleaps") - endif - - end select - ! END header ! INISTUFF From 6e2cd1eb2cb7e969409a3abafc46174784741549 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 15:43:05 +1100 Subject: [PATCH 23/28] Correct setcalendar name --- src/offline/cable_offline_driver.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/offline/cable_offline_driver.F90 b/src/offline/cable_offline_driver.F90 index 746f3ae6a..36133eacc 100644 --- a/src/offline/cable_offline_driver.F90 +++ b/src/offline/cable_offline_driver.F90 @@ -18,7 +18,7 @@ PROGRAM cable_offline_driver USE CABLE_PLUME_MIP, ONLY : PLUME_MIP_TYPE USE CABLE_CRU, ONLY : CRU_TYPE USE CABLE_site, ONLY : site_TYPE - use datetime_module, only: set_calendar + use datetime_module, only: setcalendar IMPLICIT NONE @@ -61,9 +61,9 @@ PROGRAM cable_offline_driver ! Set the calendar if (leaps) then - call set_calendar("gregorian") + call setcalendar("gregorian") else - call set_calender("noleaps") + call setcalender("noleaps") end IF (mpi_grp%size == 1) THEN From 79e5cff61a07201f187ab74f2c3c1925fd226263 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 15:45:31 +1100 Subject: [PATCH 24/28] Import leaps to the main driver --- src/offline/cable_offline_driver.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/offline/cable_offline_driver.F90 b/src/offline/cable_offline_driver.F90 index 36133eacc..5b281e5c5 100644 --- a/src/offline/cable_offline_driver.F90 +++ b/src/offline/cable_offline_driver.F90 @@ -18,6 +18,7 @@ PROGRAM cable_offline_driver USE CABLE_PLUME_MIP, ONLY : PLUME_MIP_TYPE USE CABLE_CRU, ONLY : CRU_TYPE USE CABLE_site, ONLY : site_TYPE + use cable_io_vars_module, only: leaps use datetime_module, only: setcalendar IMPLICIT NONE From 3c321a18465cf15431c76a17e3126f3d13b6b856 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 15:46:14 +1100 Subject: [PATCH 25/28] Fix type --- src/offline/cable_offline_driver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/offline/cable_offline_driver.F90 b/src/offline/cable_offline_driver.F90 index 5b281e5c5..523fdd8f7 100644 --- a/src/offline/cable_offline_driver.F90 +++ b/src/offline/cable_offline_driver.F90 @@ -65,7 +65,7 @@ PROGRAM cable_offline_driver call setcalendar("gregorian") else call setcalender("noleaps") - end + endif IF (mpi_grp%size == 1) THEN CALL serialdrv(NRRRR, dels, koffset, kend, GSWP_MID, PLUME, CRU, site) From c1cd117b664378c3025710340ee950d8dea94fcd Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 15:47:10 +1100 Subject: [PATCH 26/28] Fix typo --- src/offline/cable_offline_driver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/offline/cable_offline_driver.F90 b/src/offline/cable_offline_driver.F90 index 523fdd8f7..967c353b7 100644 --- a/src/offline/cable_offline_driver.F90 +++ b/src/offline/cable_offline_driver.F90 @@ -64,7 +64,7 @@ PROGRAM cable_offline_driver if (leaps) then call setcalendar("gregorian") else - call setcalender("noleaps") + call setcalendar("noleaps") endif IF (mpi_grp%size == 1) THEN From af4d8470d1a0f5b34733badfff62182ec990120e Mon Sep 17 00:00:00 2001 From: Whyborn Date: Fri, 9 Jan 2026 16:00:33 +1100 Subject: [PATCH 27/28] Final couple of fixes --- src/offline/cable_mpimaster.F90 | 4 ++-- src/offline/cable_mpiworker.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index ac8bb97a7..6e3a3a2f0 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -907,8 +907,8 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) CALL master_receive (ocomm, oktau, casa_ts) - IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND. & - ( IS_CASA_TIME("dwrit", ts_start, logn) THEN + IF (((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND. & + ( IS_CASA_TIME("dwrit", ts_start, logn))) THEN CALL master_receive ( ocomm, oktau, casa_dump_ts ) ENDIF diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90 index 871b508e2..eb757c914 100644 --- a/src/offline/cable_mpiworker.F90 +++ b/src/offline/cable_mpiworker.F90 @@ -292,7 +292,7 @@ SUBROUTINE mpidrv_worker (comm) ! MPI: bcast to workers so that they don't need to open the met ! file themselves CALL MPI_Bcast (dels, 1, MPI_REAL, 0, comm, ierr) - dt = timedelta(second=int(dels)) + dt = timedelta(seconds=int(dels)) ! MPI: need to know extents before creating datatypes CALL find_extents @@ -492,7 +492,7 @@ SUBROUTINE mpidrv_worker (comm) CALL MPI_Recv (MPI_BOTTOM, 1, inp_t, 0, ktau_gl, icomm, stat, ierr) ! MPI: receive casa_dump_data for this step from the master - ELSEIF ( IS_CASA_TIME("dread", curr_time, logn) ) THEN + ELSEIF ( IS_CASA_TIME("dread", ts_start, logn) ) THEN CALL MPI_Recv (MPI_BOTTOM, 1, casa_dump_t, 0, ktau_gl, icomm, stat, ierr) END IF From 94487d05e5bbce7052e376b6bdd47bfd307fafc4 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Tue, 13 Jan 2026 10:13:46 +1100 Subject: [PATCH 28/28] Update ts_start at end of timestep --- src/offline/cable_mpimaster.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 6e3a3a2f0..cd487d245 100644 --- a/src/offline/cable_mpimaster.F90 +++ b/src/offline/cable_mpimaster.F90 @@ -886,6 +886,8 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) CALL1 = .FALSE. + ts_start = ts_end + !WRITE(*,*) " ktauloop end ", ktau, CurYear