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/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 - diff --git a/src/offline/CASAONLY_LUC.F90 b/src/offline/CASAONLY_LUC.F90 index 62103fba6..704b69686 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,15 +24,14 @@ 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 - INTEGER, INTENT(IN) :: LALLOC + 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 @@ -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 @@ -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 :: 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 @@ -97,8 +96,6 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,veg,soil,casabiome,casapool, & Iw = POP%Iwood 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 @@ -131,7 +128,6 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,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) @@ -158,7 +154,7 @@ SUBROUTINE CASAONLY_LUC( dels,kstart,kend,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, & @@ -270,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("write", curr_time, logn) ) THEN ctime = ctime +1 CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, & diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90 index 591caba01..cd487d245 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,7 +331,8 @@ 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 ! END header @@ -340,16 +343,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 +457,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 +674,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 +703,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 +742,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 +760,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 +776,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 +812,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 @@ -913,6 +886,8 @@ SUBROUTINE mpidrv_master (comm, dels, koffset, kend, PLUME, CRU) CALL1 = .FALSE. + ts_start = ts_end + !WRITE(*,*) " ktauloop end ", ktau, CurYear @@ -934,9 +909,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", yyyy, oktau, kstart, & - koffset, kend, ktauday, 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 @@ -8265,7 +8239,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..eb757c914 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: @@ -224,6 +226,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 @@ -269,11 +275,12 @@ SUBROUTINE mpidrv_worker (comm) 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 +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(seconds=int(dels)) ! MPI: need to know extents before creating datatypes CALL find_extents @@ -442,6 +450,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 +468,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 +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", yyyy, ktau, kstart, koffset, & - kend, ktauday, 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 @@ -533,17 +541,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 +568,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 +602,6 @@ SUBROUTINE mpidrv_worker (comm) ENDIF - - - - - - - - - - - - - - IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)).AND. & CABLE_USER%CALL_POP) THEN @@ -620,10 +609,8 @@ SUBROUTINE mpidrv_worker (comm) ENDIF - END DO YEAR - IF (spincasa .OR. casaonly) THEN EXIT ENDIF @@ -7042,7 +7029,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 +7197,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 +7319,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/cable_offline_driver.F90 b/src/offline/cable_offline_driver.F90 index d63e4fe13..967c353b7 100644 --- a/src/offline/cable_offline_driver.F90 +++ b/src/offline/cable_offline_driver.F90 @@ -18,6 +18,8 @@ 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 @@ -58,6 +60,13 @@ PROGRAM cable_offline_driver STOP END SELECT + ! Set the calendar + if (leaps) then + call setcalendar("gregorian") + else + call setcalendar("noleaps") + endif + 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 366aa58de..4a8a5dc44 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 @@ -268,6 +270,11 @@ 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 + + ktauday = nint(d2s / dels) + ! END header ! INISTUFF @@ -280,19 +287,24 @@ 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) + + 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 +319,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 (isleapyear(YYYY).AND.kend.EQ.2920) THEN STOP 'LEAP YEAR INCOMPATIBILITY WITH INPUT MET !' ENDIF IF ( NRRRR .GT. 1 ) THEN @@ -332,23 +344,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,8 +374,8 @@ 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) + ! Set the timestep here, after any met specific modifications of 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 @@ -446,7 +447,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. @@ -478,6 +479,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: ', isoformat(ts_start) WRITE(logn,*) 'Progress -',REAL(ktau)/REAL(kend)*100.0 ! increment total timstep counter @@ -486,7 +489,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 @@ -547,14 +551,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 @@ -562,7 +558,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 (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) ENDIF @@ -590,12 +587,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("dread", yyyy, ktau, kstart, & - koffset, kend, ktauday, logn) ) THEN ! CLN READ FROM FILE INSTEAD ! + 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) ) @@ -615,7 +607,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 @@ -625,8 +618,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, & @@ -649,24 +643,27 @@ 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("write", ts_start, 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)) @@ -767,7 +764,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) @@ -777,6 +775,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. 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..b29a7642d 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, & @@ -18,7 +18,6 @@ SUBROUTINE biogeochem(ktau,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 diff --git a/src/shared/casa_ncdf.F90 b/src/shared/casa_ncdf.F90 index 593f922b9..e9a80e674 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(iotype, curr_time, 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) :: curr_time + 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(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" @@ -472,8 +444,6 @@ FUNCTION IS_CASA_TIME(iotype, yyyy, ktau, kstart, koffset, kend, ktauday, logn) END FUNCTION IS_CASA_TIME - - END MODULE casa_ncdf_module diff --git a/src/util/datetime_module.f90 b/src/util/datetime_module.f90 new file mode 100644 index 000000000..8990ea5d1 --- /dev/null +++ b/src/util/datetime_module.f90 @@ -0,0 +1,1496 @@ +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 :: setcalendar + + 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)), private :: 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 setcalendar(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 setcalendar + + 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 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 + ! `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 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 + + isNewMonth = (d%getDay() == 1 .and. d%getHour() == 0 .and.& + d%getMinute() == 0 .and. d%getSecond() == 0 .and.& + d%getMillisecond() == 0) + end function isNewMonth + + + 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 + ! `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