diff --git a/.gitmodules b/.gitmodules index 6bb16a2c..a9dff09c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,7 +20,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/ESCOMP/atmospheric_physics - fxtag = 87e76a6ce90767d1d32259ff3f11dec9c9a03b03 + fxtag = 3f5435e8a9a53eb00f14a6b14975f92580b525d4 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "rrtmgp-data"] diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index 93a49504..9512c0f9 100644 --- a/cime_config/testdefs/testlist_cam.xml +++ b/cime_config/testdefs/testlist_cam.xml @@ -188,6 +188,17 @@ + + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq_trcdata_bam_derecho/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq_trcdata_bam_derecho/shell_commands new file mode 100644 index 00000000..76f3fe15 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq_trcdata_bam_derecho/shell_commands @@ -0,0 +1,2 @@ +./xmlchange CAM_CONFIG_OPTS="--dyn none --physics-suites tracer_data_test" +./xmlchange RUN_STARTDATE="1979-01-01" diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq_trcdata_bam_derecho/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq_trcdata_bam_derecho/user_nl_cam new file mode 100644 index 00000000..63698ed5 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq_trcdata_bam_derecho/user_nl_cam @@ -0,0 +1,53 @@ +! these snapshots were constructed to specifically test whether the data read into the +! physics buffer by tracer_data in CAM will match the data read by CAM-SIMA bit-for-bit. +! this was done by +! (1) taking a FHIST_C4 snapshot (which has prescribed_ozone and aero) +! in this case a ne3pg3_fhistc4_gw_drag_cam4_oro snapshot, +! (2) using nco to zero out the prescribed ozone and aero fields from the snapshot data, +! to create the "before" snapshot. +! (3) the "after" snapshot is the original snapshot file (which has non-zero values for prescribed +! ozone and aerosols) +! the test will pass iff. the tracer_data utility in CAM-SIMA and the prescribed_ozone/aerosols schemes +! can populate the prescribed fields bit-for-bit compared to the original snapshot file. + +ncdata = '/glade/campaign/cesm/community/amwg/sima_baselines/cam_sima_test_snapshots/cam_ne3pg3_fhistc4_tracer_data_0xprescribed-ozone-aero_from_cam4_oro_snapshot_derecho_gnu_after_c20251028.nc' +ncdata_check = '/glade/campaign/cesm/community/amwg/sima_baselines/cam_sima_test_snapshots/cam_ne3pg3_fhistc4_gw_drag_cam4_oro_snapshot_derecho_gnu_after_c20250826.nc' + +! tolerances for testing +ncdata_check_err = .true. +min_difference = 2e-15 + +! cam4 vertical levels (FHIST_C4 snapshots) +pver = 26 + +prescribed_aero_datapath = "/glade/campaign/cesm/cesmdata/inputdata/atm/cam/chem/trop_mozart_aero/aero" +prescribed_aero_file = "aero_1.9x2.5_L26_1850-2005_c091112.nc" +prescribed_aero_filelist = "aero_1.9x2.5_L26_list_c070514.txt" +prescribed_aero_model = "bulk" +prescribed_aero_specifier = 'sulf:SO4', 'bcar1:CB1', 'bcar2:CB2', 'ocar1:OC1', 'ocar2:OC2', 'sslt1:SSLT01', 'sslt2:SSLT02', 'sslt3:SSLT03', 'sslt4:SSLT04', 'dust1:DST01', 'dust2:DST02', 'dust3:DST03', 'dust4:DST04' +prescribed_aero_type = "INTERP_MISSING_MONTHS" +prescribed_aero_cycle_yr = 0 +prescribed_aero_fixed_tod = 0 +prescribed_aero_fixed_ymd = 0 + +aerodep_flx_datapath = '/glade/campaign/cesm/cesmdata/inputdata/atm/cam/chem/trop_mozart_aero/aero' +aerodep_flx_file = 'aerosoldep_monthly_1849-2006_1.9x2.5_c090803.nc' +aerodep_flx_specifier = 'BCDEPWET', 'BCPHODRY', 'BCPHIDRY', 'OCDEPWET', 'OCPHODRY', 'OCPHIDRY', 'DSTX01DD', 'DSTX02DD', 'DSTX03DD', + 'DSTX04DD', 'DSTX01WD', 'DSTX02WD', 'DSTX03WD', 'DSTX04WD' +aerodep_flx_type = 'INTERP_MISSING_MONTHS' +aerodep_flx_cycle_yr = 0 + +prescribed_ozone_datapath = "/glade/campaign/cesm/cesmdata/inputdata/atm/cam/ozone" +prescribed_ozone_file = "ozone_1.9x2.5_L26_1850-2005_c091112.nc" +prescribed_ozone_filelist = "UNSET" +prescribed_ozone_name = "O3" +prescribed_ozone_type = "INTERP_MISSING_MONTHS" +prescribed_ozone_fixed_tod = 0 +prescribed_ozone_fixed_ymd = 0 +prescribed_ozone_cycle_yr = 0 + +hist_output_frequency;h1: 1*nsteps +hist_max_frames;h1: 1 +hist_add_inst_fields;h1:ozone +hist_add_inst_fields;h1:sulf_D,bcar1_D,bcar2_D,ocar1_D,ocar2_D,sslt1_D,sslt2_D,sslt3_D,sslt4_D,dust1_D,dust2_D,dust3_D,dust4_D +hist_precision;h1: REAL64 diff --git a/src/core_utils/string_core_utils.F90 b/src/core_utils/string_core_utils.F90 index 04c3cced..1970ac9f 100644 --- a/src/core_utils/string_core_utils.F90 +++ b/src/core_utils/string_core_utils.F90 @@ -9,6 +9,9 @@ module string_core_utils public :: split ! Parse a string into tokens, one at a time public :: stringify ! Convert one or more values of any intrinsic data types to a character string for pretty printing public :: tokenize ! Parse a string into tokens + public :: increment_string ! Increment a string whose ending characters are digits. + public :: last_non_digit ! Get position of last non-digit in the input string. + public :: get_last_significant_char ! Get position of last significant (non-blank, non-null) character in string. interface tokenize module procedure tokenize_into_first_last @@ -294,4 +297,112 @@ pure subroutine tokenize_into_tokens_separator(string, set, tokens, separator) end if end subroutine tokenize_into_tokens_separator + ! Increment a string whose ending characters are digits. + ! The incremented integer must be in the range [0 - (10**n)-1] + ! where n is the number of trailing digits. + ! Return values: + ! + ! 0 success + ! -1 error: no trailing digits in string + ! -2 error: incremented integer is out of range + integer function increment_string(s, inc) + integer, intent(in) :: inc ! value to increment string (may be negative) + character(len=*), intent(inout) :: s ! string with trailing digits + + integer :: & + i, & ! index + lstr, & ! number of significant characters in string + lnd, & ! position of last non-digit + ndigit, & ! number of trailing digits + ival, & ! integer value of trailing digits + pow, & ! power of ten + digit ! integer value of a single digit + + lstr = get_last_significant_char(s) + lnd = last_non_digit(s) + ndigit = lstr - lnd + + if(ndigit == 0) then + increment_string = -1 + return + end if + + ! Calculate integer corresponding to trailing digits. + ival = 0 + pow = 0 + do i = lstr,lnd+1,-1 + digit = ICHAR(s(i:i)) - ICHAR('0') + ival = ival + digit * 10**pow + pow = pow + 1 + end do + + ! Increment the integer + ival = ival + inc + if( ival < 0 .or. ival > 10**ndigit-1 ) then + increment_string = -2 + return + end if + + ! Overwrite trailing digits + pow = ndigit + do i = lnd+1,lstr + digit = MOD( ival,10**pow ) / 10**(pow-1) + s(i:i) = CHAR( ICHAR('0') + digit ) + pow = pow - 1 + end do + + increment_string = 0 + + end function increment_string + + ! Get position of last non-digit in the input string. + ! Return values: + ! > 0 => position of last non-digit + ! = 0 => token is all digits (or empty) + integer pure function last_non_digit(s) + character(len=*), intent(in) :: s + integer :: n, nn, digit + + n = get_last_significant_char(s) + if(n == 0) then ! empty string + last_non_digit = 0 + return + end if + + do nn = n,1,-1 + digit = ICHAR(s(nn:nn)) - ICHAR('0') + if( digit < 0 .or. digit > 9 ) then + last_non_digit = nn + return + end if + end do + + last_non_digit = 0 ! all characters are digits + + end function last_non_digit + + ! Get position of last significant character in string. + ! Here significant means non-blank or non-null. + ! Return values: + ! > 0 => position of last significant character + ! = 0 => no significant characters in string + integer pure function get_last_significant_char(cs) + character(len=*), intent(in) :: cs ! Input character string + integer :: l, n + + l = LEN(cs) + if( l == 0 ) then + get_last_significant_char = 0 + return + end if + + do n = l,1,-1 + if( cs(n:n) /= ' ' .and. cs(n:n) /= CHAR(0) ) then + exit + end if + end do + get_last_significant_char = n + + end function get_last_significant_char + end module string_core_utils diff --git a/src/data/registry.xml b/src/data/registry.xml index 9ac7ca92..83ba849e 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -1627,6 +1627,86 @@ pbuf_clubbtop + + + horizontal_dimension vertical_layer_dimension + sulf pbuf_sulf cnst_sulf + + + horizontal_dimension vertical_layer_dimension + bcar1 pbuf_bcar1 cnst_bcar1 + + + horizontal_dimension vertical_layer_dimension + bcar2 pbuf_bcar2 cnst_bcar2 + + + horizontal_dimension vertical_layer_dimension + ocar1 pbuf_ocar1 cnst_ocar1 + + + horizontal_dimension vertical_layer_dimension + ocar2 pbuf_ocar2 cnst_ocar2 + + + horizontal_dimension vertical_layer_dimension + sslt1 pbuf_sslt1 cnst_sslt1 + + + horizontal_dimension vertical_layer_dimension + sslt2 pbuf_sslt2 cnst_sslt2 + + + horizontal_dimension vertical_layer_dimension + sslt3 pbuf_sslt3 cnst_sslt3 + + + horizontal_dimension vertical_layer_dimension + sslt4 pbuf_sslt4 cnst_sslt4 + + + horizontal_dimension vertical_layer_dimension + dust1 pbuf_dust1 cnst_dust1 + + + horizontal_dimension vertical_layer_dimension + dust2 pbuf_dust2 cnst_dust2 + + + horizontal_dimension vertical_layer_dimension + dust3 pbuf_dust3 cnst_dust3 + + + horizontal_dimension vertical_layer_dimension + dust4 pbuf_dust4 cnst_dust4 + + max_diff(1)) then - max_diff(1) = diff - max_diff_col = col - end if - if (diff > min_difference) then - diff_count = diff_count + 1 + ! Calculate actual diffs for non-NaN values: + if (abs(current_value(col)) < min_relative_value) then + !Calculate absolute difference: + diff = abs(current_value(col) - buffer(col)) + else + !Calculate relative difference: + diff = abs(current_value(col) - buffer(col)) / & + abs(current_value(col)) + end if + ! Only update max diff if greater and not already nan: + if (diff > max_diff(1) .and. .not. has_nan) then + max_diff(1) = diff + max_diff_col = col + end if + if (diff > min_difference) then + diff_count = diff_count + 1 + end if end if end do + + !Add NaN count to total difference count: + diff_count = diff_count + nan_count + !Gather results across all nodes to get global values + call mpi_reduce(nan_count, nan_count_gl, 1, mpi_integer, & + mpi_sum, masterprocid, mpicom, ierr) call mpi_reduce(diff_count, diff_count_gl, 1, mpi_integer, & - mpi_sum, masterprocid, mpicom, ierr) + mpi_sum, masterprocid, mpicom, ierr) call mpi_allreduce(max_diff, max_diff_gl, 1, & MPI_2DOUBLE_PRECISION, & @@ -737,7 +765,8 @@ subroutine check_field_2d(file, var_names, timestep, current_value, & !Print difference stats to log file if (masterproc) then if (diff_count_gl > 0) then - call write_check_field_entry(stdname, diff_count_gl, & + call write_check_field_entry(stdname, diff_count_gl, & + nan_count_gl, & max_diff_gl(1), & int(max_diff_gl(2)), & max_diff_gl_col, is_first) @@ -764,6 +793,7 @@ subroutine check_field_3d(file, var_names, vcoord_name, timestep, & use mpi, only: mpi_maxloc, mpi_sum, mpi_status_size use mpi, only: mpi_2double_precision, mpi_integer use vert_coord, only: pver, pverp + use shr_infnan_mod, only: shr_infnan_isnan !Max possible length of variable name in file: use phys_vars_init_check, only: std_name_len @@ -800,6 +830,9 @@ subroutine check_field_3d(file, var_names, vcoord_name, timestep, & integer :: max_diff_gl_col integer :: max_diff_gl_lev integer :: diff_count_gl + integer :: nan_count ! Count of NaNs found + integer :: nan_count_gl ! Global count of NaNs + logical :: has_nan ! Flag indicating NaN was found !Initialize output variables ierr = 0 @@ -832,29 +865,52 @@ subroutine check_field_3d(file, var_names, vcoord_name, timestep, & timelevel=timestep, dim3name=trim(vcoord_name), & dim3_bnds=[1, num_levs], log_output=.false.) if (var_found) then + nan_count = 0 + has_nan = .false. + do lev = 1, num_levs do col = 1, size(buffer(:,lev)) - if (abs(current_value(col, lev)) < min_relative_value) then - !Calculate absolute difference: - diff = abs(current_value(col, lev) - buffer(col, lev)) + ! First, check if there are NaNs anywhere in the state + if (shr_infnan_isnan(current_value(col, lev))) then + nan_count = nan_count + 1 + + if (.not. has_nan) then ! First NaN found for this variable + has_nan = .true. + + ! Set max diff to NaN (if not already) to signal NaN was found + max_diff(1) = current_value(col,lev) - buffer(col,lev) ! = nan + max_diff_col = col + max_diff_lev = lev + end if else - !Calculate relative difference: - diff = abs(current_value(col, lev) - buffer(col, lev)) / & - abs(current_value(col, lev)) - end if - if (diff > max_diff(1)) then - max_diff(1) = diff - max_diff_col = col - max_diff_lev = lev - end if - !Determine if diff is large enough to be considered a "hit" - if (diff > min_difference) then - diff_count = diff_count + 1 + ! Calculate actual diffs for non-NaN values: + if (abs(current_value(col, lev)) < min_relative_value) then + !Calculate absolute difference: + diff = abs(current_value(col, lev) - buffer(col, lev)) + else + !Calculate relative difference: + diff = abs(current_value(col, lev) - buffer(col, lev)) / & + abs(current_value(col, lev)) + end if + if (diff > max_diff(1)) then + max_diff(1) = diff + max_diff_col = col + max_diff_lev = lev + end if + !Determine if diff is large enough to be considered a "hit" + if (diff > min_difference) then + diff_count = diff_count + 1 + end if end if end do end do - !Make relevant MPI calls to get global values: + !Add NaN count to total difference count: + diff_count = diff_count + nan_count + + !Gather results across all nodes to get global values + call mpi_reduce(nan_count, nan_count_gl, 1, mpi_integer, & + mpi_sum, masterprocid, mpicom, ierr) call mpi_reduce(diff_count, diff_count_gl, 1, mpi_integer, & mpi_sum, masterprocid, mpicom, ierr) @@ -889,11 +945,12 @@ subroutine check_field_3d(file, var_names, vcoord_name, timestep, & !Print difference stats to log file if (masterproc) then if (diff_count_gl > 0) then - call write_check_field_entry(stdname, diff_count_gl, & + call write_check_field_entry(stdname, diff_count_gl, & + nan_count_gl, & max_diff_gl(1), & int(max_diff_gl(2)), & - max_diff_gl_col, & - is_first, & + max_diff_gl_col, & + is_first, & max_diff_lev=max_diff_gl_lev) is_first = .false. diff_found = .true. @@ -958,6 +1015,9 @@ subroutine check_field_4d(file, var_names, vcoord_name, timestep, & integer :: max_diff_gl_lev integer :: max_diff_gl_extra_dim integer :: diff_count_gl + integer :: nan_count + integer :: nan_count_gl + logical :: has_nan !Initialize output variables ierr = 0 @@ -972,6 +1032,8 @@ subroutine check_field_4d(file, var_names, vcoord_name, timestep, & max_diff(1) = 0._kind_phys max_diff(2) = real(iam, kind_phys) !MPI rank for this task diff_found = .false. + nan_count = 0 + has_nan = .false. call cam_pio_find_var(file, var_names, found_name, vardesc, var_found) if (.not. var_found) then @@ -997,28 +1059,53 @@ subroutine check_field_4d(file, var_names, vcoord_name, timestep, & do extra_dim = 1, size(buffer, 3) do lev = 1, num_levs do col = 1, size(buffer(:,lev,extra_dim)) - if (abs(current_value(col, lev, extra_dim)) < min_relative_value) then - !Calculate absolute difference: - diff = abs(current_value(col, lev, extra_dim) - buffer(col, lev, extra_dim)) + ! Check for NaNs first + if (current_value(col, lev, extra_dim) /= current_value(col, lev, extra_dim) .or. & + buffer(col, lev, extra_dim) /= buffer(col, lev, extra_dim)) then + nan_count = nan_count + 1 + if (.not. has_nan) then + has_nan = .true. + + ! Force max_diff to NaN to signal NaN found + max_diff(1) = current_value(col, lev, extra_dim) - & + buffer(col, lev, extra_dim) + max_diff_col = col + max_diff_lev = lev + max_diff_extra_dim = extra_dim + end if else - !Calculate relative difference: - diff = abs(current_value(col, lev, extra_dim) - buffer(col, lev, extra_dim)) / & - abs(current_value(col, lev, extra_dim)) - end if - if (diff > max_diff(1)) then - max_diff(1) = diff - max_diff_col = col - max_diff_lev = lev - max_diff_extra_dim = extra_dim - end if - !Determine if diff is large enough to be considered a "hit" - if (diff > min_difference) then - diff_count = diff_count + 1 + if (abs(current_value(col, lev, extra_dim)) < min_relative_value) then + ! Absolute difference + diff = abs(current_value(col, lev, extra_dim) - & + buffer(col, lev, extra_dim)) + else + ! Relative difference + diff = abs(current_value(col, lev, extra_dim) - & + buffer(col, lev, extra_dim)) / & + abs(current_value(col, lev, extra_dim)) + end if + + if (diff > max_diff(1)) then + max_diff(1) = diff + max_diff_col = col + max_diff_lev = lev + max_diff_extra_dim = extra_dim + end if + + if (diff > min_difference) then + diff_count = diff_count + 1 + end if end if end do end do end do + !Add NaN count to total difference count: + diff_count = diff_count + nan_count + + call mpi_reduce(nan_count, nan_count_gl, 1, mpi_integer, & + mpi_sum, masterprocid, mpicom, ierr) + !Make relevant MPI calls to get global values: call mpi_reduce(diff_count, diff_count_gl, 1, mpi_integer, & mpi_sum, masterprocid, mpicom, ierr) @@ -1061,6 +1148,7 @@ subroutine check_field_4d(file, var_names, vcoord_name, timestep, & if (masterproc) then if (diff_count_gl > 0) then call write_check_field_entry(stdname, diff_count_gl, & + nan_count_gl, & max_diff_gl(1), & int(max_diff_gl(2)), & max_diff_gl_col, & @@ -1077,7 +1165,8 @@ subroutine check_field_4d(file, var_names, vcoord_name, timestep, & end subroutine check_field_4d - subroutine write_check_field_entry(stdname, diff_count, & + subroutine write_check_field_entry(stdname, & + diff_count, nan_count, & max_diff, max_diff_rank, & max_diff_col, is_first, & max_diff_lev, max_diff_extra_dim) @@ -1089,6 +1178,7 @@ subroutine write_check_field_entry(stdname, diff_count, & !Dummy variables: character(len=*), intent(in) :: stdname integer, intent(in) :: diff_count + integer, intent(in) :: nan_count real(kind_phys), intent(in) :: max_diff integer, intent(in) :: max_diff_rank !MPI rank max diff occurred on integer, intent(in) :: max_diff_col !max diff column (1st) dimension value @@ -1119,7 +1209,7 @@ subroutine write_check_field_entry(stdname, diff_count, & slen = 0 end if - !Write out difference and index valuesa: + !Write out difference and index values: if (present(max_diff_lev)) then if(present(max_diff_extra_dim)) then write(index_str, '(a,i0,a,i0,a,i0,a,i0,a)') "(",max_diff_rank,",",max_diff_col,",",max_diff_lev,",",max_diff_extra_dim,")" @@ -1132,6 +1222,12 @@ subroutine write_check_field_entry(stdname, diff_count, & write(fmt_str, '(a,i0,a)') "(1x,a,t",indent_level+1,",1x,i7,2x,e8.2,3x,a)" write(iulog, fmt_str) stdname(1:slen), diff_count, max_diff, index_str + !Print out warning if any NaN values found: + if (nan_count > 0) then + write(iulog, '(a,i0,a)') ' (!) ', nan_count, & + ' NaN values in variable!' + end if + end subroutine write_check_field_entry end module physics_data diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90 index aad15594..1496a889 100644 --- a/src/physics/utils/physics_grid.F90 +++ b/src/physics/utils/physics_grid.F90 @@ -59,7 +59,7 @@ module physics_grid ! hdim1_d and hdim2_d are dimensions of rectangular horizontal grid ! data structure, If 1D data structure, then hdim2_d == 1. integer, protected, public :: hdim1_d, hdim2_d - logical :: dycore_unstructured = .false. + logical, public :: dycore_unstructured = .false. ! Dycore name and properties character(len=8), protected, public :: dycore_name = '' @@ -131,7 +131,6 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, dycore_name_in, & real(r8), pointer :: area_d(:) real(r8) :: mem_hw_beg, mem_hw_end real(r8) :: mem_beg, mem_end - logical :: unstructured real(r8) :: temp ! For MPI integer :: ierr ! For error codes @@ -160,8 +159,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, dycore_name_in, & hdim1_d = hdim1_d_in hdim2_d = hdim2_d_in dycore_name = dycore_name_in - - unstructured = hdim2_d <= 1 + dycore_unstructured = hdim2_d <= 1 ! Calculate total number of physics columns: num_global_phys_cols = hdim1_d * hdim2_d @@ -225,7 +223,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, dycore_name_in, & ! First, create a map for the physics grid ! It's structure will depend on whether or not the physics grid is ! unstructured - if (unstructured) then + if (dycore_unstructured) then allocate(grid_map(3, columns_on_task), stat=ierr) call check_allocate(ierr, subname, 'grid_map(3, columns_on_task)', & file=__FILE__, line=__LINE__) @@ -257,7 +255,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, dycore_name_in, & end if grid_map(1, col_index) = col_index grid_map(2, col_index) = 0 ! No chunking in physics anymore - if (unstructured) then + if (dycore_unstructured) then grid_map(3, col_index) = phys_columns(col_index)%global_col_num else ! lon @@ -272,7 +270,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, dycore_name_in, & ! the physics grid ! However, these will be in the dynamics decomposition - if (unstructured) then + if (dycore_unstructured) then lon_coord => horiz_coord_create('lon', 'ncol', num_global_phys_cols, & 'longitude', 'degrees_east', 1, size(lonvals), lonvals, & map=grid_map(3,:)) @@ -314,7 +312,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, dycore_name_in, & end if call cam_grid_register('physgrid', phys_decomp, & lat_coord, lon_coord, grid_map, src_in=(/ 1, 0 /), & - unstruct=unstructured, block_indexed=.false.) + unstruct=dycore_unstructured, block_indexed=.false.) ! Copy required attributes from the dynamics array do index = 1, size(dyn_attributes) @@ -323,7 +321,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, dycore_name_in, & end do if ((.not. cam_grid_attr_exists('physgrid', 'area')) .and. & - unstructured) then + dycore_unstructured) then ! Physgrid always needs an area attribute. If we did not inherit one ! from the dycore (i.e., physics and dynamics are on different ! grids), create that attribute here (Note, a separate physics diff --git a/src/utils/horizontal_interpolate.F90 b/src/utils/horizontal_interpolate.F90 new file mode 100644 index 00000000..f3ae6bcf --- /dev/null +++ b/src/utils/horizontal_interpolate.F90 @@ -0,0 +1,352 @@ +! Horizontal interpolation from CAM chemistry utilities. +module horizontal_interpolate + + use shr_kind_mod, only: r8 => shr_kind_r8 + use physconst, only: pi + use cam_abortutils, only: endrun + + implicit none + private + + public :: xy_interp_init, xy_interp + +contains + + subroutine check_invariant(invariant_condition, message) + logical, intent(in) :: invariant_condition + character(len=*), intent(in) :: message + if (.not. (invariant_condition)) then + call endrun("Invariant check failed: "//message) + end if + end subroutine check_invariant + + real(r8) pure function normalize_lon_right(left, right) + ! Normalize the right side of a longitude interval such that + ! norm_right > left and (norm_right - left) is in (0, 360] + real(r8), intent(in) :: left, right + + normalize_lon_right = right + + do while (normalize_lon_right <= left) + normalize_lon_right = normalize_lon_right + 360.0_r8 + end do + do while (normalize_lon_right - 360.0_r8 > left) + normalize_lon_right = normalize_lon_right - 360.0_r8 + end do + end function normalize_lon_right + + real(r8) pure function lon_length(left, right) + ! Compute a longitude interval length, accouting for wrapping + ! around the globe. Input values are normalized so that the + ! return value is always in the range (0, 360]. + real(r8), intent(in) :: left, right + + lon_length = normalize_lon_right(left, right) - left + end function lon_length + + real(r8) pure function normalize_lon_value(lon) + ! Normalize a lon value to be in [0, 360) + real(r8), intent(in) :: lon + real(r8) :: norm_lon + + norm_lon = lon + do while (norm_lon < 0) + norm_lon = norm_lon + 360.0_r8 + end do + do while (norm_lon >= 360.0_r8) + norm_lon = norm_lon - 360.0_r8 + end do + + normalize_lon_value = norm_lon + end function normalize_lon_value + + ! Return the length of the overlap between the input and + ! simulation longitude ranges. Values are normalized before + ! calculation to ensure correct handling of ranges which wrap over + ! zero. + real(r8) pure function calculate_lon_overlap(input_left, input_right, sim_left, sim_right) + real(r8), intent(in) :: input_left, input_right, sim_left, sim_right + + real(r8) :: norm_input_left, norm_input_right, norm_sim_left, norm_sim_right + real(r8) :: overlap_left, overlap_right + + ! Normalize so norm_sim_left is in [0, 360) and norm_sim_left < norm_sim_right + norm_sim_left = normalize_lon_value(sim_left) + norm_sim_right = normalize_lon_right(norm_sim_left, sim_right) + + ! Normalize the input values to ensure that norm_input_left < norm_sim_left + norm_input_left = normalize_lon_value(input_left) - 360.0_r8 ! now in [-360, 0) + norm_input_right = normalize_lon_right(norm_input_left, input_right) + + ! if norm_input is strictly to the left of norm_sim, slide up by 360 + do while (norm_input_right <= norm_sim_left) + norm_input_left = norm_input_left + 360.0_r8 + norm_input_right = norm_input_right + 360.0_r8 + end do + + ! Compute overlap + overlap_left = merge(norm_input_left, norm_sim_left, norm_input_left > norm_sim_left) + overlap_right = merge(norm_input_right, norm_sim_right, norm_input_right < norm_sim_right) + if (overlap_left < overlap_right) then + calculate_lon_overlap = overlap_right - overlap_left + else + calculate_lon_overlap = 0 + end if + end function calculate_lon_overlap + + real(r8) function lon_weight(input_left, input_right, sim_left, sim_right, use_flight_distance) + ! Compute how much input data in the (input_left, input_right) + ! longitude band contributes to simulation data in the (sim_left, + ! sim_right) band. + ! + ! use_flight_distance indicates how the input data should be + ! interpreted. If true, input data is assumed to be a total value + ! for the input grid cell. If false, input data is assumed to be + ! a mixing ratio. + real(r8), intent(in) :: input_left, input_right, sim_left, sim_right + logical, intent(in) :: use_flight_distance !.true. = flight distance, .false. = all mixing ratios + + real(r8) :: overlap_len + + ! Sanity check that inputs aren't too large or too small. For + ! really huge floating point values, adding / subtracting 360 will + ! encounter roundoff errors, and could cause while loops to go + ! forever. Values much outside the [-360, 360] range are probably + ! an error anyway, so abort if we encounter anything suspicious here. + call check_invariant(((-2000.0_r8 < input_left) .and. (input_left < 2000.0_r8)), "input_left is out of range") + call check_invariant(((-2000.0_r8 < input_right) .and. (input_right < 2000.0_r8)), "input_right is out of range") + call check_invariant(((-2000.0_r8 < sim_left) .and. (sim_left < 2000.0_r8)), "sim_left is out of range") + call check_invariant(((-2000.0_r8 < sim_right) .and. (sim_right < 2000.0_r8)), "sim_right is out of range") + + overlap_len = calculate_lon_overlap(input_left, input_right, sim_left, sim_right) + + if (overlap_len == 0) then + ! No overlap; weight is zero + lon_weight = 0 + + else if (use_flight_distance) then + ! Data values are total within the grid cell. Hence, the + ! weight is just the fraction of the area of the original grid + ! cell which overlaps the new cell. + lon_weight = overlap_len/lon_length(input_left, input_right) + else + ! Data values are mixing ratios. To compute how much this + ! input grid cell contributes to the mixing ratio in the + ! simulation grid cell, we multiply by the fraction of the sim + ! grid cell occupied by the overlap. Fortunately this is just + ! the ratio of the longitude ranges, since the grid cells are + ! rectangular in lat,lon coordinates. + lon_weight = overlap_len/lon_length(sim_left, sim_right) + end if + + ! Floating point precision can wind up with lon_weight slightly + ! out of the [0, 1] range, so fix up the values if that happens. + ! If things are farther off, then crash with an error. + call check_invariant((-0.0000001_r8 <= lon_weight) .and. (lon_weight <= 1.0000001_r8), "lon_weight must be in [0, 1]") + if (lon_weight < 0.0_r8) then + lon_weight = 0.0_r8 + end if + if (lon_weight > 1.0_r8) then + lon_weight = 1.0_r8 + end if + end function lon_weight + + real(r8) pure function lat_band_weight(low, high) + ! Retuns the unit-less weight of the latitude band on the globe + ! from [low, high]. The area of a latitudinal slice of height h + ! is just 2*pi*r*h. Since we only need to get ratios of areas, + ! this just computes the difference in the y-coordinates from low + ! to high. + real(r8), intent(in) :: low, high + + lat_band_weight = sin(high*pi/180.0_r8) - sin(low*pi/180.0_r8) + end function lat_band_weight + + real(r8) function normalize_lat(lat) + ! Truncate a latitude value to be within [-90, 90]. Values much + ! outside that range are probably an error in the calling code, so + ! exit if we see those. + real(r8), intent(in) :: lat + call check_invariant(((-91.0_r8 < lat) .and. (lat < 91.0_r8)), "Lat value is out of expected range.") + if (lat < -90.0_r8) then + normalize_lat = -90.0_r8 + else if (lat > 90.0_r8) then + normalize_lat = 90.0_r8 + else + normalize_lat = lat + end if + end function normalize_lat + + real(r8) function lat_weight(input_bot, input_top, sim_bot, sim_top, use_flight_distance) + ! Compute how much input data in the (input_left, input_right) + ! longitude band contributes to simulation data in the (sim_left, + ! sim_right) band. + ! + ! use_flight_distance indicates how the input data should be + ! interpreted. If true, input data is assumed to be a total value + ! for the input grid cell. If false, input data is assumed to be + ! a mixing ratio. + real(r8), intent(in) :: input_bot, input_top, sim_bot, sim_top + logical, intent(in) :: use_flight_distance !.true. = flight distance, .false. = all mixing ratios + real(r8) :: norm_input_bot, norm_input_top, norm_sim_bot, norm_sim_top + real(r8) :: overlap_bot, overlap_top + + ! Make sure that inputs are in [-90, 90] and that bot < top + norm_input_bot = normalize_lat(input_bot) + norm_input_top = normalize_lat(input_top) + norm_sim_bot = normalize_lat(sim_bot) + norm_sim_top = normalize_lat(sim_top) + call check_invariant(norm_input_bot < norm_input_top, "must have input_bot < input_top") + call check_invariant(norm_sim_bot < norm_sim_top, "must have sim_bot < sim_top") + + overlap_bot = merge(norm_input_bot, norm_sim_bot, norm_input_bot > norm_sim_bot) + overlap_top = merge(norm_input_top, norm_sim_top, norm_input_top < norm_sim_top) + + if ((norm_input_bot == norm_input_top) .or. & + (norm_sim_bot == norm_sim_top) .or. & + (overlap_top <= overlap_bot)) then + ! No overlap + lat_weight = 0 + else if (use_flight_distance) then + ! Input values are a total for the grid cell, so the weight is + ! just the fraction of the input cell which overlaps the + ! sim cell. + lat_weight = lat_band_weight(overlap_bot, overlap_top)/lat_band_weight(norm_input_bot, norm_input_top) + call check_invariant((0 <= lat_weight) .and. (lat_weight <= 1), "dist: lat_weight must be in [0, 1]") + else + ! Input values are a mixing ratio. The amount added in the + ! overlap is spread evenly over the sim cell. + lat_weight = lat_band_weight(overlap_bot, overlap_top)/lat_band_weight(norm_sim_bot, norm_sim_top) + call check_invariant((0 <= lat_weight) .and. (lat_weight <= 1), "mixrat: lat_weight must be in [0, 1]") + end if + end function lat_weight + + ! This program computes weighting functions to map a variable of (num_input_lons,num_input_lats) resolution to (num_sim_lons,num_sim_lats) resolution + ! weight_x(num_sim_lons,num_input_lons) is the weighting function for zonal interpolation + ! weight_y(num_sim_lats,num_input_lats) is the weighting function for meridional interpolation + ! + ! Author: Chih-Chieh (Jack) Chen -- May 2010 + ! Rob von Behren (jrvb@google.com) -- Oct 2024 + subroutine xy_interp_init(num_input_lons, num_input_lats, input_lon_radians, input_lat_radians, & + num_sim_lons, num_sim_lats, weight_x, weight_y, use_flight_distance) + integer, intent(in) :: num_input_lons, num_input_lats, num_sim_lons, num_sim_lats + logical, intent(in) :: use_flight_distance ! .true. = flight distance, .false. = all mixing ratios + real(r8), intent(in) :: input_lon_radians(num_input_lons), input_lat_radians(num_input_lats) + real(r8), intent(out) :: weight_x(num_sim_lons, num_input_lons), weight_y(num_sim_lats, num_input_lats) + + real(r8) :: input_lon(num_input_lons), input_lat(num_input_lats) + real(r8) :: sim_lon(num_sim_lons), sim_lat(num_sim_lats) + real(r8) :: input_lon_edge(num_input_lons + 1), sim_lon_edge(num_sim_lons + 1) + real(r8) :: input_lat_edge(num_input_lats + 1), sim_lat_edge(num_sim_lats + 1) + real(r8) :: x1_west, x1_east, x2_west, x2_east + real(r8) :: y1_south, y1_north, y2_south, y2_north + integer :: i1, j1, i2, j2, i + + weight_x(:, :) = 0.0_r8 + weight_y(:, :) = 0.0_r8 + + ! input_lon_radians & input_lat_radians are longitude & latitude on the source mesh in radians + ! convert input_lon, input_lat from radians to degrees + input_lon(:) = input_lon_radians(:)/pi*180.0_r8 + input_lat(:) = input_lat_radians(:)/pi*180.0_r8 + + ! Set up sim_lon, sim_lat (target mesh), in CAM convention. The + ! (lon,lat) pairs are the center points of the grid cells. + do i2 = 1, num_sim_lons + sim_lon(i2) = (float(i2) - 1.0_r8)*360.0_r8/float(num_sim_lons) + end do + do j2 = 1, num_sim_lats + sim_lat(j2) = -90.0_r8 + (float(j2) - 1.0_r8)*180.0_r8/(float(num_sim_lats) - 1.0_r8) + end do + ! make sure the highest value is exactly +90, since the + ! multiplication above could give a value that is slightly off due + ! to rounding. + sim_lat(num_sim_lats) = 90.0_r8 + + ! Calculate the grid cell edges from the midpoints above. We set + ! things up so it is easy to find the boundary for grid cell i, j as: + ! lat bounds: (lat_edge(i), lat_edge(i+1)) + ! lon bounds: (lon_edge(j), lon_edge(j+1)) + ! For latitudes, we ensure that the smallest and largest edges are + ! at -90 and 90, respectively. For longitudes, the values wrap so + ! lon_edge(1) == lon_edge(num+1) - 360 + + ! Input longitude edges + do i1 = 2, num_input_lons + input_lon_edge(i1) = (input_lon(i1 - 1) + input_lon(i1))/2.0_r8 + end do + input_lon_edge(1) = (input_lon(num_input_lons) - 360_r8 + input_lon(1))/2.0_r8 + input_lon_edge(num_input_lons + 1) = input_lon_edge(1) + 360_r8 + + do i2 = 2, num_sim_lons + sim_lon_edge(i2) = (sim_lon(i2 - 1) + sim_lon(i2))/2.0_r8 + end do + sim_lon_edge(1) = (sim_lon(num_sim_lons) - 360_r8 + sim_lon(1))/2.0_r8 + sim_lon_edge(num_sim_lons + 1) = sim_lon_edge(1) + 360_r8 + + ! set up staggered lattiudes (cell edges in y) + input_lat_edge(1) = -90.0_r8 + do j1 = 2, num_input_lats + input_lat_edge(j1) = (input_lat(j1 - 1) + input_lat(j1))/2.0_r8 + end do + input_lat_edge(num_input_lats + 1) = 90.0_r8 + + sim_lat_edge(1) = -90.0_r8 + do j2 = 2, num_sim_lats + sim_lat_edge(j2) = (sim_lat(j2 - 1) + sim_lat(j2))/2.0_r8 + end do + sim_lat_edge(num_sim_lats + 1) = 90.0_r8 + + ! Compute the weight for all (input_lon, sim_lon) pairs + do i2 = 1, num_sim_lons + do i1 = 1, num_input_lons + weight_x(i2, i1) = lon_weight(input_lon_edge(i1), input_lon_edge(i1+1), sim_lon_edge(i2), sim_lon_edge(i2+1), use_flight_distance) + end do + end do + + ! Compute the weight for all (input_lat, sim_lat) pairs + do j2 = 1, num_sim_lats + do j1 = 1, num_input_lats + weight_y(j2, j1) = lat_weight(input_lat_edge(j1), input_lat_edge(j1+1), sim_lat_edge(j2), sim_lat_edge(j2+1), use_flight_distance) + end do + end do + + end subroutine xy_interp_init + + ! This program interpolates var_src(im1,jm1,km1) to var_trg(im2,jm2,km1) based on weighting functions weight_x & weight_y. + subroutine xy_interp(im1,jm1,km1,im2,jm2,pcols,ncols,weight_x,weight_y,var_src,var_trg,lons,lats,count_x,count_y,index_x,index_y) + integer, intent(in) :: im1 ! source number of longitudes + integer, intent(in) :: jm1 ! source number of latitudes + integer, intent(in) :: km1 ! source/target number of levels + integer, intent(in) :: im2 ! target number of longitudes + integer, intent(in) :: jm2 ! target number of latitudes + integer, intent(in) :: pcols + integer, intent(in) :: ncols + real(r8), intent(in) :: weight_x(im2, im1), weight_y(jm2, jm1) + real(r8), intent(in) :: var_src(im1, jm1, km1) + integer, intent(in) :: lons(pcols), lats(pcols) + integer, intent(in) :: count_x(im2), count_y(jm2) + integer, intent(in) :: index_x(im2, im1), index_y(jm2, jm1) + real(r8), intent(out) :: var_trg(pcols, km1) + integer :: n, i1, j1, k1, i2, j2, ii, jj + real(r8) :: sum_x + + var_trg(:, :) = 0.0_r8 + + do k1 = 1, km1 + do n = 1, ncols + ! interpolate in x + do jj = 1, count_y(lats(n)) + sum_x = 0.0_r8 + do ii = 1, count_x(lons(n)) + sum_x = sum_x + var_src(index_x(lons(n), ii), index_y(lats(n), jj), k1)* & + weight_x(lons(n), index_x(lons(n), ii)) + end do + var_trg(n, k1) = var_trg(n, k1) + sum_x*weight_y(lats(n), index_y(lats(n), jj)) + end do + end do + end do + + end subroutine xy_interp + +end module horizontal_interpolate diff --git a/src/utils/string_utils.F90 b/src/utils/string_utils.F90 index dbf444d8..c4e6b964 100644 --- a/src/utils/string_utils.F90 +++ b/src/utils/string_utils.F90 @@ -7,6 +7,7 @@ module string_utils use string_core_utils, only: core_int_date_to_yyyymmdd, core_int_seconds_to_hhmmss use string_core_utils, only: to_str => core_to_str use string_core_utils, only: split, stringify, tokenize + use string_core_utils, only: increment_string, last_non_digit, get_last_significant_char implicit none private @@ -20,6 +21,9 @@ module string_utils public :: to_lower ! Convert all characters in string to lower case. public :: split ! Parse a string into tokens, one at a time public :: stringify ! Convert one or more values of any intrinsic data types to a character string for pretty printing + public :: increment_string ! Increment a string whose ending characters are digits. + public :: last_non_digit ! Get position of last non-digit in the input string. + public :: get_last_significant_char ! Get position of last significant (non-blank, non-null) character in string. public :: tokenize ! Parse a string into tokens contains diff --git a/src/utils/tracer_data.F90 b/src/utils/tracer_data.F90 new file mode 100644 index 00000000..46b35911 --- /dev/null +++ b/src/utils/tracer_data.F90 @@ -0,0 +1,2752 @@ +! Module used to read (and interpolate) offline tracer data (sources and mixing ratios) +! Created by: Francis Vitt -- 2 May 2006 +! Modified by: Jim Edwards -- 10 March 2009 +! Modified by: Cheryl Craig and Chih-Chieh (Jack) Chen -- February 2010 +! Dechunkized and ported to SIMA -- Haipeng Lin, October 2025 +! +! Notes on port to SIMA: +! - chunking support has been removed +! - pbuf support is removed, so in_pbuf is no longer supported. +! the previous behavior of in_pbuf meant that if the data was +! already in pbuf, the data arrays in tracer_data were not allocated +! and would directly use the pbuf pointer. Because pbuf is retired in +! SIMA, this is no longer available. The caller should handle the data +! flow in and out of physics schemes instead of using the pbuf, by +! retrieving the data from trfld%data. +! - file removal functionality was removed. If needed, it can be a script +! that runs alongside job submission. +! - FV dycore features: latitude weighting which requires a structured dycore +! is untested as no such dycore is available in SIMA. Polar averaging is +! commented out as only used for FV. +module tracer_data + + use shr_kind_mod, only: r8 => shr_kind_r8, shr_kind_cl + use runtime_obj, only: unset_real + use pio, only: file_desc_t, var_desc_t + + implicit none + private + + public :: trfld, input3d, input2d, trfile + public :: trcdata_init + public :: advance_trcdata + public :: get_fld_data + public :: get_fld_ndx + public :: write_trc_restart + public :: read_trc_restart + public :: init_trc_restart + public :: incr_filename + + type input3d + real(r8), dimension(:, :), allocatable :: data ! ncol, lev + end type input3d + + type input2d + real(r8), dimension(:), allocatable :: data ! ncol + end type input2d + + type trfld + real(r8), dimension(:, :), allocatable :: data ! ncol, lev + type(input3d), dimension(4) :: input + character(len=32) :: srcnam + character(len=32) :: fldnam + character(len=32) :: units + type(var_desc_t) :: var_id + integer :: coords(4) ! LATDIM | LONDIM | LEVDIM | TIMDIM + integer :: order(4) ! LATDIM | LONDIM | LEVDIM | TIMDIM + logical :: srf_fld = .false. + end type trfld + + type trfile + type(input2d), dimension(4) :: ps_in + character(len=shr_kind_cl) :: pathname = ' ' + character(len=shr_kind_cl) :: curr_filename = ' ' + character(len=shr_kind_cl) :: next_filename = ' ' + type(file_desc_t) :: curr_fileid + type(file_desc_t) :: next_fileid + + type(var_desc_t), allocatable :: currfnameid ! pio restart file var id + type(var_desc_t), allocatable :: nextfnameid ! pio restart file var id + + character(len=shr_kind_cl) :: filenames_list = '' + real(r8) :: datatimem = -unset_real ! time of prv. values read in + real(r8) :: datatimep = -unset_real ! time of nxt. values read in + real(r8) :: datatimes(4) + integer :: interp_recs + real(r8), dimension(:), allocatable :: curr_data_times + real(r8), dimension(:), allocatable :: next_data_times + real(r8) :: offset_time + integer :: cyc_ndx_beg + integer :: cyc_ndx_end + integer :: cyc_yr = 0 + real(r8) :: one_yr = 0 + real(r8) :: curr_mod_time ! model time - calendar day + real(r8) :: next_mod_time ! model time - calendar day - next time step + integer :: nlon = 0 + integer :: nlat = 0 + integer :: nlev = 0 + integer :: nilev = 0 + integer :: ps_coords(3) ! LATDIM | LONDIM | TIMDIM + integer :: ps_order(3) ! LATDIM | LONDIM | TIMDIM + real(r8), dimension(:), allocatable :: lons + real(r8), dimension(:), allocatable :: lats + real(r8), dimension(:), allocatable :: levs + real(r8), dimension(:), allocatable :: ilevs + real(r8), dimension(:), allocatable :: hyam + real(r8), dimension(:), allocatable :: hybm + real(r8), dimension(:), allocatable :: hyai + real(r8), dimension(:), allocatable :: hybi + real(r8), dimension(:, :), allocatable :: weight_x, weight_y + integer, dimension(:) , allocatable :: count_x, count_y + integer, dimension(:, :), allocatable :: index_x, index_y + real(r8), dimension(:, :), allocatable :: weight0_x, weight0_y + integer, dimension(:) , allocatable :: count0_x, count0_y + integer, dimension(:, :), allocatable :: index0_x, index0_y + logical :: dist + + real(r8) :: p0 + type(var_desc_t) :: ps_id ! var id of PS variable + logical :: has_ps = .false. + logical :: zonal_ave = .false. + logical :: unstructured = .false. + logical :: alt_data = .false. + logical :: geop_alt = .false. + logical :: cyclical = .false. + logical :: cyclical_list = .false. + logical :: weight_by_lat = .false. + logical :: conserve_column = .false. + logical :: fill_in_months = .false. + logical :: fixed = .false. + logical :: initialized = .false. + logical :: top_bndry = .false. + logical :: top_layer = .false. + logical :: stepTime = .false. ! Do not interpolate in time, but use stepwise times + end type trfile + + integer, parameter :: LONDIM = 1 + integer, parameter :: LATDIM = 2 + integer, parameter :: LEVDIM = 3 + integer, parameter :: TIMDIM = 4 + + integer, parameter :: PS_TIMDIM = 3 + + integer, parameter :: ZA_LATDIM = 1 + integer, parameter :: ZA_LEVDIM = 2 + integer, parameter :: ZA_TIMDIM = 3 + + integer, parameter :: nm = 1 ! array index for previous (minus) data + integer, parameter :: np = 2 ! array index for next (plus) data + + integer, allocatable :: lon_global_grid_ndx(:) ! (ncol) + integer, allocatable :: lat_global_grid_ndx(:) ! (ncol) + +contains + +!-------------------------------------------------------------------------- +!-------------------------------------------------------------------------- + subroutine trcdata_init(specifier, filename, filelist, datapath, flds, file, & + data_cycle_yr, data_fixed_ymd, data_fixed_tod, data_type) + + use vert_coord, only: pver, pverp + use physics_grid, only: pcols => columns_on_task + + use physconst, only: pi + use string_utils, only: to_lower + use cam_abortutils, only: check_allocate, endrun + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + use time_manager, only: set_time_float_from_date + + use pio, only: pio_inquire_variable, pio_inq_vardimid, pio_inq_varid, pio_inq_dimid + use pio, only: pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR + use pio, only: pio_get_att, pio_get_var + + ! For latitude weighting functionality + !use dyn_grid, only: get_horiz_grid_int + !use physics_grid, only: get_rlat_all_p, get_rlon_all_p + !use spmd_utils, only: mpicom, masterprocid, mpi_real8, mpi_integer + !use horizontal_interpolate, only: xy_interp_init + use physics_grid, only: dycore_unstructured + use physics_grid, only: plon => hdim1_d, plat => hdim2_d + + character(len=*), intent(in) :: specifier(:) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: filelist + character(len=*), intent(in) :: datapath + type(trfld), dimension(:), pointer :: flds + type(trfile), intent(inout) :: file + integer, intent(in) :: data_cycle_yr + integer, intent(in) :: data_fixed_ymd + integer, intent(in) :: data_fixed_tod + character(len=*), intent(in) :: data_type + + character(len=*), parameter :: sub = 'trcdata_init' + real(r8), parameter :: d2r = pi/180._r8 + + integer :: f, mxnflds, astat + integer :: str_yr, str_mon, str_day + integer :: lon_dimid, lat_dimid, lev_dimid, tim_dimid, old_dimid + integer :: dimids(4), did + type(var_desc_t) :: varid + integer :: idx + integer :: ierr + integer :: errcode + real(r8) :: start_time, time1, time2 + integer :: i1, i2, j1, j2 + integer :: nvardims, vardimids(4) + + character(len=256) :: data_units + real(r8), allocatable :: lam(:), phi(:) + real(r8) :: rlats(pcols), rlons(pcols) + integer :: ncol, icol, i, j + logical :: found + integer :: aircraft_cnt + integer :: err_handling + character(len=512) :: errmsg + + call specify_fields(specifier, flds) + + file%datatimep = -unset_real + file%datatimem = -unset_real + + mxnflds = 0 + if (associated(flds)) mxnflds = size(flds) + + if (mxnflds < 1) return + + file%pathname = trim(datapath) + file%filenames_list = trim(filelist) + + file%fill_in_months = .false. + file%cyclical = .false. + file%cyclical_list = .false. + + select case (data_type) + case ('FIXED') + file%fixed = .true. + case ('INTERP_MISSING_MONTHS') + file%fill_in_months = .true. + case ('CYCLICAL') + file%cyclical = .true. + file%cyc_yr = data_cycle_yr + case ('CYCLICAL_LIST') + file%cyclical_list = .true. + file%cyc_yr = data_cycle_yr + case ('SERIAL') + ! nothing needs to be done here. + case default + write (iulog, *) sub//': invalid data type: '//trim(data_type)//' file: '//trim(filename) + write (iulog, *) sub//': valid data types: SERIAL | CYCLICAL | CYCLICAL_LIST | FIXED | INTERP_MISSING_MONTHS ' + call endrun(sub//': invalid data type: '//trim(data_type)//' file: '//trim(filename)) + end select + + if ((.not. file%fixed) .and. ((data_fixed_ymd > 0) .or. (data_fixed_tod > 0))) then + call endrun(sub//': Cannot specify data_fixed_ymd or data_fixed_tod if data type is not FIXED') + end if + if ((.not. file%cyclical) .and. (data_cycle_yr > 0)) then + call endrun(sub//': Cannot specify data_cycle_yr if data type is not CYCLICAL') + end if + + if (file%top_bndry .and. file%top_layer) then + call endrun(sub//': Cannot set both file%top_bndry and file%top_layer to TRUE.') + end if + + if (masterproc) then + write (iulog, *) sub//': data type: '//trim(data_type)//' file: '//trim(filename) + end if + + ! if there is no list of files (len_trim(file%filenames_list)<1) then + ! -> set curr_filename from namelist rather than restart data + if (len_trim(file%curr_filename) < 1 .or. len_trim(file%filenames_list) < 1 .or. file%fixed) then ! initial run + file%curr_filename = trim(filename) + + call get_model_time(file) + + if (file%fixed) then + str_yr = data_fixed_ymd/10000 + str_mon = (data_fixed_ymd - str_yr*10000)/100 + str_day = data_fixed_ymd - str_yr*10000 - str_mon*100 + call set_time_float_from_date(start_time, str_yr, str_mon, str_day, data_fixed_tod) + file%offset_time = start_time - file%curr_mod_time + else + file%offset_time = 0 + end if + end if + + call set_time_float_from_date(time2, 2, 1, 1, 0) + call set_time_float_from_date(time1, 1, 1, 1, 0) + file%one_yr = time2 - time1 + + if (file%cyclical .or. file%cyclical_list) then + file%cyc_ndx_beg = -1 + file%cyc_ndx_end = -1 + if (file%cyc_yr /= 0) then + call set_time_float_from_date(time1, file%cyc_yr, 1, 1, 0) + call set_time_float_from_date(time2, file%cyc_yr + 1, 1, 1, 0) + file%one_yr = time2 - time1 + end if + + if (file%cyclical) then + call open_trc_datafile(file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times, & + cyc_ndx_beg=file%cyc_ndx_beg, cyc_ndx_end=file%cyc_ndx_end, cyc_yr=file%cyc_yr) + else + call open_trc_datafile(file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times) + end if + else + call open_trc_datafile(file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times) + file%curr_data_times = file%curr_data_times - file%offset_time + end if + + call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling) + ierr = pio_inq_dimid(file%curr_fileid, 'ncol', idx) + file%unstructured = (ierr == PIO_NOERR) + if (.not. file%unstructured) then + ierr = pio_inq_dimid(file%curr_fileid, 'lon', idx) + file%zonal_ave = (ierr /= PIO_NOERR) + end if + call pio_seterrorhandling(File%curr_fileid, err_handling) + + if (file%zonal_ave) then + file%nlon = 1 + else + if (.not. file%unstructured) then + call get_dimension(file%curr_fileid, 'lon', file%nlon, dimid=old_dimid, data=file%lons) + + file%lons = file%lons*d2r + + lon_dimid = old_dimid + end if + end if + + ierr = pio_inq_dimid(file%curr_fileid, 'time', old_dimid) + + if (.not. file%unstructured) then + ! Hack to work with weird netCDF and old gcc or NAG bug. + tim_dimid = old_dimid + + call get_dimension(file%curr_fileid, 'lat', file%nlat, dimid=old_dimid, data=file%lats) + file%lats = file%lats*d2r + + lat_dimid = old_dimid + end if + + call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling) + ierr = pio_inq_varid(file%curr_fileid, 'PS', file%ps_id) + file%has_ps = (ierr == PIO_NOERR) + ierr = pio_inq_dimid(file%curr_fileid, 'altitude', idx) + file%alt_data = (ierr == PIO_NOERR) + + call pio_seterrorhandling(File%curr_fileid, err_handling) + + if (file%has_ps .and. .not. file%unstructured) then + if (file%zonal_ave) then + ierr = pio_inq_vardimid(file%curr_fileid, file%ps_id, dimids(1:2)) + do did = 1, 2 + if (dimids(did) == lat_dimid) then + file%ps_coords(LATDIM) = did + file%ps_order(did) = LATDIM + else if (dimids(did) == tim_dimid) then + file%ps_coords(PS_TIMDIM) = did + file%ps_order(did) = PS_TIMDIM + end if + end do + else + ierr = pio_inq_vardimid(file%curr_fileid, file%ps_id, dimids(1:3)) + do did = 1, 3 + if (dimids(did) == lon_dimid) then + file%ps_coords(LONDIM) = did + file%ps_order(did) = LONDIM + else if (dimids(did) == lat_dimid) then + file%ps_coords(LATDIM) = did + file%ps_order(did) = LATDIM + else if (dimids(did) == tim_dimid) then + file%ps_coords(PS_TIMDIM) = did + file%ps_order(did) = PS_TIMDIM + end if + end do + end if + end if + + if (masterproc) then + write (iulog, *) sub//': file%has_ps = ', file%has_ps + end if + + if (file%alt_data) then + call get_dimension(file%curr_fileid, 'altitude_int', file%nilev, data=file%ilevs) + call get_dimension(file%curr_fileid, 'altitude', file%nlev, dimid=old_dimid, data=file%levs) + else + call get_dimension(file%curr_fileid, 'lev', file%nlev, dimid=old_dimid, data=file%levs) + if (old_dimid > 0) then + file%levs = file%levs*100._r8 ! mbar->pascals + end if + end if + + ! For some bizarre reason, netCDF with older gcc is keeping a pointer to the dimid, and overwriting it later! + ! Hackish workaround is to make a copy... + lev_dimid = old_dimid + + if (file%has_ps) then + + allocate (file%hyam(file%nlev), file%hybm(file%nlev), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'file%hyam(file%nlev), file%hybm(file%nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + allocate (file%hyai(file%nlev + 1), file%hybi(file%nlev + 1), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'file%hyai(file%nlev + 1), file%hybi(file%nlev + 1)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling) + ierr = pio_inq_varid(file%curr_fileid, 'P0', varid) + call pio_seterrorhandling(File%curr_fileid, err_handling) + + if (ierr == PIO_NOERR) then + ierr = pio_get_var(file%curr_fileid, varid, file%p0) + else + file%p0 = 100000._r8 + end if + ierr = pio_inq_varid(file%curr_fileid, 'hyam', varid) + ierr = pio_get_var(file%curr_fileid, varid, file%hyam) + ierr = pio_inq_varid(file%curr_fileid, 'hybm', varid) + ierr = pio_get_var(file%curr_fileid, varid, file%hybm) + if (file%conserve_column) then + ierr = pio_inq_varid(file%curr_fileid, 'hyai', varid) + ierr = pio_get_var(file%curr_fileid, varid, file%hyai) + ierr = pio_inq_varid(file%curr_fileid, 'hybi', varid) + ierr = pio_get_var(file%curr_fileid, varid, file%hybi) + end if + + allocate (file%ps_in(1)%data(pcols), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'file%ps_in(1)%data(pcols)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + allocate (file%ps_in(2)%data(pcols), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'file%ps_in(2)%data(pcols)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + if (file%fill_in_months) then + allocate (file%ps_in(3)%data(pcols), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'file%ps_in(3)%data(pcols)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + allocate (file%ps_in(4)%data(pcols), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'file%ps_in(4)%data(pcols)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + end if + end if + + call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling) + + flds_loop: do f = 1, mxnflds + + ! get netcdf variable id for the field + ierr = pio_inq_varid(file%curr_fileid, flds(f)%srcnam, flds(f)%var_id) + if (ierr /= pio_noerr) then + call endrun(sub//': Cannot find var "'//trim(flds(f)%srcnam)// & + '" in file "'//trim(file%curr_filename)//'"') + end if + + ! determine if the field has a vertical dimension + + if (lev_dimid > 0) then + ierr = pio_inquire_variable(file%curr_fileid, flds(f)%var_id, ndims=nvardims) + ierr = pio_inquire_variable(file%curr_fileid, flds(f)%var_id, dimids=vardimids(:nvardims)) + flds(f)%srf_fld = .not. any(vardimids(:nvardims) == lev_dimid) + else + flds(f)%srf_fld = .true. + end if + + ! allocate memory for data in container. + if (flds(f)%srf_fld .or. file%top_bndry .or. file%top_layer) then + ! surface/top boundary/top layer field. + allocate (flds(f)%data(pcols, 1), stat=astat, errmsg=errmsg) + else + allocate (flds(f)%data(pcols, pver), stat=astat, errmsg=errmsg) + end if + call check_allocate(astat, sub, 'flds(f)%data(pcols, nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + if (flds(f)%srf_fld) then + allocate (flds(f)%input(1)%data(pcols, 1), stat=astat, errmsg=errmsg) + else + allocate (flds(f)%input(1)%data(pcols, file%nlev), stat=astat, errmsg=errmsg) + end if + call check_allocate(astat, sub, 'flds(f)%input(1)%data(pcols, file%nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + if (flds(f)%srf_fld) then + allocate (flds(f)%input(2)%data(pcols, 1), stat=astat, errmsg=errmsg) + else + allocate (flds(f)%input(2)%data(pcols, file%nlev), stat=astat, errmsg=errmsg) + end if + call check_allocate(astat, sub, 'flds(f)%input(2)%data(pcols, file%nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + if (file%fill_in_months) then + if (flds(f)%srf_fld) then + allocate (flds(f)%input(3)%data(pcols, 1), stat=astat, errmsg=errmsg) + else + allocate (flds(f)%input(3)%data(pcols, file%nlev), stat=astat, errmsg=errmsg) + end if + call check_allocate(astat, sub, 'flds(f)%input(3)%data(pcols, file%nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + if (flds(f)%srf_fld) then + allocate (flds(f)%input(4)%data(pcols, 1), stat=astat, errmsg=errmsg) + else + allocate (flds(f)%input(4)%data(pcols, file%nlev), stat=astat, errmsg=errmsg) + end if + call check_allocate(astat, sub, 'flds(f)%input(4)%data(pcols, file%nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + end if + + if (file%zonal_ave) then + ierr = pio_inq_vardimid(file%curr_fileid, flds(f)%var_id, dimids(1:3)) + do did = 1, 3 + if (dimids(did) == lat_dimid) then + flds(f)%coords(ZA_LATDIM) = did + flds(f)%order(did) = ZA_LATDIM + else if (dimids(did) == lev_dimid) then + flds(f)%coords(ZA_LEVDIM) = did + flds(f)%order(did) = ZA_LEVDIM + else if (dimids(did) == tim_dimid) then + flds(f)%coords(ZA_TIMDIM) = did + flds(f)%order(did) = ZA_TIMDIM + end if + end do + else if (flds(f)%srf_fld .and. .not. file%unstructured) then + ierr = pio_inq_vardimid(file%curr_fileid, flds(f)%var_id, dimids(1:3)) + do did = 1, 3 + if (dimids(did) == lon_dimid) then + flds(f)%coords(LONDIM) = did + flds(f)%order(did) = LONDIM + else if (dimids(did) == lat_dimid) then + flds(f)%coords(LATDIM) = did + flds(f)%order(did) = LATDIM + else if (dimids(did) == tim_dimid) then + flds(f)%coords(PS_TIMDIM) = did + flds(f)%order(did) = PS_TIMDIM + end if + end do + else if (.not. file%unstructured) then + ierr = pio_inq_vardimid(file%curr_fileid, flds(f)%var_id, dimids) + do did = 1, 4 + if (dimids(did) == lon_dimid) then + flds(f)%coords(LONDIM) = did + flds(f)%order(did) = LONDIM + else if (dimids(did) == lat_dimid) then + flds(f)%coords(LATDIM) = did + flds(f)%order(did) = LATDIM + else if (dimids(did) == lev_dimid) then + flds(f)%coords(LEVDIM) = did + flds(f)%order(did) = LEVDIM + else if (dimids(did) == tim_dimid) then + flds(f)%coords(TIMDIM) = did + flds(f)%order(did) = TIMDIM + end if + end do + end if + + ! retrieve units from file (used by downstream code to potentially + ! perform unit conversions on read data) + ! + ! convert units to lowercase to facilitate comparisons + ierr = pio_get_att(file%curr_fileid, flds(f)%var_id, 'units', data_units) + flds(f)%units = trim(to_lower(data_units(1:32))) + + end do flds_loop + + call pio_seterrorhandling(File%curr_fileid, err_handling) + + ! if weighting by latitude, compute weighting for horizontal interpolation + if (file%weight_by_lat) then + if (dycore_unstructured) then + call endrun(sub//': weighting by latitude not implemented for unstructured grids') + end if + + call endrun(sub//': weighting by latitude (used by aircraft emis) is untested in SIMA; uncomment this line for testing.') + ! WARNING: in SIMA, currently implemented dycores are unstructured. + ! The below code has been ported to the best of ability, + ! but is completely untested. (hplin, 10/9/25) + + ! allocate (lam(plon), phi(plat)) + ! call get_horiz_grid_int(plat, clat_d_out=phi) + ! call get_horiz_grid_int(plon, clon_d_out=lam) + + ! if (.not. allocated(lon_global_grid_ndx)) allocate (lon_global_grid_ndx(pcols)) + ! if (.not. allocated(lat_global_grid_ndx)) allocate (lat_global_grid_ndx(pcols)) + ! lon_global_grid_ndx = -huge(1) + ! lat_global_grid_ndx = -huge(1) + + ! ncol = pcols ! active columns + ! call get_rlat_all_p(ncol, rlats(:ncol)) + ! call get_rlon_all_p(ncol, rlons(:ncol)) + ! do icol = 1, ncol + ! found = .false. + ! find_col: do j = 1, plat + ! do i = 1, plon + ! if (rlats(icol) == phi(j) .and. rlons(icol) == lam(i)) then + ! found = .true. + ! exit find_col + ! end if + ! end do + ! end do find_col + + ! if (.not. found) call endrun(sub//': not able to find physics column coordinate') + ! lon_global_grid_ndx(icol) = i + ! lat_global_grid_ndx(icol) = j + ! end do + + ! deallocate (phi, lam) + + ! ! weight_x & weight_y are weighting function for x & y interpolation + ! allocate (file%weight_x(plon, file%nlon), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%weight_x allocation error = ', astat + ! call endrun(sub//': failed to allocate weight_x array') + ! end if + ! allocate (file%weight_y(plat, file%nlat), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%weight_y allocation error = ', astat + ! call endrun(sub//': failed to allocate weight_y array') + ! end if + ! allocate (file%count_x(plon), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%count_x allocation error = ', astat + ! call endrun(sub//': failed to allocate count_x array') + ! end if + ! allocate (file%count_y(plat), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%count_y allocation error = ', astat + ! call endrun(sub//': failed to allocate count_y array') + ! end if + ! allocate (file%index_x(plon, file%nlon), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%index_x allocation error = ', astat + ! call endrun(sub//': failed to allocate index_x array') + ! end if + ! allocate (file%index_y(plat, file%nlat), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%index_y allocation error = ', astat + ! call endrun(sub//': failed to allocate index_y array') + ! end if + ! file%weight_x(:, :) = 0.0_r8 + ! file%weight_y(:, :) = 0.0_r8 + ! file%count_x(:) = 0 + ! file%count_y(:) = 0 + ! file%index_x(:, :) = 0 + ! file%index_y(:, :) = 0 + + ! if (file%dist) then + ! allocate (file%weight0_x(plon, file%nlon), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%weight0_x allocation error = ', astat + ! call endrun(sub//': failed to allocate weight0_x array') + ! end if + ! allocate (file%weight0_y(plat, file%nlat), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%weight0_y allocation error = ', astat + ! call endrun(sub//': failed to allocate weight0_y array') + ! end if + ! allocate (file%count0_x(plon), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%count0_x allocation error = ', astat + ! call endrun(sub//': failed to allocate count0_x array') + ! end if + ! allocate (file%count0_y(plat), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%count0_y allocation error = ', astat + ! call endrun(sub//': failed to allocate count0_y array') + ! end if + ! allocate (file%index0_x(plon, file%nlon), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%index0_x allocation error = ', astat + ! call endrun(sub//': failed to allocate index0_x array') + ! end if + ! allocate (file%index0_y(plat, file%nlat), stat=astat) + ! if (astat /= 0) then + ! write (iulog, *) sub//': file%index0_y allocation error = ', astat + ! call endrun(sub//': failed to allocate index0_y array') + ! end if + ! file%weight0_x(:, :) = 0.0_r8 + ! file%weight0_y(:, :) = 0.0_r8 + ! file%count0_x(:) = 0 + ! file%count0_y(:) = 0 + ! file%index0_x(:, :) = 0 + ! file%index0_y(:, :) = 0 + ! end if + + ! if (masterproc) then + ! ! compute weighting. NOTE: we always set + ! ! use_flight_distance=.false. for this path since these + ! ! weights are used to inerpolate field values like PS even + ! ! when the file contains other data which should be treated + ! ! as per-cell totals. + ! call xy_interp_init(file%nlon, file%nlat, file%lons, file%lats, & + ! plon, plat, file%weight_x, file%weight_y, .false.) + + ! do i2 = 1, plon + ! file%count_x(i2) = 0 + ! do i1 = 1, file%nlon + ! if (file%weight_x(i2, i1) > 0.0_r8) then + ! file%count_x(i2) = file%count_x(i2) + 1 + ! file%index_x(i2, file%count_x(i2)) = i1 + ! end if + ! end do + ! end do + + ! do j2 = 1, plat + ! file%count_y(j2) = 0 + ! do j1 = 1, file%nlat + ! if (file%weight_y(j2, j1) > 0.0_r8) then + ! file%count_y(j2) = file%count_y(j2) + 1 + ! file%index_y(j2, file%count_y(j2)) = j1 + ! end if + ! end do + ! end do + + ! if (file%dist) then + ! call xy_interp_init(file%nlon, file%nlat, file%lons, file%lats, & + ! plon, plat, file%weight0_x, file%weight0_y, .true.) + + ! do i2 = 1, plon + ! file%count0_x(i2) = 0 + ! do i1 = 1, file%nlon + ! if (file%weight0_x(i2, i1) > 0.0_r8) then + ! file%count0_x(i2) = file%count0_x(i2) + 1 + ! file%index0_x(i2, file%count0_x(i2)) = i1 + ! end if + ! end do + ! end do + + ! do j2 = 1, plat + ! file%count0_y(j2) = 0 + ! do j1 = 1, file%nlat + ! if (file%weight0_y(j2, j1) > 0.0_r8) then + ! file%count0_y(j2) = file%count0_y(j2) + 1 + ! file%index0_y(j2, file%count0_y(j2)) = j1 + ! end if + ! end do + ! end do + + ! end if + ! end if + + ! call mpi_bcast(file%weight_x, plon*file%nlon, mpi_real8, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_x") + ! call mpi_bcast(file%weight_y, plat*file%nlat, mpi_real8, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_y") + ! call mpi_bcast(file%count_x, plon, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_x") + ! call mpi_bcast(file%count_y, plat, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_y") + ! call mpi_bcast(file%index_x, plon*file%nlon, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index_x") + ! call mpi_bcast(file%index_y, plat*file%nlat, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index_y") + ! if (file%dist) then + ! call mpi_bcast(file%weight0_x, plon*file%nlon, mpi_real8, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_x") + ! call mpi_bcast(file%weight0_y, plat*file%nlat, mpi_real8, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_y") + ! call mpi_bcast(file%count0_x, plon, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_x") + ! call mpi_bcast(file%count0_y, plat, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_y") + ! call mpi_bcast(file%index0_x, plon*file%nlon, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_x") + ! call mpi_bcast(file%index0_y, plat*file%nlat, mpi_integer, masterprocid, mpicom, ierr) + ! if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_y") + ! end if + + end if + end subroutine trcdata_init + + !----------------------------------------------------------------------- + ! Reads more data if needed and interpolates data to current model time + !----------------------------------------------------------------------- + subroutine advance_trcdata( & + flds, file, & + pmid, pint, phis, zi) + + use ccpp_kinds, only: kind_phys + use perf_mod, only: t_startf, t_stopf + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + + ! dimensions of the grid can be retrieved directly + use vert_coord, only: pver, pverp + use physics_grid, only: ncol => columns_on_task + + type(trfile), intent(inout) :: file + type(trfld), intent(inout) :: flds(:) + + ! state inputs for interpolating to model grid. + real(kind_phys), intent(in) :: pmid(:,:) ! air pressure [Pa] + real(kind_phys), intent(in) :: pint(:,:) ! air pressure at interfaces [Pa] + real(kind_phys), intent(in) :: phis(:) ! surface geopotential [m2 s-2] + real(kind_phys), intent(in) :: zi(:,:) ! height above surface, interfaces [m] + + real(r8) :: data_time + + call t_startf('advance_trcdata') + if (.not. (file%fixed .and. file%initialized)) then + + call get_model_time(file) + + data_time = file%datatimep + + if (file%cyclical .or. file%cyclical_list) then + ! wrap around + if ((file%datatimep < file%datatimem) .and. (file%curr_mod_time > file%datatimem)) then + data_time = data_time + file%one_yr + end if + end if + + ! For stepTime need to advance if the times are equal + ! Should not impact other runs? + if (file%curr_mod_time >= data_time) then + call t_startf('read_next_trcdata') + call read_next_trcdata(flds, file) + call t_stopf('read_next_trcdata') + if (masterproc) write (iulog, *) 'READ_NEXT_TRCDATA: ', flds%fldnam + end if + + file%initialized = .true. + + end if + + ! need to interpolate the data, regardless + ! each mpi task needs to interpolate + call t_startf('interpolate_trcdata') + call interpolate_trcdata(ncol = ncol, & + pver = pver, & + pverp = pverp, & + pmid = pmid(:ncol,:pver), & + pint = pint(:ncol,:pverp), & + phis = phis(:ncol), & + zi = zi(:ncol, :pverp), & + flds = flds(:), & + file = file) + call t_stopf('interpolate_trcdata') + + call t_stopf('advance_trcdata') + + end subroutine advance_trcdata + + pure subroutine get_fld_data(flds, field_name, data, ncol) + type(trfld), intent(inout) :: flds(:) + character(len=*), intent(in) :: field_name + real(r8), intent(out) :: data(:, :) + integer, intent(in) :: ncol + + integer :: f, nflds + + data(:, :) = 0._r8 + nflds = size(flds) + + do f = 1, nflds + if (trim(flds(f)%fldnam) == trim(field_name)) then + data(:ncol, :) = flds(f)%data(:ncol, :) + end if + end do + + end subroutine get_fld_data + + pure subroutine get_fld_ndx(flds, field_name, idx) + + type(trfld), intent(in) :: flds(:) + character(len=*), intent(in) :: field_name + integer, intent(out) :: idx + integer :: f, nflds + + idx = -1 + nflds = size(flds) + + do f = 1, nflds + if (trim(flds(f)%fldnam) == trim(field_name)) then + idx = f + return + end if + end do + + end subroutine get_fld_ndx + + subroutine get_model_time(file) + use time_manager, only: set_time_float_from_date + use time_manager, only: get_curr_date + use time_manager, only: get_step_size + + type(trfile), intent(inout) :: file + + integer :: yr, mon, day, ncsec ! components of a date + + call get_curr_date(yr, mon, day, ncsec) + + if (file%cyclical .or. file%cyclical_list) yr = file%cyc_yr + call set_time_float_from_date(file%curr_mod_time, yr, mon, day, ncsec) + file%next_mod_time = file%curr_mod_time + get_step_size()/86400._r8 + + end subroutine get_model_time + + subroutine check_files(file, fids, itms, times_found) + type(trfile), intent(inout) :: file + type(file_desc_t), intent(out) :: fids(2) ! ids of files that contains these recs + integer, optional, intent(out) :: itms(2) + logical, optional, intent(inout) :: times_found + + logical :: list_cycled + + list_cycled = .false. + + !----------------------------------------------------------------------- + ! If next time beyond the end of the time list, + ! then increment the filename and move on to the next file + !----------------------------------------------------------------------- + if ((file%next_mod_time > file%curr_data_times(size(file%curr_data_times))) .or. file%cyclical_list) then + if (file%cyclical_list) then + if (allocated(file%next_data_times)) then + if ((file%curr_mod_time > file%datatimep)) then + call advance_file(file) + end if + end if + end if + + if (.not. allocated(file%next_data_times)) then + ! open next file if not already opened... + if (file%cyclical_list) then + file%next_filename = incr_filename(file%curr_filename, filenames_list=file%filenames_list, datapath=file%pathname, & + cyclical_list=file%cyclical_list, list_cycled=list_cycled) + else + file%next_filename = incr_filename(file%curr_filename, filenames_list=file%filenames_list, datapath=file%pathname) + end if + call open_trc_datafile(file%next_filename, file%pathname, file%next_fileid, file%next_data_times) + file%next_data_times = file%next_data_times - file%offset_time + end if + end if + + !----------------------------------------------------------------------- + ! If using next_data_times and the current is greater than or equal to the next, then + ! close the current file, and set up for next file. + !----------------------------------------------------------------------- + if (allocated(file%next_data_times)) then + if (file%cyclical_list .and. list_cycled) then ! special case - list cycled + file%datatimem = file%curr_data_times(size(file%curr_data_times)) + itms(1) = size(file%curr_data_times) + fids(1) = file%curr_fileid + + file%datatimep = file%next_data_times(1) + itms(2) = 1 + fids(2) = file%next_fileid + + times_found = .true. + else if (file%curr_mod_time >= file%next_data_times(1)) then + call advance_file(file) + end if + end if + + end subroutine check_files + + ! Increment or decrement a date string withing a filename + ! the filename date section is assumed to be of the form yyyy-dd-mm + function incr_filename(filename, filenames_list, datapath, cyclical_list, list_cycled, abort) + use spmd_utils, only: masterproc + use cam_abortutils, only: endrun + use cam_logfile, only: iulog + + use string_utils, only: increment_string + use shr_kind_mod, only: shr_kind_cm + + character(len=*), intent(in) :: filename ! present dynamical dataset filename + character(len=*), optional, intent(in) :: filenames_list + character(len=*), optional, intent(in) :: datapath + logical, optional, intent(in) :: cyclical_list ! If true, allow list to cycle + logical, optional, intent(out) :: list_cycled + logical, optional, intent(in) :: abort + + character(len=shr_kind_cl) :: incr_filename ! next filename in the sequence + + !----------------------------------------------------------------------- + ! ... local variables + !----------------------------------------------------------------------- + integer :: pos, istat + character(len=shr_kind_cl) :: fn_new, line, filepath + integer :: ios, unitnumber + logical :: abort_run + character(len=shr_kind_cm) :: errmsg + + character(len=*), parameter :: sub = 'incr_filename' + + if (present(abort)) then + abort_run = abort + else + abort_run = .true. + end if + + if (present(list_cycled)) list_cycled = .false. + + if ((.not. present(filenames_list)) .or. (len_trim(filenames_list) == 0)) then + !----------------------------------------------------------------------- + ! ... ccm type filename + !----------------------------------------------------------------------- + pos = len_trim(filename) + fn_new = filename(:pos) + if (masterproc) write (iulog, *) sub//': old filename = ', trim(fn_new) + if (fn_new(pos - 2:) == '.nc') then + pos = pos - 3 + end if + istat = increment_string(fn_new(:pos), 1) + if (istat /= 0) then + write (iulog, *) sub//': increment_string returned ', istat + write (iulog, *) ' while trying to decrement ', trim(fn_new) + call endrun(sub//': increment_string failure') + end if + + else + + !------------------------------------------------------------------- + ! ... open filenames_list + !------------------------------------------------------------------- + if (masterproc) write(iulog, *) sub//': old filename = ', trim(filename) + if (masterproc) write(iulog, *) sub//': open filenames_list : ', trim(filenames_list) + + if (present(datapath)) then + filepath = trim(datapath)//'/'//trim(filenames_list) + else + filepath = trim(filenames_list) + end if + + open (newunit=unitnumber, file=filepath, iostat=ios, iomsg=errmsg, status="OLD") + if (ios /= 0) then + call endrun('not able to open file: '//trim(filepath)//'; error = '// trim(errmsg)) + end if + + !------------------------------------------------------------------- + ! ... read file names + !------------------------------------------------------------------- + read (unit=unitnumber, fmt='(A)', iostat=ios, iomsg=errmsg) line + if (ios /= 0) then + if (abort_run) then + call endrun(sub//': not able to increment file name from filenames_list file: '//trim(filenames_list)//'; error = ' // trim(errmsg)) + else + fn_new = 'NOT_FOUND' + incr_filename = trim(fn_new) + return + end if + end if + + !------------------------------------------------------------------- + ! If current filename is '', then initialize with the first filename read in + ! and skip this section. + !------------------------------------------------------------------- + if (filename /= '') then + + !------------------------------------------------------------------- + ! otherwise read until find current filename + !------------------------------------------------------------------- + do while (trim(line) /= trim(filename)) + read (unit=unitnumber, fmt='(A)', iostat=ios, iomsg=errmsg) line + if (ios /= 0) then + if (abort_run) then + call endrun(sub//': not able to increment file name from filenames_list file: '//trim(filenames_list)//'; error = ' // trim(errmsg)) + else + fn_new = 'NOT_FOUND' + incr_filename = trim(fn_new) + return + end if + end if + end do + + !------------------------------------------------------------------- + ! Read next filename + !------------------------------------------------------------------- + read (unit=unitnumber, fmt='(A)', iostat=ios, iomsg=errmsg) line + + !--------------------------------------------------------------------------------- + ! If cyclical_list, then an end of file is not an error, but rather + ! a signal to rewind and start over + !--------------------------------------------------------------------------------- + + if (ios /= 0) then + if (present(cyclical_list)) then + if (cyclical_list) then + list_cycled = .true. + rewind (unitnumber) + read (unit=unitnumber, fmt='(A)', iostat=ios, iomsg=errmsg) line + ! Error here should never happen, but check just in case + if (ios /= 0) then + call endrun(sub//': not able to increment file name from filenames_list file: '//trim(filenames_list)//'; error = '//trim(errmsg)) + end if + else + call endrun(sub//': not able to increment file name from filenames_list file: '//trim(filenames_list)//'; error = ' //trim(errmsg)) + end if + else + if (abort_run) then + call endrun(sub//': not able to increment file name from filenames_list file: '//trim(filenames_list)//'; error = '//trim(errmsg)) + else + fn_new = 'NOT_FOUND' + incr_filename = trim(fn_new) + return + end if + end if + end if + + end if + + !--------------------------------------------------------------------------------- + ! Assign the current filename and close the filelist + !--------------------------------------------------------------------------------- + fn_new = trim(line) + + close (unit=unitnumber) + end if + + !--------------------------------------------------------------------------------- + ! return the current filename + !--------------------------------------------------------------------------------- + incr_filename = trim(fn_new) + if (masterproc) write (iulog, *) sub//': new filename = ', trim(incr_filename) + + end function incr_filename + + subroutine find_times(itms, fids, time, file, datatimem, datatimep, times_found) + use cam_abortutils, only: check_allocate, endrun + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + use shr_kind_mod, only: shr_kind_cm + + integer, intent(out) :: itms(2) ! record numbers that bracket time + type(file_desc_t), intent(out) :: fids(2) ! ids of files that contains these recs + real(r8), intent(in) :: time ! time of interest + type(trfile), intent(in) :: file + real(r8), intent(out) :: datatimem + real(r8), intent(out) :: datatimep + logical, intent(inout) :: times_found + + integer :: np1 ! current forward time index of dataset + integer :: n, i ! + integer :: curr_tsize, next_tsize, all_tsize + integer :: astat + integer :: cyc_tsize + + real(r8), allocatable, dimension(:) :: all_data_times + + character(len=shr_kind_cm) :: errmsg + character(len=*), parameter :: subname = "find_times" + + curr_tsize = size(file%curr_data_times) + next_tsize = 0 + if (allocated(file%next_data_times)) next_tsize = size(file%next_data_times) + + all_tsize = curr_tsize + next_tsize + + allocate (all_data_times(all_tsize), stat=astat, errmsg=errmsg) + call check_allocate(astat, subname, 'all_data_times(all_tsize)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + all_data_times(:curr_tsize) = file%curr_data_times(:) + if (next_tsize > 0) all_data_times(curr_tsize + 1:all_tsize) = file%next_data_times(:) + + if (.not. file%cyclical) then + if (all(all_data_times(:) > time)) then + write (iulog, *) subname//': ALL data times are after ', time + write (iulog, *) subname//': file: ', trim(file%curr_filename) + write (iulog, *) subname//': time: ', time + call endrun(subname//': all(all_data_times(:) > time) '//trim(file%curr_filename)) + end if + + ! find bracketing times + find_times_loop: do n = 1, all_tsize - 1 + np1 = n + 1 + datatimem = all_data_times(n) !+ file%offset_time + datatimep = all_data_times(np1) !+ file%offset_time + ! When stepTime, datatimep may not equal the time (as only datatimem is used) + ! Should not break other runs? + if ((time >= datatimem) .and. (time < datatimep)) then + times_found = .true. + exit find_times_loop + end if + end do find_times_loop + + else ! file%cyclical + + cyc_tsize = file%cyc_ndx_end - file%cyc_ndx_beg + 1 + + if (cyc_tsize > 1) then + + call findplb(all_data_times(file%cyc_ndx_beg:file%cyc_ndx_end), cyc_tsize, time, n) + + if (n == cyc_tsize) then + np1 = 1 + else + np1 = n + 1 + end if + + datatimem = all_data_times(n + file%cyc_ndx_beg - 1) + datatimep = all_data_times(np1 + file%cyc_ndx_beg - 1) + times_found = .true. + + end if + end if + + if (.not. times_found) then + if (masterproc) then + write (iulog, *) subname//': Failed to find dates bracketing desired time =', time + write (iulog, *) 'filename = '//trim(file%curr_filename) + write (iulog, *) ' datatimem = ', file%datatimem + write (iulog, *) ' datatimep = ', file%datatimep + end if + return + end if + + deallocate (all_data_times, stat=astat) + if (astat /= 0) then + write (iulog, *) subname//': failed to deallocate all_data_times array; error = ', astat + call endrun(subname//': failed to deallocate all_data_times array') + end if + + if (.not. file%cyclical) then + itms(1) = n + itms(2) = np1 + else + itms(1) = n + file%cyc_ndx_beg - 1 + itms(2) = np1 + file%cyc_ndx_beg - 1 + end if + + fids(:) = file%curr_fileid + + do i = 1, 2 + if (itms(i) > curr_tsize) then + itms(i) = itms(i) - curr_tsize + fids(i) = file%next_fileid + end if + end do + + end subroutine find_times + + subroutine read_next_trcdata(flds, file) + use cam_abortutils, only: endrun + + use time_manager, only: set_time_float_from_date, set_date_from_time_float + use time_manager, only: get_curr_date + + type(trfile), intent(inout) :: file + type(trfld), intent(inout) :: flds(:) + + integer :: recnos(4), i, f, nflds + integer :: cnt4(4) ! array of counts for each dimension + integer :: strt4(4) ! array of starting indices + integer :: cnt3(3) ! array of counts for each dimension + integer :: strt3(3) ! array of starting indices + type(file_desc_t) :: fids(4) + logical :: times_found + + integer :: cur_yr, cur_mon, cur_day, cur_sec, yr1, yr2, mon, date, sec + real(r8) :: series1_time, series2_time + type(file_desc_t) :: fid1, fid2 + + nflds = size(flds) + times_found = .false. + + do while (.not. times_found) + call find_times(recnos, fids, file%curr_mod_time, file, file%datatimem, file%datatimep, times_found) + if (.not. times_found) then + call check_files(file, fids, recnos, times_found) + end if + end do + + !-------------------------------------------------------------- + ! If stepTime, then no time interpolation is to be done + !-------------------------------------------------------------- + if (file%stepTime) then + file%interp_recs = 1 + else + file%interp_recs = 2 + end if + + if (file%fill_in_months) then + + if (file%datatimep - file%datatimem > file%one_yr) then + + call get_curr_date(cur_yr, cur_mon, cur_day, cur_sec) + + call set_date_from_time_float(file%datatimem, yr1, mon, date, sec) + call set_date_from_time_float(file%datatimep, yr2, mon, date, sec) + + call set_time_float_from_date(series1_time, yr1, cur_mon, cur_day, cur_sec) + call set_time_float_from_date(series2_time, yr2, cur_mon, cur_day, cur_sec) + + fid1 = fids(1) + fid2 = fids(2) + file%cyclical = .true. + call set_cycle_indices(fid1, file%cyc_ndx_beg, file%cyc_ndx_end, yr1) + call find_times(recnos(1:2), fids(1:2), series1_time, file, file%datatimes(1), file%datatimes(2), times_found) + + if (.not. times_found) then + call endrun('read_next_trcdata: time not found for series1_time') + end if + call set_cycle_indices(fid2, file%cyc_ndx_beg, file%cyc_ndx_end, yr2) + + if (fid1%fh /= fid2%fh) then + file%cyc_ndx_beg = file%cyc_ndx_beg + size(file%curr_data_times) + file%cyc_ndx_end = file%cyc_ndx_end + size(file%curr_data_times) + end if + call find_times(recnos(3:4), fids(3:4), series2_time, file, file%datatimes(3), file%datatimes(4), times_found) + if (.not. times_found) then + call endrun('read_next_trcdata: time not found for series2_time') + end if + file%cyclical = .false. + file%interp_recs = 4 + + call set_date_from_time_float(file%datatimes(1), yr1, mon, date, sec) + call set_time_float_from_date(file%datatimem, cur_yr, mon, date, sec) + if (file%datatimes(1) > file%datatimes(2)) then ! wrap around + if (cur_mon == 1) then + call set_time_float_from_date(file%datatimem, cur_yr - 1, mon, date, sec) + end if + end if + + call set_date_from_time_float(file%datatimes(2), yr1, mon, date, sec) + call set_time_float_from_date(file%datatimep, cur_yr, mon, date, sec) + if (file%datatimes(1) > file%datatimes(2)) then ! wrap around + if (cur_mon == 12) then + call set_time_float_from_date(file%datatimep, cur_yr + 1, mon, date, sec) + end if + end if + + end if + + end if + + ! + ! Set up hyperslab corners + ! + do i = 1, file%interp_recs + + strt4(:) = 1 + strt3(:) = 1 + + do f = 1, nflds + if (file%zonal_ave) then + cnt3(flds(f)%coords(ZA_LATDIM)) = file%nlat + if (flds(f)%srf_fld) then + cnt3(flds(f)%coords(ZA_LEVDIM)) = 1 + else + cnt3(flds(f)%coords(ZA_LEVDIM)) = file%nlev + end if + cnt3(flds(f)%coords(ZA_TIMDIM)) = 1 + strt3(flds(f)%coords(ZA_TIMDIM)) = recnos(i) + call read_za_trc(fids(i), flds(f)%var_id, flds(f)%input(i)%data, strt3, cnt3, file, & + (/flds(f)%order(ZA_LATDIM), flds(f)%order(ZA_LEVDIM)/)) + else if (flds(f)%srf_fld) then + if (file%unstructured) then + ! read data directly onto the unstructureed phys grid -- assumes input data is on same grid as phys + call read_physgrid_2d(fids(i), flds(f)%fldnam, recnos(i), flds(f)%input(i)%data(:, 1)) + else + cnt3(flds(f)%coords(LONDIM)) = file%nlon + cnt3(flds(f)%coords(LATDIM)) = file%nlat + cnt3(flds(f)%coords(PS_TIMDIM)) = 1 + strt3(flds(f)%coords(PS_TIMDIM)) = recnos(i) + call read_2d_trc(fids(i), flds(f)%var_id, flds(f)%input(i)%data(:, 1), strt3, cnt3, file, & + (/flds(f)%order(LONDIM), flds(f)%order(LATDIM)/)) + end if + else + if (file%unstructured) then + ! read data directly onto the unstructureed phys grid -- assumes input data is on same grid as phys + if (file%alt_data) then + call read_physgrid_3d(fids(i), flds(f)%fldnam, 'altitude', file%nlev, recnos(i), flds(f)%input(i)%data(:, :)) + else + call read_physgrid_3d(fids(i), flds(f)%fldnam, 'lev', file%nlev, recnos(i), flds(f)%input(i)%data(:, :)) + end if + else + cnt4(flds(f)%coords(LONDIM)) = file%nlon + cnt4(flds(f)%coords(LATDIM)) = file%nlat + cnt4(flds(f)%coords(LEVDIM)) = file%nlev + cnt4(flds(f)%coords(TIMDIM)) = 1 + strt4(flds(f)%coords(TIMDIM)) = recnos(i) + call read_3d_trc(fids(i), flds(f)%var_id, flds(f)%input(i)%data, strt4, cnt4, file, & + (/flds(f)%order(LONDIM), flds(f)%order(LATDIM), flds(f)%order(LEVDIM)/)) + end if + + end if + + end do + + if (file%has_ps) then + if (file%unstructured) then + call read_physgrid_2d(fids(i), 'PS', recnos(i), file%ps_in(i)%data) + else + cnt3 = 1 + strt3 = 1 + if (.not. file%zonal_ave) then + cnt3(file%ps_coords(LONDIM)) = file%nlon + end if + cnt3(file%ps_coords(LATDIM)) = file%nlat + cnt3(file%ps_coords(PS_TIMDIM)) = 1 + strt3(file%ps_coords(PS_TIMDIM)) = recnos(i) + if (file%zonal_ave) then + call read_2d_trc(fids(i), file%ps_id, file%ps_in(i)%data, strt3(1:2), cnt3(1:2), file, & + (/1, 2/)) + else + call read_2d_trc(fids(i), file%ps_id, file%ps_in(i)%data, strt3, cnt3, file, & + (/file%ps_order(LONDIM), file%ps_order(LATDIM)/)) + end if + end if + end if + + end do + + end subroutine read_next_trcdata + +!------------------------------------------------------------------------ + + subroutine read_2d_trc(fid, vid, loc_arr, strt, cnt, file, order) + use physics_grid, only: pcols => columns_on_task + use physics_grid, only: plon => hdim1_d, plat => hdim2_d + + use interpolate_data, only: lininterp_init, lininterp, interp_type, lininterp_finish + use horizontal_interpolate, only: xy_interp + + use physics_grid, only: get_rlat_all_p, get_rlon_all_p + use physconst, only: pi + use cam_abortutils, only: check_allocate + use perf_mod, only: t_startf, t_stopf + + use pio, only: pio_get_var + + type(file_desc_t), intent(in) :: fid + type(var_desc_t), intent(in) :: vid + integer, intent(in) :: strt(:), cnt(:), order(2) + real(r8), intent(out) :: loc_arr(:) ! (ncol) + type(trfile), intent(in) :: file + + real(r8) :: to_lats(pcols), to_lons(pcols) + real(r8), allocatable, target :: wrk2d(:, :) ! (cnt(1), cnt(2)) + real(r8), pointer :: wrk2d_in(:, :) ! (file%nlon, file%nlat) + + integer :: ierr + real(r8), parameter :: zero = 0_r8, twopi = 2_r8*pi + type(interp_type) :: lon_wgts, lat_wgts + integer :: lons(pcols), lats(pcols) + real(r8) :: file_lats(file%nlat) + + character(len=512) :: errmsg + character(len=*), parameter :: subname = "read_2d_trc" + + nullify (wrk2d_in) + allocate (wrk2d(cnt(1), cnt(2)), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'wrk2d(cnt(1), cnt(2))', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + if (order(1) /= 1 .or. order(2) /= 2 .or. cnt(1) /= file%nlon .or. cnt(2) /= file%nlat) then + allocate (wrk2d_in(file%nlon, file%nlat), stat=ierr, errmsg=errmsg) + + call check_allocate(ierr, subname, 'wrk2d_in(file%nlon, file%nlat)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + end if + + ierr = pio_get_var(fid, vid, strt, cnt, wrk2d) + if (associated(wrk2d_in)) then + wrk2d_in = reshape(wrk2d(:, :), (/file%nlon, file%nlat/), order=order) + deallocate (wrk2d) + else + wrk2d_in => wrk2d + end if + + ! PGI 13.9 bug workaround. + file_lats = file%lats + + ! For zonal average, only interpolate along latitude. + if (file%zonal_ave) then + call get_rlat_all_p(pcols, to_lats) + + call lininterp_init(file_lats, file%nlat, to_lats, pcols, 1, lat_wgts) + + call lininterp(wrk2d_in(1, :), file%nlat, loc_arr(1:pcols), pcols, lat_wgts) + + call lininterp_finish(lat_wgts) + else + ! if weighting by latitude, the perform horizontal interpolation by using weight_x, weight_y + + if (file%weight_by_lat) then + + call t_startf('xy_interp') + + lons(:pcols) = lon_global_grid_ndx(:pcols) + lats(:pcols) = lat_global_grid_ndx(:pcols) + + ! NOTE: This uses weight_[xy] instead of weight0_[xy] and + ! hence treats the values as a field rather than per-cell + ! totals. When file%dist == TRUE, this path only appears + ! to be used to interpolate PS, which is probably the + ! correct behavior. + ! + ! @reviewers: The control flow is convoluted here, so + ! this merits some additional scrutiny. + ! + ! in SIMA: pcols (size of array) and ncols (loop dim) here now equal + call xy_interp(file%nlon, file%nlat, 1, plon, plat, pcols, pcols, & + file%weight_x, file%weight_y, wrk2d_in, loc_arr(:), & + lons, lats, file%count_x, file%count_y, file%index_x, file%index_y) + + call t_stopf('xy_interp') + + else + call get_rlat_all_p(pcols, to_lats) + call get_rlon_all_p(pcols, to_lons) + + call lininterp_init(file%lons, file%nlon, to_lons, pcols, 2, lon_wgts, zero, twopi) + call lininterp_init(file%lats, file%nlat, to_lats, pcols, 1, lat_wgts) + + call lininterp(wrk2d_in, file%nlon, file%nlat, loc_arr(1:pcols), pcols, lon_wgts, lat_wgts) + + call lininterp_finish(lon_wgts) + call lininterp_finish(lat_wgts) + end if + + end if + + if (allocated(wrk2d)) then + deallocate(wrk2d) + else + deallocate(wrk2d_in) + end if + + ! FV only: commented out for SIMA + !call polar_average(loc_arr) + end subroutine read_2d_trc + +!------------------------------------------------------------------------ + + ! Read zonal average data + subroutine read_za_trc(fid, vid, loc_arr, strt, cnt, file, order) + use physics_grid, only: pcols => columns_on_task + + use physics_grid, only: get_rlat_all_p + use cam_abortutils, only: check_allocate + + use pio, only: pio_get_var + + use interpolate_data, only: lininterp_init, lininterp, interp_type, lininterp_finish + + type(file_desc_t), intent(in) :: fid + type(var_desc_t), intent(in) :: vid + integer, intent(in) :: strt(:), cnt(:) + integer, intent(in) :: order(2) + real(r8), intent(out):: loc_arr(:, :) + type(trfile), intent(in) :: file + + type(interp_type) :: lat_wgts + real(r8) :: to_lats(pcols), wrk(pcols) + real(r8), allocatable, target :: wrk2d(:, :) + real(r8), pointer :: wrk2d_in(:, :) + integer :: k, ierr + character(len=512) :: errmsg + character(len=*), parameter :: subname = "read_za_trc" + + nullify (wrk2d_in) + allocate (wrk2d(cnt(1), cnt(2)), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'wrk2d(cnt(1), cnt(2))', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + if (order(1) /= 1 .or. order(2) /= 2 .or. cnt(1) /= file%nlat .or. cnt(2) /= file%nlev) then + allocate (wrk2d_in(file%nlat, file%nlev), stat=ierr, errmsg=errmsg) + + call check_allocate(ierr, subname, 'wrk2d_in(file%nlat, file%nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + end if + + ierr = pio_get_var(fid, vid, strt, cnt, wrk2d) + if (associated(wrk2d_in)) then + wrk2d_in = reshape(wrk2d(:, :), (/file%nlat, file%nlev/), order=order) + deallocate (wrk2d) + else + wrk2d_in => wrk2d + end if + + call get_rlat_all_p(pcols, to_lats) + + call lininterp_init(file%lats, file%nlat, to_lats, pcols, 1, lat_wgts) + do k = 1, file%nlev + call lininterp(wrk2d_in(:, k), file%nlat, wrk(1:pcols), pcols, lat_wgts) + loc_arr(1:pcols, k) = wrk(1:pcols) + end do + call lininterp_finish(lat_wgts) + + if (allocated(wrk2d)) then + deallocate (wrk2d) + else + deallocate (wrk2d_in) + end if + end subroutine read_za_trc + + ! this assumes the input data is gridded to match the physics grid + subroutine read_physgrid_2d(ncid, varname, recno, data) + use physics_grid, only: pcols => columns_on_task + + use cam_field_read, only: cam_read_field + + use cam_abortutils, only: endrun + + type(file_desc_t), intent(inout) :: ncid + character(len=*), intent(in) :: varname + integer, intent(in) :: recno + real(r8), intent(out) :: data(1:pcols) + + logical :: found + + call cam_read_field(varname=varname, ncid=ncid, & + field=data(:pcols), readvar=found, & + gridname='physgrid', & + timelevel=recno) + + if (.not. found) then + call endrun('tracer_data::read_physgrid_2d: Could not find '//trim(varname)//' field in input datafile') + end if + + end subroutine read_physgrid_2d + + ! this assumes the input data is gridded to match the physics grid + subroutine read_physgrid_3d(ncid, varname, vrt_coord_name, nlevs, recno, data) + use physics_grid, only: pcols => columns_on_task + + use cam_field_read, only: cam_read_field + + use cam_abortutils, only: endrun + + type(file_desc_t), intent(inout) :: ncid + character(len=*), intent(in) :: varname + character(len=*), intent(in) :: vrt_coord_name + integer, intent(in) :: nlevs + integer, intent(in) :: recno + real(r8), intent(out) :: data(1:pcols, 1:nlevs) + + logical :: found + + call cam_read_field(varname=varname, ncid=ncid, & + field=data(:pcols,:nlevs), readvar=found, & + gridname='physgrid', & + timelevel=recno, & + dim3name=vrt_coord_name, dim3_bnds=[1, nlevs]) + + if (.not. found) then + call endrun('tracer_data::read_physgrid_3d: Could not find '//trim(varname)//' field in input datafile') + end if + + end subroutine read_physgrid_3d + + !------------------------------------------------------------------------ + + subroutine read_3d_trc(fid, vid, loc_arr, strt, cnt, file, order) + use physics_grid, only: pcols => columns_on_task + use physics_grid, only: get_rlat_all_p, get_rlon_all_p + use physics_grid, only: plon => hdim1_d, plat => hdim2_d + + use physconst, only: pi + + use cam_abortutils, only: check_allocate, endrun + use cam_logfile, only: iulog + use perf_mod, only: t_startf, t_stopf + + use pio, only: pio_get_var + + ! Interpolation utils + use interpolate_data, only: lininterp_init, lininterp, interp_type, lininterp_finish + use horizontal_interpolate, only: xy_interp + + type(file_desc_t), intent(in) :: fid + type(var_desc_t), intent(in) :: vid + real(r8), intent(out) :: loc_arr(:, :) + integer, intent(in) :: strt(:) + integer, intent(in) :: cnt(:) + type(trfile), intent(in) :: file + integer, intent(in) :: order(3) + + integer :: lons(pcols), lats(pcols) + integer :: ierr + + real(r8), allocatable, target :: wrk3d(:, :, :) + real(r8), pointer :: wrk3d_in(:, :, :) + real(r8) :: to_lons(pcols), to_lats(pcols) + real(r8), parameter :: zero = 0_r8, twopi = 2_r8*pi + type(interp_type) :: lon_wgts, lat_wgts + + character(len=512) :: errmsg + character(len=*), parameter :: subname = "read_3d_trc" + + loc_arr(:, :) = 0._r8 + nullify (wrk3d_in) + allocate (wrk3d(cnt(1), cnt(2), cnt(3)), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'wrk3d(cnt(1), cnt(2), cnt(3))', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + ierr = pio_get_var(fid, vid, strt, cnt, wrk3d) + + if (order(1) /= 1 .or. order(2) /= 2 .or. order(3) /= 3 .or. & + cnt(1) /= file%nlon .or. cnt(2) /= file%nlat .or. cnt(3) /= file%nlev) then + allocate (wrk3d_in(file%nlon, file%nlat, file%nlev), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'wrk3d_in(file%nlon, file%nlat, file%nlev)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + wrk3d_in = reshape(wrk3d(:, :, :), (/file%nlon, file%nlat, file%nlev/), order=order) + deallocate (wrk3d) + else + wrk3d_in => wrk3d + end if + + ! If weighting by latitude, then perform horizontal interpolation by using weight_x, weight_y + if (file%weight_by_lat) then + + call t_startf('xy_interp') + if (file%dist) then + lons(:pcols) = lon_global_grid_ndx(:pcols) + lats(:pcols) = lat_global_grid_ndx(:pcols) + + call xy_interp(file%nlon, file%nlat, file%nlev, plon, plat, pcols, pcols, & + file%weight0_x, file%weight0_y, wrk3d_in, loc_arr(:, :), & + lons, lats, file%count0_x, file%count0_y, file%index0_x, file%index0_y) + else + lons(:pcols) = lon_global_grid_ndx(:pcols) + lats(:pcols) = lat_global_grid_ndx(:pcols) + + call xy_interp(file%nlon, file%nlat, file%nlev, plon, plat, pcols, pcols, & + file%weight_x, file%weight_y, wrk3d_in, loc_arr(:, :), & + lons, lats, file%count_x, file%count_y, file%index_x, file%index_y) + end if + call t_stopf('xy_interp') + + else + call get_rlat_all_p(pcols, to_lats) + call get_rlon_all_p(pcols, to_lons) + + call lininterp_init(file%lons, file%nlon, to_lons(1:pcols), pcols, 2, lon_wgts, zero, twopi) + call lininterp_init(file%lats, file%nlat, to_lats(1:pcols), pcols, 1, lat_wgts) + + call lininterp(wrk3d_in, file%nlon, file%nlat, file%nlev, loc_arr(:, :), pcols, pcols, lon_wgts, lat_wgts) + + call lininterp_finish(lon_wgts) + call lininterp_finish(lat_wgts) + end if + + if (allocated(wrk3d)) then + deallocate(wrk3d, stat=ierr) + else + deallocate(wrk3d_in, stat=ierr) + end if + if (ierr /= 0) then + write(iulog, *) subname//': failed to deallocate wrk3d array; error = ', ierr + call endrun(subname//': failed to deallocate wrk3d array') + end if + + ! FV only: commented out for SIMA + !call polar_average(file%nlev, loc_arr) + end subroutine read_3d_trc + +!------------------------------------------------------------------------------ + + subroutine interpolate_trcdata( & + ncol, pver, pverp, & + pmid, pint, phis, zi, & + flds, file) + + use ccpp_kinds, only: kind_phys + use physconst, only: cday, rga + + ! dependency on atmos_phys to_be_ccppized code: + use ccpp_tuvx_utils, only: rebin + + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: pverp + ! state variables used for interpolation + real(kind_phys), intent(in) :: pmid(:, :) + real(kind_phys), intent(in) :: pint(:, :) + real(kind_phys), intent(in) :: phis(:) + real(kind_phys), intent(in) :: zi(:, :) + + type(trfld), intent(inout) :: flds(:) + type(trfile), intent(inout) :: file + + real(r8) :: fact1, fact2 + real(r8) :: deltat + integer :: f, nflds, i, k + real(r8) :: ps(ncol) + real(r8) :: datain(ncol, file%nlev) + real(r8) :: pin(ncol, file%nlev) + real(r8) :: model_z(pverp) + real(r8) :: data_col(pver) + + real(r8), parameter :: m2km = 1.e-3_r8 + + nflds = size(flds) + + if (file%interp_recs == 4) then + deltat = file%datatimes(3) - file%datatimes(1) + fact1 = (file%datatimes(3) - file%datatimem)/deltat + fact2 = 1._r8 - fact1 + + if (file%has_ps) then + file%ps_in(1)%data(:ncol) = fact1*file%ps_in(1)%data(:ncol) + fact2*file%ps_in(3)%data(:ncol) + end if + do f = 1, nflds + flds(f)%input(1)%data(:ncol, :) = fact1*flds(f)%input(1)%data(:ncol, :) + fact2*flds(f)%input(3)%data(:ncol, :) + end do + + deltat = file%datatimes(4) - file%datatimes(2) + fact1 = (file%datatimes(4) - file%datatimep)/deltat + fact2 = 1._r8 - fact1 + + if (file%has_ps) then + file%ps_in(2)%data(:ncol) = fact1*file%ps_in(2)%data(:ncol) + fact2*file%ps_in(4)%data(:ncol) + end if + do f = 1, nflds + flds(f)%input(2)%data(:ncol, :) = fact1*flds(f)%input(2)%data(:ncol, :) + fact2*flds(f)%input(4)%data(:ncol, :) + end do + + end if + !------------------------------------------------------------------------- + ! If file%interp_recs=1 then no time interpolation -- set + ! fact1=1 and fact2=0 and will just use first value unmodified + !------------------------------------------------------------------------- + + if (file%interp_recs == 1) then + fact1 = 1._r8 + fact2 = 0._r8 + else + file%interp_recs = 2 + + deltat = file%datatimep - file%datatimem + + if (file%cyclical .and. (deltat < 0._r8)) then + deltat = deltat + file%one_yr + if (file%datatimep >= file%curr_mod_time) then + fact1 = (file%datatimep - file%curr_mod_time)/deltat + else + fact1 = (file%datatimep + file%one_yr - file%curr_mod_time)/deltat + end if + else + fact1 = (file%datatimep - file%curr_mod_time)/deltat + end if + + ! this assures that FIXED data are b4b on restarts + if (file%fixed) then + fact1 = dble(int(fact1*cday + .5_r8))/dble(cday) + end if + fact2 = 1._r8 - fact1 + end if + + fld_loop: do f = 1, nflds + + if (file%alt_data) then + if (fact2 == 0) then ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps) + datain(:ncol, :) = fact1*flds(f)%input(nm)%data(:ncol, :) + else + datain(:ncol, :) = fact1*flds(f)%input(nm)%data(:ncol, :) + fact2*flds(f)%input(np)%data(:ncol, :) + end if + do i = 1, ncol + model_z(1:pverp) = m2km*zi(i, pverp:1:-1) + if (file%geop_alt) then + model_z(1:pverp) = model_z(1:pverp) + m2km*phis(i)*rga + end if + if (file%conserve_column) then + call interpz_conserve(file%nlev, pver, file%ilevs, model_z, datain(i, :), data_col(:)) + else + call rebin(file%nlev, pver, file%ilevs, model_z, datain(i, :), data_col(:)) + end if + flds(f)%data(i, :) = data_col(pver:1:-1) + end do + else ! .not. alt_data + if (file%nlev > 1) then + if (file%has_ps) then + if (fact2 == 0) then ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps) + ps(:ncol) = fact1*file%ps_in(nm)%data(:ncol) + else + ps(:ncol) = fact1*file%ps_in(nm)%data(:ncol) + fact2*file%ps_in(np)%data(:ncol) + end if + do i = 1, ncol + do k = 1, file%nlev + pin(i, k) = file%p0*file%hyam(k) + ps(i)*file%hybm(k) + end do + end do + else + do k = 1, file%nlev + pin(:, k) = file%levs(k) + end do + end if + end if + + if (flds(f)%srf_fld) then + do i = 1, ncol + if (fact2 == 0) then ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps) + flds(f)%data(i, 1) = & + fact1*flds(f)%input(nm)%data(i, 1) + else + flds(f)%data(i, 1) = & + fact1*flds(f)%input(nm)%data(i, 1) + fact2*flds(f)%input(np)%data(i, 1) + end if + end do + else + if (fact2 == 0) then ! This needed as %data is not set if fact2=0 (and lahey compiler core dumps) + datain(:ncol, :) = fact1*flds(f)%input(nm)%data(:ncol, :) + else + datain(:ncol, :) = fact1*flds(f)%input(nm)%data(:ncol, :) + fact2*flds(f)%input(np)%data(:ncol, :) + end if + if (file%top_bndry) then + call vert_interp_ub(ncol, file%nlev, file%levs, datain(:ncol, :), flds(f)%data(:ncol, 1)) + else if (file%top_layer) then + call vert_interp_ub_var(ncol, file%nlev, file%levs, pmid(:ncol, 1), datain(:ncol, :), flds(f)%data(:ncol, 1)) + else if (file%conserve_column) then + call vert_interp_mixrat(ncol, file%nlev, pver, pint, & + datain, flds(f)%data(:, :), & + file%p0, ps, file%hyai, file%hybi, file%dist) + else + call vert_interp(ncol, file%nlev, pin, pmid, datain, flds(f)%data(:, :)) + end if + end if + + end if + + end do fld_loop + + end subroutine interpolate_trcdata + + subroutine get_dimension(fid, dname, dsize, dimid, data) + use cam_abortutils, only: check_allocate, endrun + use cam_logfile, only: iulog + use shr_kind_mod, only: shr_kind_cm + + use pio, only: pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR + use pio, only: pio_inq_dimid, pio_inq_dimlen, pio_inq_varid + use pio, only: pio_get_var + + type(file_desc_t), intent(inout) :: fid + character(*), intent(in) :: dname + integer, intent(out) :: dsize + + integer, optional, intent(out) :: dimid + real(r8), optional, allocatable, dimension(:) :: data + + integer :: vid, ierr, id + integer :: err_handling + + character(len=shr_kind_cm) :: errmsg + + character(len=*), parameter :: sub = 'get_dimension' + + call pio_seterrorhandling(fid, PIO_BCAST_ERROR, oldmethod=err_handling) + ierr = pio_inq_dimid(fid, dname, id) + call pio_seterrorhandling(fid, err_handling) + + if (ierr == PIO_NOERR) then + ierr = pio_inq_dimlen(fid, id, dsize) + + if (present(dimid)) then + dimid = id + end if + + if (present(data)) then + if (allocated(data)) then + deallocate (data, stat=ierr) + if (ierr /= 0) then + write (iulog, *) sub//': data deallocation error = ', ierr + call endrun(sub//': failed to deallocate data array') + end if + end if + allocate (data(dsize), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, sub, 'data(dsize)', file=__FILE__, line=__LINE__, errmsg=errmsg) + + ierr = pio_inq_varid(fid, dname, vid) + ierr = pio_get_var(fid, vid, data) + end if + else + dsize = 1 + if (present(dimid)) then + dimid = -1 + end if + end if + + end subroutine get_dimension + + subroutine set_cycle_indices(fileid, cyc_ndx_beg, cyc_ndx_end, cyc_yr) + use cam_abortutils, only: check_allocate, endrun + + use pio, only: pio_get_var, pio_inq_varid + + type(file_desc_t), intent(inout) :: fileid + integer, intent(out) :: cyc_ndx_beg + integer, intent(out) :: cyc_ndx_end + integer, intent(in) :: cyc_yr + + character(len=512) :: errmsg + character(len=*), parameter :: subname = "set_cycle_indices" + + integer, allocatable, dimension(:) :: dates + integer :: timesize, i, errflg, year, ierr + type(var_desc_t) :: dateid + call get_dimension(fileid, 'time', timesize) + cyc_ndx_beg = -1 + + allocate (dates(timesize), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'dates(timesize)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + ierr = pio_inq_varid(fileid, 'date', dateid) + ierr = pio_get_var(fileid, dateid, dates) + + do i = 1, timesize + year = dates(i)/10000 + if (year == cyc_yr) then + if (cyc_ndx_beg < 0) then + cyc_ndx_beg = i + end if + cyc_ndx_end = i + end if + end do + deallocate (dates, stat=errflg) + if (errflg /= 0) then + call endrun(subname // ': failed to deallocate dates array') + end if + if (cyc_ndx_beg < 0) then + write (*, *) subname // ': cycle year not found : ', cyc_yr + call endrun(subname // ': cycle year not found') + end if + + end subroutine set_cycle_indices + + subroutine open_trc_datafile(fname, path, piofile, times, cyc_ndx_beg, cyc_ndx_end, cyc_yr) + use ioFileMod, only: cam_get_file + use cam_pio_utils, only: cam_pio_openfile + use cam_abortutils, only: check_allocate, endrun + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + + use time_manager, only: set_time_float_from_date, set_date_from_time_float + + use pio, only: PIO_NOWRITE + use pio, only: pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR + use pio, only: pio_inq_varid + use pio, only: pio_get_var + + character(*), intent(in) :: fname + character(*), intent(in) :: path + type(file_desc_t), intent(inout) :: piofile + real(r8), allocatable, intent(inout) :: times(:) + + integer, optional, intent(out) :: cyc_ndx_beg + integer, optional, intent(out) :: cyc_ndx_end + integer, optional, intent(in) :: cyc_yr + + character(len=shr_kind_cl) :: filen, filepath + integer :: year, month, day, i, timesize + integer :: dateid, secid + integer, allocatable, dimension(:) :: dates, datesecs + integer :: ierr + logical :: need_first_ndx + integer :: err_handling + + character(len=512) :: errmsg + character(len=*), parameter :: subname = "open_trc_datafile" + + if (len_trim(path) == 0) then + filepath = trim(fname) + else + filepath = trim(path)//'/'//trim(fname) + end if + ! + ! open file and get fileid + ! + call cam_get_file(filepath, filen, allow_fail=.false.) + call cam_pio_openfile(piofile, filen, PIO_NOWRITE) + if (masterproc) write (iulog, *) 'open_trc_datafile: ', trim(filen) + + call get_dimension(piofile, 'time', timesize) + + if (allocated(times)) then + deallocate (times, stat=ierr) + if (ierr /= 0) then + write (iulog, *) 'open_trc_datafile: data deallocation error = ', ierr + call endrun('open_trc_datafile: failed to deallocate data array') + end if + end if + + allocate (times(timesize), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'times(timesize)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + allocate (dates(timesize), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'dates(timesize)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + allocate (datesecs(timesize), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'datesecs(timesize)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + ierr = pio_inq_varid(piofile, 'date', dateid) + call pio_seterrorhandling(piofile, PIO_BCAST_ERROR, oldmethod=err_handling) + ierr = pio_inq_varid(piofile, 'datesec', secid) + call pio_seterrorhandling(piofile, err_handling) + + if (ierr == PIO_NOERR) then + ierr = pio_get_var(piofile, secid, datesecs) + else + datesecs = 0 + end if + + ierr = pio_get_var(piofile, dateid, dates) + need_first_ndx = .true. + + do i = 1, timesize + year = dates(i)/10000 + month = mod(dates(i), 10000)/100 + day = mod(dates(i), 100) + call set_time_float_from_date(times(i), year, month, day, datesecs(i)) + if (present(cyc_yr)) then + if (year == cyc_yr) then + if (present(cyc_ndx_beg) .and. need_first_ndx) then + cyc_ndx_beg = i + need_first_ndx = .false. + end if + if (present(cyc_ndx_end)) then + cyc_ndx_end = i + end if + end if + end if + end do + + deallocate (dates, stat=ierr) + if (ierr /= 0) then + if (masterproc) write (iulog, *) subname //': failed to deallocate dates array; error = ', ierr + call endrun(subname //': failed to deallocate dates array') + end if + deallocate (datesecs, stat=ierr) + if (ierr /= 0) then + if (masterproc) write (iulog, *) subname //': failed to deallocate datesec array; error = ', ierr + call endrun(subname //': failed to deallocate datesec array') + end if + + if (present(cyc_yr) .and. present(cyc_ndx_beg)) then + if (cyc_ndx_beg < 0) then + write (iulog, *) subname //': cycle year not found : ', cyc_yr + call endrun(subname //': cycle year not found '//trim(filepath)) + end if + end if + + end subroutine open_trc_datafile + + subroutine specify_fields(specifier, fields) + use cam_abortutils, only: check_allocate + use shr_kind_mod, only: shr_kind_cm + + character(len=*), intent(in) :: specifier(:) + type(trfld), pointer, dimension(:) :: fields + + integer :: fld_cnt, astat + integer :: i, j + character(len=shr_kind_cl) :: str1, str2 + character(len=32), allocatable, dimension(:) :: fld_name, src_name + integer :: nflds + + character(len=shr_kind_cm) :: errmsg + character(len=*), parameter :: sub = 'specify_fields' + + nflds = size(specifier) + + allocate (fld_name(nflds), src_name(nflds), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'fld_name(nflds), src_name(nflds)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + fld_cnt = 0 + + count_cnst: do i = 1, nflds + + if (len_trim(specifier(i)) == 0) then + exit count_cnst + end if + + j = scan(specifier(i), ':') + + if (j > 0) then + str1 = trim(adjustl(specifier(i) (:j - 1))) + str2 = trim(adjustl(specifier(i) (j + 1:))) + fld_name(i) = trim(adjustl(str1)) + src_name(i) = trim(adjustl(str2)) + else + fld_name(i) = trim(adjustl(specifier(i))) + src_name(i) = trim(adjustl(specifier(i))) + end if + + fld_cnt = fld_cnt + 1 + + end do count_cnst + + if (fld_cnt < 1) then + nullify (fields) + return + end if + + !----------------------------------------------------------------------- + ! ... allocate field type array + !----------------------------------------------------------------------- + allocate (fields(fld_cnt), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'fields(fld_cnt)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + do i = 1, fld_cnt + fields(i)%fldnam = fld_name(i) + fields(i)%srcnam = src_name(i) + end do + + deallocate (fld_name, src_name) + + end subroutine specify_fields + + ! This routine advances to the next file + subroutine advance_file(file) + use cam_abortutils, only: check_allocate, endrun + use cam_logfile, only: iulog + use shr_kind_mod, only: shr_kind_cm + + use pio, only: pio_closefile + + type(trfile), intent(inout) :: file + + !----------------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------------- + character(len=shr_kind_cl) :: ctmp + character(len=shr_kind_cl) :: loc_fname + integer :: astat + + character(len=shr_kind_cm) :: errmsg + character(len=*), parameter :: sub = 'advance_file' + + !----------------------------------------------------------------------- + ! close current file ... + !----------------------------------------------------------------------- + call pio_closefile(file%curr_fileid) + + !----------------------------------------------------------------------- + ! Advance the filename and file id + !----------------------------------------------------------------------- + file%curr_filename = file%next_filename + file%curr_fileid = file%next_fileid + + !----------------------------------------------------------------------- + ! Advance the curr_data_times + !----------------------------------------------------------------------- + deallocate (file%curr_data_times, stat=astat) + if (astat /= 0) then + write (iulog, *) sub//': failed to deallocate file%curr_data_times array; error = ', astat + call endrun(sub//': failed to deallocate file%curr_data_times array') + end if + + allocate (file%curr_data_times(size(file%next_data_times)), stat=astat, errmsg=errmsg) + call check_allocate(astat, sub, 'file%curr_data_times(size(file%next_data_times))', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + file%curr_data_times(:) = file%next_data_times(:) + + !----------------------------------------------------------------------- + ! delete information about next file (as was just assigned to current) + !----------------------------------------------------------------------- + file%next_filename = '' + + deallocate (file%next_data_times, stat=astat) + if (astat /= 0) then + write (iulog, *) sub//': failed to deallocate file%next_data_times array; error = ', astat + call endrun(sub//': failed to deallocate file%next_data_times array') + end if + + end subroutine advance_file + +!------------------------------------------------------------------------------ + + subroutine init_trc_restart(whence, piofile, tr_file) + use pio, only: pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR + use pio, only: pio_char + use pio, only: pio_inq_dimid, pio_def_dim + use pio, only: pio_put_att, pio_def_var + + character(len=*), intent(in) :: whence + type(file_desc_t), intent(inout) :: piofile + type(trfile), intent(inout) :: tr_file + + character(len=32) :: name + integer :: ioerr, mcdimid, maxlen + integer :: err_handling + + ! Dimension should already be defined in restart file + call pio_seterrorhandling(pioFile, PIO_BCAST_ERROR, oldmethod=err_handling) + ioerr = pio_inq_dimid(pioFile, 'max_chars', mcdimid) + call pio_seterrorhandling(pioFile, err_handling) + ! but define it if nessasary + if (ioerr /= PIO_NOERR) then + ioerr = pio_def_dim(pioFile, 'max_chars', SHR_KIND_CL, mcdimid) + end if + + if (len_trim(tr_file%curr_filename) > 1) then + allocate (tr_file%currfnameid) + name = trim(whence)//'_curr_fname' + ioerr = pio_def_var(pioFile, name, pio_char, (/mcdimid/), tr_file%currfnameid) + ioerr = pio_put_att(pioFile, tr_file%currfnameid, 'offset_time', tr_file%offset_time) + maxlen = len_trim(tr_file%curr_filename) + ioerr = pio_put_att(pioFile, tr_file%currfnameid, 'actual_len', maxlen) + else + deallocate(tr_file%currfnameid) + end if + + if (len_trim(tr_file%next_filename) > 1) then + allocate (tr_file%nextfnameid) + name = trim(whence)//'_next_fname' + ioerr = pio_def_var(pioFile, name, pio_char, (/mcdimid/), tr_file%nextfnameid) + maxlen = len_trim(tr_file%next_filename) + ioerr = pio_put_att(pioFile, tr_file%nextfnameid, 'actual_len', maxlen) + else + deallocate(tr_file%nextfnameid) + end if + end subroutine init_trc_restart + + ! writes file names to restart file + subroutine write_trc_restart(piofile, tr_file) + use pio, only: pio_put_var + + type(file_desc_t), intent(inout) :: piofile + type(trfile), intent(inout) :: tr_file + + integer :: ioerr ! error status + if (allocated(tr_file%currfnameid)) then + ioerr = pio_put_var(pioFile, tr_file%currfnameid, tr_file%curr_filename) + deallocate (tr_file%currfnameid) + end if + if (allocated(tr_file%nextfnameid)) then + ioerr = pio_put_var(pioFile, tr_file%nextfnameid, tr_file%next_filename) + deallocate (tr_file%nextfnameid) + end if + + end subroutine write_trc_restart + + ! reads file names from restart file + subroutine read_trc_restart(whence, piofile, tr_file) + use pio, only: pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR + use pio, only: pio_inq_varid, pio_get_att, pio_get_var + + character(len=*), intent(in) :: whence + type(file_desc_t), intent(inout) :: piofile + type(trfile), intent(inout) :: tr_file + type(var_desc_t) :: vdesc + character(len=64) :: name + integer :: ioerr ! error status + integer :: slen + integer :: err_handling + + call pio_seterrorhandling(piofile, PIO_BCAST_ERROR, oldmethod=err_handling) + name = trim(whence)//'_curr_fname' + ioerr = pio_inq_varid(piofile, name, vdesc) + if (ioerr == PIO_NOERR) then + tr_file%curr_filename = ' ' + ioerr = pio_get_att(piofile, vdesc, 'offset_time', tr_file%offset_time) + ioerr = pio_get_att(piofile, vdesc, 'actual_len', slen) + ioerr = pio_get_var(piofile, vdesc, tr_file%curr_filename) + if (slen < SHR_KIND_CL) tr_file%curr_filename(slen + 1:) = ' ' + end if + + name = trim(whence)//'_next_fname' + ioerr = pio_inq_varid(piofile, name, vdesc) + if (ioerr == PIO_NOERR) then + tr_file%next_filename = ' ' + ioerr = pio_get_att(piofile, vdesc, 'actual_len', slen) + ioerr = pio_get_var(piofile, vdesc, tr_file%next_filename) + if (slen < SHR_KIND_CL) tr_file%next_filename(slen + 1:) = ' ' + end if + call pio_seterrorhandling(piofile, err_handling) + + end subroutine read_trc_restart + + !------------------------------------------------------------------------------ + ! Various utility subroutines below: + !------------------------------------------------------------------------------ + pure subroutine interpz_conserve(nsrc, ntrg, src_x, trg_x, src, trg) + + integer, intent(in) :: nsrc ! dimension source array + integer, intent(in) :: ntrg ! dimension target array + real(r8), intent(in) :: src_x(nsrc + 1) ! source coordinates + real(r8), intent(in) :: trg_x(ntrg + 1) ! target coordinates + real(r8), intent(in) :: src(nsrc) ! source array + real(r8), intent(out) :: trg(ntrg) ! target array + + !--------------------------------------------------------------- + ! ... local variables + !--------------------------------------------------------------- + integer :: i, j + integer :: sil + real(r8) :: tl, y + real(r8) :: bot, top + + do i = 1, ntrg + tl = trg_x(i) + if ((tl < src_x(nsrc + 1)) .and. (trg_x(i + 1) > src_x(1))) then + do sil = 1, nsrc + if ((tl - src_x(sil))*(tl - src_x(sil + 1)) <= 0.0_r8) then + exit + end if + end do + + if (tl < src_x(1)) sil = 1 + + y = 0.0_r8 + bot = max(tl, src_x(1)) + top = trg_x(i + 1) + do j = sil, nsrc + if (top > src_x(j + 1)) then + y = y + (src_x(j + 1) - bot)*src(j)/(src_x(j + 1) - src_x(j)) + bot = src_x(j + 1) + else + y = y + (top - bot)*src(j)/(src_x(j + 1) - src_x(j)) + exit + end if + end do + trg(i) = y + else + trg(i) = 0.0_r8 + end if + end do + + if (trg_x(1) > src_x(1)) then + top = trg_x(1) + bot = src_x(1) + y = 0.0_r8 + do j = 1, nsrc + if (top > src_x(j + 1)) then + y = y + (src_x(j + 1) - bot)*src(j)/(src_x(j + 1) - src_x(j)) + bot = src_x(j + 1) + else + y = y + (top - bot)*src(j)/(src_x(j + 1) - src_x(j)) + exit + end if + end do + trg(1) = trg(1) + y + end if + + end subroutine interpz_conserve + + subroutine vert_interp_mixrat(ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi, use_flight_distance) + use cam_abortutils, only: endrun + + integer, intent(in) :: ncol + integer, intent(in) :: nsrc ! dimension source array + integer, intent(in) :: ntrg ! dimension target array + real(r8), intent(in) :: trg_x(ncol, ntrg + 1) ! target coordinates + real(r8), intent(in) :: src(ncol, nsrc) ! source array + real(r8), intent(out) :: trg(ncol, ntrg) ! target array + real(r8), intent(in) :: ps(ncol) ! surface pressure + real(r8), intent(in) :: p0 + real(r8), intent(in) :: hyai(nsrc + 1) + real(r8), intent(in) :: hybi(nsrc + 1) + logical, intent(in) :: use_flight_distance ! .true. = flight distance, .false. = mixing ratio + + !--------------------------------------------------------------- + ! ... local variables + !--------------------------------------------------------------- + integer :: i, j, n + real(r8) :: y, trg_lo, trg_hi, src_lo, src_hi, overlap, outside + real(r8) :: src_x(nsrc + 1) ! source coordinates + + character(len=*), parameter :: sub = 'vert_interp_mixrat' + + do n = 1, ncol ! loop over columns + + trg(n, :) = 0.0_r8 ! probably not needed + + ! calculate source pressure levels from hybrid coords + do i = 1, nsrc + 1 + src_x(i) = p0*hyai(i) + ps(n)*hybi(i) + end do + + ! Check the invariant that src_x and trg_x values are + ! ascending. This could also be checked at an earlier stage to + ! avoid doing so for every interpolation call. + if (.not. ALL(src_x(1:nsrc) < src_x(2:nsrc + 1))) then + call endrun(sub//': src_x values are not ascending') + end if + if (.not. ALL(trg_x(n, 1:ntrg) < trg_x(n, 2:ntrg + 1))) then + call endrun(sub//': trg_x values are not ascending') + end if + + do i = 1, ntrg + y = 0.0_r8 + + trg_lo = trg_x(n, i) + trg_hi = trg_x(n, i + 1) + + do j = 1, nsrc + src_lo = src_x(j) + src_hi = src_x(j + 1) + + overlap = min(src_hi, trg_hi) - max(src_lo, trg_lo) + if (overlap > 0.0_r8) then + if (use_flight_distance) then + ! add input based on the overlap fraction + y = y + src(n, j)*overlap/(src_hi - src_lo) + else + ! convert to mass by multiplying by dp + y = y + src(n, j)*overlap + end if + end if + end do + trg(n, i) = y + end do + + ! Handle source values outside the target range. Since we want + ! to preserve the total amount, add these to the first/last + ! target bucket. + trg_lo = trg_x(n, 1) + y = 0.0_r8 + do j = 1, nsrc + src_lo = src_x(j) + src_hi = src_x(j + 1) + + if (src_lo < trg_lo) then + if (src_hi <= trg_lo) then + ! whole source interval is outside the target range + outside = src_hi - src_lo + else + ! There was some overlap, which would have been added + ! previously. Only add the parts outside the target + ! range. + outside = trg_lo - src_lo + end if + if (use_flight_distance) then + ! add the input scaled by the fraction outside + y = y + src(n, j)*outside/(src_hi - src_lo) + else + ! convert to mass by multiplying by dp + y = y + src(n, j)*outside + end if + else + exit + end if + end do + trg(n, 1) = trg(n, 1) + y + + trg_hi = trg_x(n, ntrg + 1) + y = 0.0_r8 + do j = nsrc, 1, -1 + src_lo = src_x(j) + src_hi = src_x(j + 1) + + if (src_hi > trg_hi) then + if (src_lo >= trg_hi) then + ! whole source interval is outside the target range + outside = src_hi - src_lo + else + ! There was some overlap, which would have been added + ! previously. Only add the parts outside the target + ! range. + outside = src_hi - trg_hi + end if + if (use_flight_distance) then + ! add the full input + y = y + src(n, j)*outside/(src_hi - src_lo) + else + ! convert to mass by multiplying by dp + y = y + src(n, j)*outside + end if + else + exit + end if + end do + trg(n, ntrg) = trg(n, ntrg) + y + + ! turn mass into mixing ratio + if (.not. use_flight_distance) then + do i = 1, ntrg + trg(n, i) = trg(n, i)/(trg_x(n, i + 1) - trg_x(n, i)) + end do + end if + + end do + + end subroutine vert_interp_mixrat + + ! Interpolate data from current time-interpolated values to model levels + pure subroutine vert_interp(ncol, levsiz, pin, pmid, datain, dataout) + use vert_coord, only: pver + + integer, intent(in) :: ncol + integer, intent(in) :: levsiz + real(r8), intent(in) :: pin(ncol, levsiz) + real(r8), intent(in) :: pmid(ncol, pver) + real(r8), intent(in) :: datain(ncol, levsiz) + real(r8), intent(out) :: dataout(ncol, pver) + + ! local storage + integer :: i ! longitude index + integer :: k, kk, kkstart ! level indices + integer :: kupper(ncol) ! Level indices for interpolation + real(r8) :: dpu ! upper level pressure difference + real(r8) :: dpl ! lower level pressure difference + + !-------------------------------------------------------------------------- + ! + ! Initialize index array + ! + do i = 1, ncol + kupper(i) = 1 + end do + + do k = 1, pver + ! + ! Top level we need to start looking is the top level for the previous k + ! for all column points + ! + kkstart = levsiz + do i = 1, ncol + kkstart = min0(kkstart, kupper(i)) + end do + ! + ! Store level indices for interpolation + ! + do kk = kkstart, levsiz - 1 + do i = 1, ncol + if (pin(i, kk) < pmid(i, k) .and. pmid(i, k) <= pin(i, kk + 1)) then + kupper(i) = kk + end if + end do + end do + ! interpolate or extrapolate... + do i = 1, ncol + if (pmid(i, k) < pin(i, 1)) then + dataout(i, k) = datain(i, 1)*pmid(i, k)/pin(i, 1) + else if (pmid(i, k) > pin(i, levsiz)) then + dataout(i, k) = datain(i, levsiz) + else + dpu = pmid(i, k) - pin(i, kupper(i)) + dpl = pin(i, kupper(i) + 1) - pmid(i, k) + dataout(i, k) = (datain(i, kupper(i))*dpl + & + datain(i, kupper(i) + 1)*dpu)/(dpl + dpu) + end if + end do + end do + + end subroutine vert_interp + + ! Interpolate data from current time-interpolated values to top interface pressure + pure subroutine vert_interp_ub(ncol, nlevs, plevs, datain, dataout) + use ref_pres, only: ptop_ref + + integer, intent(in) :: ncol + integer, intent(in) :: nlevs + real(r8), intent(in) :: plevs(nlevs) + real(r8), intent(in) :: datain(ncol, nlevs) + real(r8), intent(out) :: dataout(ncol) + + ! + ! local variables + ! + integer :: i, ku, kl, kk + real(r8) :: pinterp, delp + + pinterp = ptop_ref + + if (pinterp <= plevs(1)) then + kl = 1 + ku = 1 + delp = 0._r8 + else if (pinterp >= plevs(nlevs)) then + kl = nlevs + ku = nlevs + delp = 0._r8 + else + + do kk = 2, nlevs + if (pinterp <= plevs(kk)) then + ku = kk + kl = kk - 1 + delp = log(pinterp/plevs(kk))/log(plevs(kk - 1)/plevs(kk)) + exit + end if + end do + + end if + + do i = 1, ncol + dataout(i) = datain(i, kl) + delp*(datain(i, ku) - datain(i, kl)) + end do + + end subroutine vert_interp_ub + + ! Interpolate data from current time-interpolated values to press + pure subroutine vert_interp_ub_var(ncol, nlevs, plevs, press, datain, dataout) + + integer, intent(in) :: ncol + integer, intent(in) :: nlevs + real(r8), intent(in) :: plevs(nlevs) + real(r8), intent(in) :: press(ncol) + real(r8), intent(in) :: datain(ncol, nlevs) + real(r8), intent(out) :: dataout(ncol) + + ! + ! local variables + ! + integer :: i, k + integer :: ku, kl + real(r8) :: delp + + do i = 1, ncol + if (press(i) <= plevs(1)) then + kl = 1 + ku = 1 + delp = 0._r8 + else if (press(i) >= plevs(nlevs)) then + kl = nlevs + ku = nlevs + delp = 0._r8 + else + do k = 2, nlevs + if (press(i) <= plevs(k)) then + ku = k + kl = k - 1 + delp = log(press(i)/plevs(k))/log(plevs(k - 1)/plevs(k)) + exit + end if + end do + end if + + dataout(i) = datain(i, kl) + delp*(datain(i, ku) - datain(i, kl)) + end do + + end subroutine vert_interp_ub_var + + ! Purpose: "find periodic lower bound" + ! Search the input array for the lower bound of the interval that + ! contains the input value. The returned index satifies: + ! x(index) .le. xval .lt. x(index+1) + ! Assume the array represents values in one cycle of a periodic coordinate. + ! So, if xval .lt. x(1), or xval .ge. x(nx), then the index returned is nx. + ! + ! Author: B. Eaton + pure subroutine findplb(x, nx, xval, index) + integer, intent(in) :: nx ! size of x + real(r8), intent(in) :: x(nx) ! strictly increasing array + real(r8), intent(in) :: xval ! value to be searched for in x + + integer, intent(out) :: index + + ! Local variables: + integer i + !----------------------------------------------------------------------- + + if (xval .lt. x(1) .or. xval .ge. x(nx)) then + index = nx + return + end if + + do i = 2, nx + if (xval .lt. x(i)) then + index = i - 1 + return + end if + end do + + end subroutine findplb + +end module tracer_data diff --git a/test/unit/fortran/src/core_utils/test_string_core_utils.pf b/test/unit/fortran/src/core_utils/test_string_core_utils.pf index 7db43ce6..1dc68e0d 100644 --- a/test/unit/fortran/src/core_utils/test_string_core_utils.pf +++ b/test/unit/fortran/src/core_utils/test_string_core_utils.pf @@ -651,3 +651,203 @@ subroutine test_tokenize_into_tokens_separator_by_known_string_set() @assertTrue(allocated(separator)) @assertEqual(expected_separator_comma_space_separated, separator, whitespace=keep_all) end subroutine test_tokenize_into_tokens_separator_by_known_string_set + +@test +subroutine test_get_last_significant_char_empty_string() + use funit + use string_core_utils, only: get_last_significant_char + + character(len=10) :: test_string + + test_string = '' + @assertEqual(0, get_last_significant_char(test_string)) +end subroutine test_get_last_significant_char_empty_string + +@test +subroutine test_get_last_significant_char_all_blanks() + use funit + use string_core_utils, only: get_last_significant_char + + character(len=10) :: test_string + + test_string = ' ' + @assertEqual(0, get_last_significant_char(test_string)) +end subroutine test_get_last_significant_char_all_blanks + +@test +subroutine test_get_last_significant_char_trailing_spaces() + use funit + use string_core_utils, only: get_last_significant_char + + character(len=20) :: test_string + + test_string = 'Mercury ' + @assertEqual(7, get_last_significant_char(test_string)) +end subroutine test_get_last_significant_char_trailing_spaces + +@test +subroutine test_get_last_significant_char_no_trailing_spaces() + use funit + use string_core_utils, only: get_last_significant_char + + character(len=7) :: test_string + + test_string = 'Mercury' + @assertEqual(7, get_last_significant_char(test_string)) +end subroutine test_get_last_significant_char_no_trailing_spaces + +@test +subroutine test_last_non_digit_empty_string() + use funit + use string_core_utils, only: last_non_digit + + character(len=10) :: test_string + + test_string = '' + @assertEqual(0, last_non_digit(test_string)) +end subroutine test_last_non_digit_empty_string + +@test +subroutine test_last_non_digit_all_digits() + use funit + use string_core_utils, only: last_non_digit + + character(len=10) :: test_string + + test_string = '12345' + @assertEqual(0, last_non_digit(test_string)) +end subroutine test_last_non_digit_all_digits + +@test +subroutine test_last_non_digit_with_trailing_digits() + use funit + use string_core_utils, only: last_non_digit + + character(len=20) :: test_string + + test_string = 'test_file_123' + @assertEqual(10, last_non_digit(test_string)) +end subroutine test_last_non_digit_with_trailing_digits + +@test +subroutine test_last_non_digit_no_trailing_digits() + use funit + use string_core_utils, only: last_non_digit + + character(len=20) :: test_string + + test_string = 'test_file_abc' + @assertEqual(13, last_non_digit(test_string)) +end subroutine test_last_non_digit_no_trailing_digits + +@test +subroutine test_last_non_digit_single_digit() + use funit + use string_core_utils, only: last_non_digit + + character(len=10) :: test_string + + test_string = 'file5' + @assertEqual(4, last_non_digit(test_string)) +end subroutine test_last_non_digit_single_digit + +@test +subroutine test_last_non_digit_with_leading_zeros() + use funit + use string_core_utils, only: last_non_digit + + character(len=15) :: test_string + + test_string = 'snapshot_00042' + @assertEqual(9, last_non_digit(test_string)) +end subroutine test_last_non_digit_with_leading_zeros + +@test +subroutine test_increment_string_no_trailing_digits() + use funit + use string_core_utils, only: increment_string + + character(len=20) :: test_string + integer :: result + + test_string = 'test_file' + result = increment_string(test_string, 1) + + @assertEqual(-1, result) + @assertEqual('test_file', test_string) +end subroutine test_increment_string_no_trailing_digits + +@test +subroutine test_increment_string_single_digit_increment() + use funit + use string_core_utils, only: increment_string + + character(len=20) :: test_string + integer :: result + + test_string = 'file5' + result = increment_string(test_string, 1) + + @assertEqual(0, result) + @assertEqual('file6', test_string) +end subroutine test_increment_string_single_digit_increment + +@test +subroutine test_increment_string_multiple_digits() + use funit + use string_core_utils, only: increment_string + + character(len=20) :: test_string + integer :: result + + test_string = 'snapshot_042' + result = increment_string(test_string, 1) + + @assertEqual(0, result) + @assertEqual('snapshot_043', test_string) +end subroutine test_increment_string_multiple_digits + +@test +subroutine test_increment_string_with_rollover() + use funit + use string_core_utils, only: increment_string + + character(len=20) :: test_string + integer :: result + + test_string = 'file09' + result = increment_string(test_string, 1) + + @assertEqual(0, result) + @assertEqual('file10', test_string) +end subroutine test_increment_string_with_rollover + +@test +subroutine test_increment_string_with_large_increment() + use funit + use string_core_utils, only: increment_string + + character(len=20) :: test_string + integer :: result + + test_string = 'test_000' + result = increment_string(test_string, 42) + + @assertEqual(0, result) + @assertEqual('test_042', test_string) +end subroutine test_increment_string_with_large_increment + +@test +subroutine test_increment_string_positive_overflow() + use funit + use string_core_utils, only: increment_string + + character(len=20) :: test_string + integer :: result + + test_string = 'file99' + result = increment_string(test_string, 1) + + @assertEqual(-2, result) +end subroutine test_increment_string_positive_overflow +