From ab7fb1d0edc492f125b83ce1874bdcfb7846539a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 26 Jan 2025 17:20:44 -0500 Subject: [PATCH 1/9] added matlab reader for tile file (ASCII and nc4) --- .../util/shared/matlab/read_tile_file_ASCII.m | 375 ++++++++++++++++++ .../util/shared/matlab/read_tile_file_nc4.m | 42 ++ 2 files changed, 417 insertions(+) create mode 100644 GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m create mode 100644 GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m b/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m new file mode 100644 index 0000000..3dab7c3 --- /dev/null +++ b/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m @@ -0,0 +1,375 @@ +function [ tile_coord ] = read_til_file_ASCII( fname_til, fname_catchmentdef ) + +% read tile coordinates from "til" file (EASEv2 and cube-sphere) and +% put data into tile_coord structure (similar to that defined in GEOSldas) +% +% inputs: +% fname_til : Name (with full path) of ASCII tile file (see examples below). +% fname_catchmentdef : OPTIONAL, name of matching catchment.def file +% Do NOT specify if tile file is from bcs in new layout; will be +% assembled from fname_til. +% +% fname_til examples: +% '/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/v12/geometry/EASEv2_M36/EASEv2_M36_964x406.til' +% '/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/v12/geometry/CF0360x6C_DE1440xPE0720/CF0360x6C_DE1440xPE0720-Pfafstetter.til' +% +% +% IMPORTANT: +% This reader is designed to work for bcs versions NL5, v11, v12, and newer (beginning ~2022). +% It should be backward-compatible for most older bcs versions, but there is no +% expectation that it will work universally across all old and new versions. +% +% reichle, 24 Jan 2025 +% +% ------------------------------------------------------------- + +% determine fname_catchmentdef + +%%%if ~exist( fname_catchmentdef ) + + ind=strfind(fname_til,'/'); + + % verify that tile file is from bcs in new layout and naming convention + + if ~strcmp( fname_til(ind(end-2)+1:ind(end-1)-1), 'geometry' ) + + error('cannot derive fname_catchmentdef') + + end + + resolution = fname_til(ind(end-1)+1:ind(end)-1); + + fname_catchmentdef = [ fname_til(1:ind(end-2)), 'land/', resolution, '/clsm/catchment.def' ]; + +%%%end + +% ------------------------------------------------------------- +% +% read tile file + +disp(['reading from ', fname_til]) + +ifp = fopen( fname_til, 'r' ); + +% ------------------------------------------------------------- +% +% EASE grid tile space? + +isEASE = 0; + +if ~isempty(findstr('EASE',fname_til)), isEASE=1; end + +% -------------- +% +% determine number of data columns + +if isEASE + + % NOTE: very old EASE versions have N_data_col=8 + % + % between SMAP bcs v001 (through Vv3030) and SMAP bcs v003[a] (from Vv4030), + % a new (9th) column was added in the *.til file + + N_data_col = 9; + +else + + N_data_col = 12; + +end + +% ------------------------------------------------------------- +% +% Number of integers in first header line depends on bcs version: + +N_headerline1 = 4; % for bcs versions after ~2022 + +old_versions = { 'Fortuna', ... + 'Ganymed', ... + 'Heracles', ... + 'Icarus', ... + 'Jason', ... + 'DYAMOND2', ... + 'FV3-1_0', ... + 'NL3', ... + 'NL4', ... + }; + +for kk=1:length(old_versions) + + if ~isempty( strfind( fname_til, old_versions{kk} )) + + N_headerline1 = 3; + + disp([ 'detected old bcs version ', old_versions{kk} ]) + + end +end + +% make exception for select "new" versions (stored in legacy bcs directory) that would otherwise +% be classified as "old" per the loop above + +new_versions = { 'Icarus-NLv5', ... + 'DYAMOND', ... + }; + +for kk=1:length(new_versions) + + if ~isempty( strfind( fname_til, new_versions{kk} )) + + N_headerline1 = 4; + + disp([ 'detected exception for new bcs version ', new_versions{kk} ]) + + end +end + +% echo final value of N_headerline1 + +disp(['using N_headerline1 = ', num2str(N_headerline1) ]) + +% ------------------------------------------------------------- +% +% read header + +% header line 1: + +tmpdata = fscanf( ifp, '%f', N_headerline1 ); + +N_tile_tmp = tmpdata(1); + +if N_headerline1==4 + + tile_coord.N_PfafCat = tmpdata(2); + tile_coord.raster_nx = tmpdata(3); + tile_coord.raster_ny = tmpdata(4); + +elseif N_headerline1==3 + + tile_coord.N_PfafCat = NaN; + tile_coord.raster_nx = tmpdata(2); + tile_coord.raster_ny = tmpdata(3); + +end + +% -------------- +% +% header line 2: + +tile_coord.N_grids = fscanf( ifp, '%f', 1 ); + +% verify N_grids (should be 1 for EASE and 2 for non-EASE) + +if isEASE & tile_coord.N_grids~=1, error(['unexpected N_grids for EASE tile space: ', num2str(tile_coord.N_grids)]), end + +if ~isEASE & tile_coord.N_grids~=2, error(['unexpected N_grids for non-EASE tile space: ', num2str(tile_coord.N_grids)]), end + +% Deal with older EASE bcs versions having (useless) header lines for grid 2, +% despite having the correct value of N_grids=1 in header line 2. + +if isEASE & ( ~isempty(findstr('NL3',fname_til)) | ~isempty(findstr('NL4',fname_til)) ) + + tmp_n_grids = 2; + +else + + tmp_n_grids = tile_coord.N_grids; + +end + +% -------------- +% +% header lines 3-5: + +tile_coord.Grid_Name = fscanf( ifp, '%s', 1 ); +tile_coord.IM = fscanf( ifp, '%f', 1 ); +tile_coord.JM = fscanf( ifp, '%f', 1 ); + +% -------------- +% +% header lines 6-8 (if present): + +if tmp_n_grids==2 + + % NOTE: EASE NL3 and NL4 contain additional header lines for second grid, despite (correct) N_grids=1 + + tile_coord.Grid2_Name = fscanf( ifp, '%s', 1 ); + tile_coord.IM2 = fscanf( ifp, '%f', 1 ); + tile_coord.JM2 = fscanf( ifp, '%f', 1 ); + +end + +% -------------- +% +% read data + +tmpdata = fscanf( ifp, '%f' ); + +fclose(ifp); + +disp('done reading "til" file') + +% verify N_tile*N_data_col against number of data read + +if length(tmpdata)~=N_tile_tmp*N_data_col, error('something wrong with N_tile*N_data_col'), end + +% convert into 2d array + +tmpdata = transpose( reshape(tmpdata, [N_data_col, N_tile_tmp]) ); + +% -------------------------------------------------- +% +% assign tile_id +% +% by convention, tile_id is equal to index number of tile in *.til file + +tmp_tileid = transpose(1:N_tile_tmp); + +% -------------------------------------------------- +% +% subset (if requested) + +%ind = find(tmpdata(:,1)==100); % keep only land +%ind = find(tmpdata(:,1)== 19); % keep only lakes +%ind = find(tmpdata(:,1)== 20); % keep only ice +%ind = find(tmpdata(:,1)~= 0); % keep land, lake, landice + +%tmpdata = tmpdata( ind,:); + +%tmp_tileid = tmp_tileid( ind, 1); + +% -------------------------------------------------- +% +% copy data into tile_coord structure + +tile_coord.N_tile = size(tmpdata,1); % number of tiles (assign here in case of subsetting above) + +% the following are universal (for EASE and non-EASE, all tile types) + +tile_coord.tile_id = tmp_tileid; % tile ID + +tile_coord.typ = tmpdata(:, 1); % tile type + +tile_coord.com_lon = tmpdata(:, 3); % center-of-mass longitude of tile +tile_coord.com_lat = tmpdata(:, 4); % center-of-mass latitude of tile + +tile_coord.i_indg = tmpdata(:, 5); % i index of tile on global "atm" (or EASE) grid +tile_coord.j_indg = tmpdata(:, 6); % j index of tile on global "atm" (or EASE) grid +tile_coord.frac_cell = tmpdata(:, 7); % area fraction of "atm" (or EASE) grid cell + +% initialize remaining fields (to be filled below) + +tmpNaN = NaN*ones(tile_coord.N_tile,1); + +tile_coord.min_lon = tmpNaN; % min longitude of tile +tile_coord.max_lon = tmpNaN; % max longitude of tile +tile_coord.min_lat = tmpNaN; % min latitude of tile +tile_coord.max_lat = tmpNaN; % max latitude of tile + +tile_coord.elev = tmpNaN; % elevation of tile + +tile_coord.area = tmpNaN; % area of "atm" grid cell +tile_coord.pfaf = tmpNaN; % index of (hydrological) Pfafstetter catchment +tile_coord.frac_pfaf = tmpNaN; % area fraction of Pfafstetter catchment + +if ~isEASE + + tile_coord.i_indg2 = tmpNaN; % i index of tile on global "ocean" grid (grid 2) + tile_coord.j_indg2 = tmpNaN; % j index of tile on global "ocean" grid (grid 2) + tile_coord.frac_cell2 = tmpNaN; % area fraction of "ocean" grid cell (grid 2) + +end + +% -------------------------------------------------- +% +% get indices for tile types + +ind_land = find( tile_coord.typ == 100); +ind_lake = find( tile_coord.typ == 19); +ind_landice = find( tile_coord.typ == 20); +ind_ocean = find( tile_coord.typ == 0); + +ind_NOTland = find( tile_coord.typ ~= 100); +ind_NOTocean = find( tile_coord.typ ~= 0); + +if isEASE + + col_area = NaN; + col_pfaf = 2; % land tiles only + col_frac_pfaf = NaN; % land tiles only + +else + + col_area = 2; + col_pfaf = 9; % land tiles only + col_frac_pfaf = 11; + +end + + +if ~isnan(col_area) + + tile_coord.area = tmpdata(:, col_area); % tile area + +end + +tile_coord.pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_pfaf); + +% ------------------------------------------------- +% +% EASE grid tile file: +% +% column 8: Pfafstetter index (same as column 2) +% column 9: Pfafstetter ID +% +% non-EASE grid tile file: +% +% column 9: i_indg2 [for ocean tiles] *OR* pfaf [for non-ocean tiles] +% column 10: j_indg2 +% column 11: frac_cell2 [for ocean tiles] *OR* frac_pfaf [for non-ocean tiles] +% column 12: dummy_index2 + +if ~isEASE + + tile_coord.frac_pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_frac_pfaf); + + tile_coord.i_indg2( ind_ocean) = tmpdata(ind_ocean, 9); + tile_coord.j_indg2( ind_ocean) = tmpdata(ind_ocean, 10); + tile_coord.frac_cell2(ind_ocean) = tmpdata(ind_ocean, 11); + +end + +% ------------------------------------------------- +% +% get additional info from (land-only) "catchment.def" file + +tctmp = read_catchmentdef( fname_catchmentdef ); + +% verify consistency of catchment.def file and tile file + +if tctmp.N_tile~=length(ind_land) + + error('mismatch between tile file and catchment.def file: N_tile') + +end + +if ( any(tile_coord.tile_id(ind_land) - tctmp.tile_id) | ... + any(tile_coord.pfaf( ind_land) - tctmp.pfaf ) ... + ) + + error('mismatch between tile file and catchment.def file: tile_id or pfaf') + +end + +% fill tile_coord with relevant info from catchment.def + +tile_coord.min_lon(ind_land) = tctmp.min_lon; +tile_coord.max_lon(ind_land) = tctmp.max_lon; +tile_coord.min_lat(ind_land) = tctmp.min_lat; +tile_coord.max_lat(ind_land) = tctmp.max_lat; + +tile_coord.elev( ind_land) = tctmp.elev; + + +% ========== EOF =================================================== + diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m b/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m new file mode 100644 index 0000000..e28b343 --- /dev/null +++ b/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m @@ -0,0 +1,42 @@ +function [ tile_coord ] = read_til_file_nc4( fname ) + +% read tile coordinates from nc4 "til" file (EASEv2 and cube-sphere) +% +% reichle, 24 Jan 2025 +% +% ------------------------------------------------------------- + +disp(['reading from ', fname]) + +s = ncinfo( fname ); + +tile_coord.N_tile = s.Dimensions.Length; + +varnames = {s.Variables.Name} + +for kk=1:length(varnames) + + this_varname = varnames{kk}; + + tmpdata = ncread( fname, this_varname ); + + % rename some fields + + if strcmp(this_varname,'pfaf_index' ), this_varname = 'pfaf'; end + if strcmp(this_varname,'i_indg1' ), this_varname = 'i_indg'; end + if strcmp(this_varname,'j_indg1' ), this_varname = 'j_indg'; end + if strcmp(this_varname,'frac_cell1' ), this_varname = 'frac_cell'; end + if strcmp(this_varname,'dummy_index1'), this_varname = 'dummy_index'; end + + cmd = ['tile_coord.', this_varname, ' = tmpdata;']; + + %disp(cmd) + eval(cmd) + +end + + + + + +% ========== EOF =================================================== From be03a284455f379d153dae7b042f17e44e0c710f Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 26 Jan 2025 17:22:27 -0500 Subject: [PATCH 2/9] rename tile_coord%pfaf to tile_coord%pfaf_index --- GEOSldas_App/preprocess_ldas_routines.F90 | 20 ++-- .../util/shared/matlab/read_tilecoord.m | 100 +++++++++--------- LDAS_Shared/LDAS_TileCoordType.F90 | 34 +++--- LDAS_Shared/LDAS_ensdrv_mpi.F90 | 2 +- 4 files changed, 78 insertions(+), 78 deletions(-) diff --git a/GEOSldas_App/preprocess_ldas_routines.F90 b/GEOSldas_App/preprocess_ldas_routines.F90 index 70af1d7..ea34807 100644 --- a/GEOSldas_App/preprocess_ldas_routines.F90 +++ b/GEOSldas_App/preprocess_ldas_routines.F90 @@ -409,7 +409,7 @@ subroutine domain_setup( & ! locals - integer :: n, this_tileid, this_catpfaf, N_exclude, N_include, indomain, rc + integer :: n, this_tileid, N_exclude, N_include, indomain, rc integer, dimension(N_cat_global) :: ExcludeList, IncludeList, tmp_d2g @@ -2964,7 +2964,7 @@ subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_la read (tmpline,*) & tile_coord(i)%typ, & ! 1 - tile_coord(i)%pfaf, & ! 2 + tile_coord(i)%pfaf_index, & ! 2 tile_coord(i)%com_lon, & ! 3 tile_coord(i)%com_lat, & ! 4 tile_coord(i)%i_indg, & ! 5 @@ -2985,7 +2985,7 @@ subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_la read (tmpline,*) & tile_coord(i)%typ, & ! 1 - tile_coord(i)%pfaf, & ! 2 * + tile_coord(i)%pfaf_index, & ! 2 * tile_coord(i)%com_lon, & ! 3 tile_coord(i)%com_lat, & ! 4 tile_coord(i)%i_indg, & ! 5 @@ -3009,7 +3009,7 @@ subroutine LDAS_read_til_file( tile_file, catch_file, tile_grid_g, tile_coord_la tile_coord(i)%j_indg, & ! 6 tile_coord(i)%frac_cell, & ! 7 tmpint1, & ! 8 - tile_coord(i)%pfaf, & ! 9 * + tile_coord(i)%pfaf_index, & ! 9 * tmpint2, & ! 10 tile_coord(i)%frac_pfaf, & ! 11 tmpint3 ! 12 * (previously "tile_id") @@ -3250,7 +3250,7 @@ subroutine read_catchment_def( catchment_def_file, N_tile, tile_coord ) ! ! Header line: N_tile ! - ! Columns: tile_id, Pfaf, min_lon, max_lon, min_lat, max_lat, [elev] + ! Columns: tile_id, pfaf_index, min_lon, max_lon, min_lat, max_lat, [elev] ! ! Elevation [m] is ONLY available for EASE grid tile definitions @@ -3266,7 +3266,7 @@ subroutine read_catchment_def( catchment_def_file, N_tile, tile_coord ) integer :: i, istat, tmpint1, sweep - integer, dimension(N_tile) :: tmp_tileid, tmp_pfaf + integer, dimension(N_tile) :: tmp_tileid, tmp_pfafindex character(len=*), parameter :: Iam = 'read_catchment_def' character(len=400) :: err_msg @@ -3305,7 +3305,7 @@ subroutine read_catchment_def( catchment_def_file, N_tile, tile_coord ) ! read 7 columns, avoid using exact format specification - read (10,*, iostat=istat) tmp_tileid(i), tmp_pfaf(i), & + read (10,*, iostat=istat) tmp_tileid(i), tmp_pfafindex(i), & tile_coord(i)%min_lon, & tile_coord(i)%max_lon, & tile_coord(i)%min_lat, & @@ -3316,7 +3316,7 @@ subroutine read_catchment_def( catchment_def_file, N_tile, tile_coord ) ! read 6 columns, avoid using exact format specification - read (10,*, iostat=istat) tmp_tileid(i), tmp_pfaf(i), & + read (10,*, iostat=istat) tmp_tileid(i), tmp_pfafindex(i), & tile_coord(i)%min_lon, & tile_coord(i)%max_lon, & tile_coord(i)%min_lat, & @@ -3358,8 +3358,8 @@ subroutine read_catchment_def( catchment_def_file, N_tile, tile_coord ) end do ! loop through sweeps - if ( any(tile_coord(1:N_tile)%tile_id/=tmp_tileid) .or. & - any(tile_coord(1:N_tile)%pfaf /=tmp_pfaf) ) then + if ( any(tile_coord(1:N_tile)%tile_id /=tmp_tileid) .or. & + any(tile_coord(1:N_tile)%pfaf_index/=tmp_pfafindex) ) then err_msg = 'tile_coord_file and catchment_def_file mismatch. (2)' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) diff --git a/GEOSldas_App/util/shared/matlab/read_tilecoord.m b/GEOSldas_App/util/shared/matlab/read_tilecoord.m index 422f6fb..61cb646 100644 --- a/GEOSldas_App/util/shared/matlab/read_tilecoord.m +++ b/GEOSldas_App/util/shared/matlab/read_tilecoord.m @@ -76,71 +76,71 @@ ifp = fopen( fname, 'r', machfmt); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.N_tile = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.N_tile = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); Nt = tile_coord.N_tile; - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.tile_id = fread( ifp, Nt, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.tile_id = fread( ifp, Nt, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.typ = fread( ifp, Nt, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.typ = fread( ifp, Nt, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.pfaf = fread( ifp, Nt, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.pfaf_index = fread( ifp, Nt, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.com_lon = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.com_lon = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.com_lat = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.com_lat = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.min_lon = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.min_lon = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.max_lon = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.max_lon = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.min_lat = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.min_lat = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.max_lat = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.max_lat = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.i_indg = fread( ifp, Nt, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.i_indg = fread( ifp, Nt, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.j_indg = fread( ifp, Nt, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.j_indg = fread( ifp, Nt, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.frac_cell = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.frac_cell = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.frac_pfaf = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.frac_pfaf = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.area = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.area = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); - fortran_tag = fread( ifp, 1, int_precision ); - tile_coord.elev = fread( ifp, Nt, float_precision ); - fortran_tag = fread( ifp, 1, int_precision ); + fortran_tag = fread( ifp, 1, int_precision ); + tile_coord.elev = fread( ifp, Nt, float_precision ); + fortran_tag = fread( ifp, 1, int_precision ); % if requested, convert to ASCII (txt) file @@ -166,7 +166,7 @@ '%5i%5i%13.6f%13.6f%13.4f%13.4f\n'], ... [tile_coord.tile_id(ii), ... tile_coord.typ(ii), ... - tile_coord.pfaf(ii), ... + tile_coord.pfaf_index(ii), ... tile_coord.com_lon(ii), ... tile_coord.com_lat(ii), ... tile_coord.min_lon(ii), ... @@ -208,7 +208,7 @@ tile_coord.tile_id = tmpdata(:, 1); tile_coord.typ = tmpdata(:, 2); - tile_coord.pfaf = tmpdata(:, 3); + tile_coord.pfaf_index = tmpdata(:, 3); tile_coord.com_lon = tmpdata(:, 4); tile_coord.com_lat = tmpdata(:, 5); tile_coord.min_lon = tmpdata(:, 6); diff --git a/LDAS_Shared/LDAS_TileCoordType.F90 b/LDAS_Shared/LDAS_TileCoordType.F90 index 4861ddc..f89aa5e 100644 --- a/LDAS_Shared/LDAS_TileCoordType.F90 +++ b/LDAS_Shared/LDAS_TileCoordType.F90 @@ -50,7 +50,7 @@ module LDAS_TileCoordType integer :: tile_id ! unique tile ID integer :: f_num ! full domain ID integer :: typ ! (0=MAPL_Ocean, 100=MAPL_Land, 19=MAPL_Lake, 20=MAPL_LandIce) - integer :: pfaf ! Pfafstetter number (for land tiles, NOT unique) + integer :: pfaf_index ! index of Pfafstetter catchment (for land tiles, NOT unique) real :: com_lon ! center-of-mass longitude real :: com_lat ! center-of-mass latitude real :: min_lon ! minimum longitude (bounding box for tile) @@ -398,21 +398,21 @@ subroutine io_tile_coord_type( action, unitnum, N_tile, tile_coord ) write (unitnum) N_tile - write (unitnum) (tile_coord(n)%tile_id, n=1,N_tile) - write (unitnum) (tile_coord(n)%typ, n=1,N_tile) - write (unitnum) (tile_coord(n)%pfaf, n=1,N_tile) - write (unitnum) (tile_coord(n)%com_lon, n=1,N_tile) - write (unitnum) (tile_coord(n)%com_lat, n=1,N_tile) - write (unitnum) (tile_coord(n)%min_lon, n=1,N_tile) - write (unitnum) (tile_coord(n)%max_lon, n=1,N_tile) - write (unitnum) (tile_coord(n)%min_lat, n=1,N_tile) - write (unitnum) (tile_coord(n)%max_lat, n=1,N_tile) - write (unitnum) (tile_coord(n)%i_indg, n=1,N_tile) - write (unitnum) (tile_coord(n)%j_indg, n=1,N_tile) - write (unitnum) (tile_coord(n)%frac_cell, n=1,N_tile) - write (unitnum) (tile_coord(n)%frac_pfaf, n=1,N_tile) - write (unitnum) (tile_coord(n)%area, n=1,N_tile) - write (unitnum) (tile_coord(n)%elev, n=1,N_tile) + write (unitnum) (tile_coord(n)%tile_id, n=1,N_tile) + write (unitnum) (tile_coord(n)%typ, n=1,N_tile) + write (unitnum) (tile_coord(n)%pfaf_index, n=1,N_tile) + write (unitnum) (tile_coord(n)%com_lon, n=1,N_tile) + write (unitnum) (tile_coord(n)%com_lat, n=1,N_tile) + write (unitnum) (tile_coord(n)%min_lon, n=1,N_tile) + write (unitnum) (tile_coord(n)%max_lon, n=1,N_tile) + write (unitnum) (tile_coord(n)%min_lat, n=1,N_tile) + write (unitnum) (tile_coord(n)%max_lat, n=1,N_tile) + write (unitnum) (tile_coord(n)%i_indg, n=1,N_tile) + write (unitnum) (tile_coord(n)%j_indg, n=1,N_tile) + write (unitnum) (tile_coord(n)%frac_cell, n=1,N_tile) + write (unitnum) (tile_coord(n)%frac_pfaf, n=1,N_tile) + write (unitnum) (tile_coord(n)%area, n=1,N_tile) + write (unitnum) (tile_coord(n)%elev, n=1,N_tile) case ('r','R') @@ -437,7 +437,7 @@ subroutine io_tile_coord_type( action, unitnum, N_tile, tile_coord ) read (unitnum, iostat=istat) tmp_int; if (istat>0) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) tile_coord(:)%typ = tmp_int(:) read (unitnum, iostat=istat) tmp_int; if (istat>0) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - tile_coord(:)%pfaf = tmp_int(:) + tile_coord(:)%pfaf_index = tmp_int(:) read (unitnum, iostat=istat) tmp_real; if (istat>0) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) tile_coord(:)%com_lon = tmp_real(:) diff --git a/LDAS_Shared/LDAS_ensdrv_mpi.F90 b/LDAS_Shared/LDAS_ensdrv_mpi.F90 index 07cf17b..1da40cd 100644 --- a/LDAS_Shared/LDAS_ensdrv_mpi.F90 +++ b/LDAS_Shared/LDAS_ensdrv_mpi.F90 @@ -129,7 +129,7 @@ subroutine init_MPI_types() ! integer :: tile_id ! unique tile ID ! integer :: f_num ! full domain ID ! integer :: typ ! (0=MAPL_Ocean, 100=MAPL_Land, 19=MAPL_Lake, 20=MAPL_LandIce) - ! integer :: pfaf ! Pfafstetter number (for land tiles, NOT unique) + ! integer :: pfaf_index ! index of Pfafstetter catchment (for land tiles, NOT unique) ! real :: com_lon ! center-of-mass longitude ! real :: com_lat ! center-of-mass latitude ! real :: min_lon ! minimum longitude (bounding box for tile) From 8f370ac44af6c1b342899125a54abe7feba86e17 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 26 Jan 2025 18:00:09 -0500 Subject: [PATCH 3/9] fix optional input arg in matlab ASCII tile file reader (read_tile_file_ASCII.m) --- GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m b/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m index 3dab7c3..a092e6a 100644 --- a/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m +++ b/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m @@ -25,7 +25,7 @@ % determine fname_catchmentdef -%%%if ~exist( fname_catchmentdef ) +if ~exist('fname_catchmentdef','var') ind=strfind(fname_til,'/'); @@ -41,7 +41,7 @@ fname_catchmentdef = [ fname_til(1:ind(end-2)), 'land/', resolution, '/clsm/catchment.def' ]; -%%%end +end % ------------------------------------------------------------- % From d7a9a7513255a78cd3f9e194da05ccf3f9c40667 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 26 Jan 2025 18:15:16 -0500 Subject: [PATCH 4/9] add reading of attributes to matlab nc4 tile file reader (read_tile_file_nc4.m) --- .../util/shared/matlab/read_tile_file_nc4.m | 36 ++++++++++++++++--- 1 file changed, 31 insertions(+), 5 deletions(-) diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m b/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m index e28b343..31fd1d3 100644 --- a/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m +++ b/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m @@ -10,9 +10,39 @@ s = ncinfo( fname ); +% dimensions + tile_coord.N_tile = s.Dimensions.Length; -varnames = {s.Variables.Name} +% --------------------------------- + +% read tile space attributes + +attnames = {s.Attributes.Name}; + +for kk=1:length(attnames) + + this_attname = attnames{kk}; + + % skip select attributes + + if strcmp(this_attname,'NCO' ), continue, end + if strcmp(this_attname,'history'), continue, end + + tmpdata = ncreadatt( fname, '/', this_attname ); + + cmd = ['tile_coord.', this_attname, ' = tmpdata;']; + + %disp(cmd) + eval(cmd) + +end + +% --------------------------------- + +% read tile variables + +varnames = {s.Variables.Name}; for kk=1:length(varnames) @@ -35,8 +65,4 @@ end - - - - % ========== EOF =================================================== From ca60151b43c0938a578b02a78f25ca0a89ed406b Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 27 Jan 2025 07:09:33 -0500 Subject: [PATCH 5/9] renamed matlab tile file reader (single reader for nc4 and ASCII) --- ...ead_tile_file_ASCII.m => read_tile_file.m} | 0 .../util/shared/matlab/read_tile_file_nc4.m | 68 ------------------- 2 files changed, 68 deletions(-) rename GEOSldas_App/util/shared/matlab/{read_tile_file_ASCII.m => read_tile_file.m} (100%) delete mode 100644 GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m b/GEOSldas_App/util/shared/matlab/read_tile_file.m similarity index 100% rename from GEOSldas_App/util/shared/matlab/read_tile_file_ASCII.m rename to GEOSldas_App/util/shared/matlab/read_tile_file.m diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m b/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m deleted file mode 100644 index 31fd1d3..0000000 --- a/GEOSldas_App/util/shared/matlab/read_tile_file_nc4.m +++ /dev/null @@ -1,68 +0,0 @@ -function [ tile_coord ] = read_til_file_nc4( fname ) - -% read tile coordinates from nc4 "til" file (EASEv2 and cube-sphere) -% -% reichle, 24 Jan 2025 -% -% ------------------------------------------------------------- - -disp(['reading from ', fname]) - -s = ncinfo( fname ); - -% dimensions - -tile_coord.N_tile = s.Dimensions.Length; - -% --------------------------------- - -% read tile space attributes - -attnames = {s.Attributes.Name}; - -for kk=1:length(attnames) - - this_attname = attnames{kk}; - - % skip select attributes - - if strcmp(this_attname,'NCO' ), continue, end - if strcmp(this_attname,'history'), continue, end - - tmpdata = ncreadatt( fname, '/', this_attname ); - - cmd = ['tile_coord.', this_attname, ' = tmpdata;']; - - %disp(cmd) - eval(cmd) - -end - -% --------------------------------- - -% read tile variables - -varnames = {s.Variables.Name}; - -for kk=1:length(varnames) - - this_varname = varnames{kk}; - - tmpdata = ncread( fname, this_varname ); - - % rename some fields - - if strcmp(this_varname,'pfaf_index' ), this_varname = 'pfaf'; end - if strcmp(this_varname,'i_indg1' ), this_varname = 'i_indg'; end - if strcmp(this_varname,'j_indg1' ), this_varname = 'j_indg'; end - if strcmp(this_varname,'frac_cell1' ), this_varname = 'frac_cell'; end - if strcmp(this_varname,'dummy_index1'), this_varname = 'dummy_index'; end - - cmd = ['tile_coord.', this_varname, ' = tmpdata;']; - - %disp(cmd) - eval(cmd) - -end - -% ========== EOF =================================================== From dd70d8234d691395ffab043771dd7c937642fae8 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 27 Jan 2025 07:10:36 -0500 Subject: [PATCH 6/9] merged matlab tile file readers for ASCII and nc4 into single matlab script (read_tile_file.m) --- .../util/shared/matlab/read_tile_file.m | 736 ++++++++++-------- 1 file changed, 407 insertions(+), 329 deletions(-) diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file.m b/GEOSldas_App/util/shared/matlab/read_tile_file.m index a092e6a..2312229 100644 --- a/GEOSldas_App/util/shared/matlab/read_tile_file.m +++ b/GEOSldas_App/util/shared/matlab/read_tile_file.m @@ -1,375 +1,453 @@ -function [ tile_coord ] = read_til_file_ASCII( fname_til, fname_catchmentdef ) +function [ tile_coord ] = read_til_file( fname_til, fname_catchmentdef ) -% read tile coordinates from "til" file (EASEv2 and cube-sphere) and +% read tile coordinates from tile file ("til") file (nc4 or ASCII, EASEv2 or cube-sphere) and % put data into tile_coord structure (similar to that defined in GEOSldas) % % inputs: -% fname_til : Name (with full path) of ASCII tile file (see examples below). -% fname_catchmentdef : OPTIONAL, name of matching catchment.def file -% Do NOT specify if tile file is from bcs in new layout; will be -% assembled from fname_til. +% fname_til : Name (with full path) of nc4 or ASCII tile file (see examples below). +% +% fname_catchmentdef : Name of matching catchment.def file. +% OPTIONAL, for use with ASCII tile file from bcs in legacy layout +% (assembled from fname_til if tile file is from bcs in new layout). % % fname_til examples: % '/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/v12/geometry/EASEv2_M36/EASEv2_M36_964x406.til' % '/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/v12/geometry/CF0360x6C_DE1440xPE0720/CF0360x6C_DE1440xPE0720-Pfafstetter.til' % -% % IMPORTANT: % This reader is designed to work for bcs versions NL5, v11, v12, and newer (beginning ~2022). % It should be backward-compatible for most older bcs versions, but there is no % expectation that it will work universally across all old and new versions. % -% reichle, 24 Jan 2025 +% reichle, 27 Jan 2025 % % ------------------------------------------------------------- -% determine fname_catchmentdef - -if ~exist('fname_catchmentdef','var') - - ind=strfind(fname_til,'/'); - - % verify that tile file is from bcs in new layout and naming convention - - if ~strcmp( fname_til(ind(end-2)+1:ind(end-1)-1), 'geometry' ) - - error('cannot derive fname_catchmentdef') - - end - - resolution = fname_til(ind(end-1)+1:ind(end)-1); - - fname_catchmentdef = [ fname_til(1:ind(end-2)), 'land/', resolution, '/clsm/catchment.def' ]; - -end - -% ------------------------------------------------------------- -% -% read tile file - disp(['reading from ', fname_til]) -ifp = fopen( fname_til, 'r' ); +% detect nc4 vs ASCII -% ------------------------------------------------------------- -% -% EASE grid tile space? +is_nc4 = 1; -isEASE = 0; - -if ~isempty(findstr('EASE',fname_til)), isEASE=1; end - -% -------------- -% -% determine number of data columns - -if isEASE - - % NOTE: very old EASE versions have N_data_col=8 - % - % between SMAP bcs v001 (through Vv3030) and SMAP bcs v003[a] (from Vv4030), - % a new (9th) column was added in the *.til file - - N_data_col = 9; - -else - - N_data_col = 12; - -end +if ~strcmp(fname_til(end-2:end),'nc4'), is_nc4=0; end % ------------------------------------------------------------- -% -% Number of integers in first header line depends on bcs version: - -N_headerline1 = 4; % for bcs versions after ~2022 -old_versions = { 'Fortuna', ... - 'Ganymed', ... - 'Heracles', ... - 'Icarus', ... - 'Jason', ... - 'DYAMOND2', ... - 'FV3-1_0', ... - 'NL3', ... - 'NL4', ... - }; +if is_nc4 -for kk=1:length(old_versions) + % read nc4 file - if ~isempty( strfind( fname_til, old_versions{kk} )) + s = ncinfo( fname_til ); - N_headerline1 = 3; + % dimensions - disp([ 'detected old bcs version ', old_versions{kk} ]) + tile_coord.N_tile = s.Dimensions.Length; + + % --------------------------------- + + % read tile space attributes + + attnames = {s.Attributes.Name}; + + for kk=1:length(attnames) + + this_attname = attnames{kk}; + + % skip select attributes + + if strcmp(this_attname,'NCO' ), continue, end + if strcmp(this_attname,'history'), continue, end + + tmpdata = ncreadatt( fname_til, '/', this_attname ); + + cmd = ['tile_coord.', this_attname, ' = tmpdata;']; + + %disp(cmd) + eval(cmd) end -end -% make exception for select "new" versions (stored in legacy bcs directory) that would otherwise -% be classified as "old" per the loop above - -new_versions = { 'Icarus-NLv5', ... - 'DYAMOND', ... - }; + % --------------------------------- + + % read tile variables + + varnames = {s.Variables.Name}; + + for kk=1:length(varnames) + + this_varname = varnames{kk}; + + tmpdata = ncread( fname_til, this_varname ); + + % rename some fields + + if strcmp(this_varname,'pfaf_index' ), this_varname = 'pfaf'; end + if strcmp(this_varname,'i_indg1' ), this_varname = 'i_indg'; end + if strcmp(this_varname,'j_indg1' ), this_varname = 'j_indg'; end + if strcmp(this_varname,'frac_cell1' ), this_varname = 'frac_cell'; end + if strcmp(this_varname,'dummy_index1'), this_varname = 'dummy_index'; end + + cmd = ['tile_coord.', this_varname, ' = tmpdata;']; -for kk=1:length(new_versions) + %disp(cmd) + eval(cmd) + + end - if ~isempty( strfind( fname_til, new_versions{kk} )) +else - N_headerline1 = 4; + % ------------------------------------------------------------- + % + % read ASCII file - disp([ 'detected exception for new bcs version ', new_versions{kk} ]) + % determine fname_catchmentdef + + if ~exist('fname_catchmentdef','var') + + ind=strfind(fname_til,'/'); + + % verify that tile file is from bcs in new layout and naming convention + + if ~strcmp( fname_til(ind(end-2)+1:ind(end-1)-1), 'geometry' ) + error('cannot derive fname_catchmentdef; try using optional input argument') + + end + + resolution = fname_til(ind(end-1)+1:ind(end)-1); + + fname_catchmentdef = [ fname_til(1:ind(end-2)), 'land/', resolution, '/clsm/catchment.def' ]; + end -end - -% echo final value of N_headerline1 - -disp(['using N_headerline1 = ', num2str(N_headerline1) ]) - -% ------------------------------------------------------------- -% -% read header - -% header line 1: - -tmpdata = fscanf( ifp, '%f', N_headerline1 ); - -N_tile_tmp = tmpdata(1); - -if N_headerline1==4 - - tile_coord.N_PfafCat = tmpdata(2); - tile_coord.raster_nx = tmpdata(3); - tile_coord.raster_ny = tmpdata(4); - -elseif N_headerline1==3 - - tile_coord.N_PfafCat = NaN; - tile_coord.raster_nx = tmpdata(2); - tile_coord.raster_ny = tmpdata(3); - -end - -% -------------- -% -% header line 2: - -tile_coord.N_grids = fscanf( ifp, '%f', 1 ); - -% verify N_grids (should be 1 for EASE and 2 for non-EASE) - -if isEASE & tile_coord.N_grids~=1, error(['unexpected N_grids for EASE tile space: ', num2str(tile_coord.N_grids)]), end - -if ~isEASE & tile_coord.N_grids~=2, error(['unexpected N_grids for non-EASE tile space: ', num2str(tile_coord.N_grids)]), end - -% Deal with older EASE bcs versions having (useless) header lines for grid 2, -% despite having the correct value of N_grids=1 in header line 2. - -if isEASE & ( ~isempty(findstr('NL3',fname_til)) | ~isempty(findstr('NL4',fname_til)) ) - - tmp_n_grids = 2; - -else - - tmp_n_grids = tile_coord.N_grids; - -end - -% -------------- -% -% header lines 3-5: - -tile_coord.Grid_Name = fscanf( ifp, '%s', 1 ); -tile_coord.IM = fscanf( ifp, '%f', 1 ); -tile_coord.JM = fscanf( ifp, '%f', 1 ); - -% -------------- -% -% header lines 6-8 (if present): - -if tmp_n_grids==2 - - % NOTE: EASE NL3 and NL4 contain additional header lines for second grid, despite (correct) N_grids=1 + + % ------------------------------------------------------------- + % + % read tile file + + disp(['reading from ', fname_til]) + + ifp = fopen( fname_til, 'r' ); + + % ------------------------------------------------------------- + % + % EASE grid tile space? + + isEASE = 0; + + if ~isempty(findstr('EASE',fname_til)), isEASE=1; end + + % -------------- + % + % determine number of data columns + + if isEASE + + % NOTE: very old EASE versions have N_data_col=8 + % + % between SMAP bcs v001 (through Vv3030) and SMAP bcs v003[a] (from Vv4030), + % a new (9th) column was added in the *.til file + + N_data_col = 9; + + else + + N_data_col = 12; + + end + + % ------------------------------------------------------------- + % + % Number of integers in first header line depends on bcs version: + + N_headerline1 = 4; % for bcs versions after ~2022 + + old_versions = { 'Fortuna', ... + 'Ganymed', ... + 'Heracles', ... + 'Icarus', ... + 'Jason', ... + 'DYAMOND2', ... + 'FV3-1_0', ... + 'NL3', ... + 'NL4', ... + }; + + for kk=1:length(old_versions) + + if ~isempty( strfind( fname_til, old_versions{kk} )) + + N_headerline1 = 3; + + disp([ 'detected old bcs version ', old_versions{kk} ]) + + end + end + + % make exception for select "new" versions (stored in legacy bcs directory) that would otherwise + % be classified as "old" per the loop above + + new_versions = { 'Icarus-NLv5', ... + 'DYAMOND', ... + }; + + for kk=1:length(new_versions) + + if ~isempty( strfind( fname_til, new_versions{kk} )) + + N_headerline1 = 4; + + disp([ 'detected exception for new bcs version ', new_versions{kk} ]) + + end + end + + % echo final value of N_headerline1 + + disp(['using N_headerline1 = ', num2str(N_headerline1) ]) + + % ------------------------------------------------------------- + % + % read header + + % header line 1: + + tmpdata = fscanf( ifp, '%f', N_headerline1 ); + + N_tile_tmp = tmpdata(1); - tile_coord.Grid2_Name = fscanf( ifp, '%s', 1 ); - tile_coord.IM2 = fscanf( ifp, '%f', 1 ); + if N_headerline1==4 + + tile_coord.N_PfafCat = tmpdata(2); + tile_coord.raster_nx = tmpdata(3); + tile_coord.raster_ny = tmpdata(4); + + elseif N_headerline1==3 + + tile_coord.N_PfafCat = NaN; + tile_coord.raster_nx = tmpdata(2); + tile_coord.raster_ny = tmpdata(3); + + end + + % -------------- + % + % header line 2: + + tile_coord.N_grids = fscanf( ifp, '%f', 1 ); + + % verify N_grids (should be 1 for EASE and 2 for non-EASE) + + if isEASE & tile_coord.N_grids~=1, error(['unexpected N_grids for EASE tile space: ', num2str(tile_coord.N_grids)]), end + + if ~isEASE & tile_coord.N_grids~=2, error(['unexpected N_grids for non-EASE tile space: ', num2str(tile_coord.N_grids)]), end + + % Deal with older EASE bcs versions having (useless) header lines for grid 2, + % despite having the correct value of N_grids=1 in header line 2. + + if isEASE & ( ~isempty(findstr('NL3',fname_til)) | ~isempty(findstr('NL4',fname_til)) ) + + tmp_n_grids = 2; + + else + + tmp_n_grids = tile_coord.N_grids; + + end + + % -------------- + % + % header lines 3-5: + + tile_coord.Grid_Name = fscanf( ifp, '%s', 1 ); + tile_coord.IM = fscanf( ifp, '%f', 1 ); + tile_coord.JM = fscanf( ifp, '%f', 1 ); + + % -------------- + % + % header lines 6-8 (if present): + + if tmp_n_grids==2 + + % NOTE: EASE NL3 and NL4 contain additional header lines for second grid, despite (correct) N_grids=1 + + tile_coord.Grid2_Name = fscanf( ifp, '%s', 1 ); + tile_coord.IM2 = fscanf( ifp, '%f', 1 ); tile_coord.JM2 = fscanf( ifp, '%f', 1 ); -end - -% -------------- -% -% read data - -tmpdata = fscanf( ifp, '%f' ); - -fclose(ifp); - -disp('done reading "til" file') - -% verify N_tile*N_data_col against number of data read - -if length(tmpdata)~=N_tile_tmp*N_data_col, error('something wrong with N_tile*N_data_col'), end - -% convert into 2d array - -tmpdata = transpose( reshape(tmpdata, [N_data_col, N_tile_tmp]) ); - -% -------------------------------------------------- -% -% assign tile_id -% -% by convention, tile_id is equal to index number of tile in *.til file - -tmp_tileid = transpose(1:N_tile_tmp); - -% -------------------------------------------------- -% -% subset (if requested) - -%ind = find(tmpdata(:,1)==100); % keep only land -%ind = find(tmpdata(:,1)== 19); % keep only lakes -%ind = find(tmpdata(:,1)== 20); % keep only ice -%ind = find(tmpdata(:,1)~= 0); % keep land, lake, landice - -%tmpdata = tmpdata( ind,:); - -%tmp_tileid = tmp_tileid( ind, 1); - -% -------------------------------------------------- -% -% copy data into tile_coord structure - -tile_coord.N_tile = size(tmpdata,1); % number of tiles (assign here in case of subsetting above) - -% the following are universal (for EASE and non-EASE, all tile types) - -tile_coord.tile_id = tmp_tileid; % tile ID - -tile_coord.typ = tmpdata(:, 1); % tile type - -tile_coord.com_lon = tmpdata(:, 3); % center-of-mass longitude of tile -tile_coord.com_lat = tmpdata(:, 4); % center-of-mass latitude of tile - -tile_coord.i_indg = tmpdata(:, 5); % i index of tile on global "atm" (or EASE) grid -tile_coord.j_indg = tmpdata(:, 6); % j index of tile on global "atm" (or EASE) grid -tile_coord.frac_cell = tmpdata(:, 7); % area fraction of "atm" (or EASE) grid cell - -% initialize remaining fields (to be filled below) - -tmpNaN = NaN*ones(tile_coord.N_tile,1); - -tile_coord.min_lon = tmpNaN; % min longitude of tile -tile_coord.max_lon = tmpNaN; % max longitude of tile -tile_coord.min_lat = tmpNaN; % min latitude of tile -tile_coord.max_lat = tmpNaN; % max latitude of tile - -tile_coord.elev = tmpNaN; % elevation of tile - -tile_coord.area = tmpNaN; % area of "atm" grid cell -tile_coord.pfaf = tmpNaN; % index of (hydrological) Pfafstetter catchment -tile_coord.frac_pfaf = tmpNaN; % area fraction of Pfafstetter catchment - -if ~isEASE - - tile_coord.i_indg2 = tmpNaN; % i index of tile on global "ocean" grid (grid 2) - tile_coord.j_indg2 = tmpNaN; % j index of tile on global "ocean" grid (grid 2) - tile_coord.frac_cell2 = tmpNaN; % area fraction of "ocean" grid cell (grid 2) - -end - -% -------------------------------------------------- -% -% get indices for tile types - -ind_land = find( tile_coord.typ == 100); -ind_lake = find( tile_coord.typ == 19); -ind_landice = find( tile_coord.typ == 20); -ind_ocean = find( tile_coord.typ == 0); - -ind_NOTland = find( tile_coord.typ ~= 100); -ind_NOTocean = find( tile_coord.typ ~= 0); - -if isEASE - - col_area = NaN; - col_pfaf = 2; % land tiles only - col_frac_pfaf = NaN; % land tiles only - -else - - col_area = 2; - col_pfaf = 9; % land tiles only - col_frac_pfaf = 11; + end -end - - -if ~isnan(col_area) - - tile_coord.area = tmpdata(:, col_area); % tile area + % -------------- + % + % read data -end - -tile_coord.pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_pfaf); - -% ------------------------------------------------- -% -% EASE grid tile file: -% -% column 8: Pfafstetter index (same as column 2) -% column 9: Pfafstetter ID -% -% non-EASE grid tile file: -% -% column 9: i_indg2 [for ocean tiles] *OR* pfaf [for non-ocean tiles] -% column 10: j_indg2 -% column 11: frac_cell2 [for ocean tiles] *OR* frac_pfaf [for non-ocean tiles] -% column 12: dummy_index2 - -if ~isEASE - - tile_coord.frac_pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_frac_pfaf); - - tile_coord.i_indg2( ind_ocean) = tmpdata(ind_ocean, 9); - tile_coord.j_indg2( ind_ocean) = tmpdata(ind_ocean, 10); - tile_coord.frac_cell2(ind_ocean) = tmpdata(ind_ocean, 11); - -end - -% ------------------------------------------------- -% -% get additional info from (land-only) "catchment.def" file - -tctmp = read_catchmentdef( fname_catchmentdef ); - -% verify consistency of catchment.def file and tile file - -if tctmp.N_tile~=length(ind_land) - - error('mismatch between tile file and catchment.def file: N_tile') - -end - -if ( any(tile_coord.tile_id(ind_land) - tctmp.tile_id) | ... - any(tile_coord.pfaf( ind_land) - tctmp.pfaf ) ... - ) - - error('mismatch between tile file and catchment.def file: tile_id or pfaf') - -end - -% fill tile_coord with relevant info from catchment.def - -tile_coord.min_lon(ind_land) = tctmp.min_lon; -tile_coord.max_lon(ind_land) = tctmp.max_lon; -tile_coord.min_lat(ind_land) = tctmp.min_lat; -tile_coord.max_lat(ind_land) = tctmp.max_lat; - -tile_coord.elev( ind_land) = tctmp.elev; + tmpdata = fscanf( ifp, '%f' ); + + fclose(ifp); + + disp('done reading "til" file') + + % verify N_tile*N_data_col against number of data read + + if length(tmpdata)~=N_tile_tmp*N_data_col, error('something wrong with N_tile*N_data_col'), end + + % convert into 2d array + + tmpdata = transpose( reshape(tmpdata, [N_data_col, N_tile_tmp]) ); + + % -------------------------------------------------- + % + % assign tile_id + % + % by convention, tile_id is equal to index number of tile in *.til file + + tmp_tileid = transpose(1:N_tile_tmp); + + % -------------------------------------------------- + % + % subset (if requested) + + %ind = find(tmpdata(:,1)==100); % keep only land + %ind = find(tmpdata(:,1)== 19); % keep only lakes + %ind = find(tmpdata(:,1)== 20); % keep only ice + %ind = find(tmpdata(:,1)~= 0); % keep land, lake, landice + + %tmpdata = tmpdata( ind,:); + + %tmp_tileid = tmp_tileid( ind, 1); + + % -------------------------------------------------- + % + % copy data into tile_coord structure + + tile_coord.N_tile = size(tmpdata,1); % number of tiles (assign here in case of subsetting above) + + % the following are universal (for EASE and non-EASE, all tile types) + + tile_coord.tile_id = tmp_tileid; % tile ID + + tile_coord.typ = tmpdata(:, 1); % tile type + + tile_coord.com_lon = tmpdata(:, 3); % center-of-mass longitude of tile + tile_coord.com_lat = tmpdata(:, 4); % center-of-mass latitude of tile + + tile_coord.i_indg = tmpdata(:, 5); % i index of tile on global "atm" (or EASE) grid + tile_coord.j_indg = tmpdata(:, 6); % j index of tile on global "atm" (or EASE) grid + tile_coord.frac_cell = tmpdata(:, 7); % area fraction of "atm" (or EASE) grid cell + + % initialize remaining fields (to be filled below) + + tmpNaN = NaN*ones(tile_coord.N_tile,1); + + tile_coord.min_lon = tmpNaN; % min longitude of tile + tile_coord.max_lon = tmpNaN; % max longitude of tile + tile_coord.min_lat = tmpNaN; % min latitude of tile + tile_coord.max_lat = tmpNaN; % max latitude of tile + + tile_coord.elev = tmpNaN; % elevation of tile + + tile_coord.area = tmpNaN; % area of "atm" grid cell + tile_coord.pfaf = tmpNaN; % index of (hydrological) Pfafstetter catchment + tile_coord.frac_pfaf = tmpNaN; % area fraction of Pfafstetter catchment + + if ~isEASE + + tile_coord.i_indg2 = tmpNaN; % i index of tile on global "ocean" grid (grid 2) + tile_coord.j_indg2 = tmpNaN; % j index of tile on global "ocean" grid (grid 2) + tile_coord.frac_cell2 = tmpNaN; % area fraction of "ocean" grid cell (grid 2) + + end + + % -------------------------------------------------- + % + % get indices for tile types + + ind_land = find( tile_coord.typ == 100); + ind_lake = find( tile_coord.typ == 19); + ind_landice = find( tile_coord.typ == 20); + ind_ocean = find( tile_coord.typ == 0); + + ind_NOTland = find( tile_coord.typ ~= 100); + ind_NOTocean = find( tile_coord.typ ~= 0); + + if isEASE + + col_area = NaN; + col_pfaf = 2; % land tiles only + col_frac_pfaf = NaN; % land tiles only + + else + + col_area = 2; + col_pfaf = 9; % land tiles only + col_frac_pfaf = 11; + + end + + + if ~isnan(col_area) + + tile_coord.area = tmpdata(:, col_area); % tile area + + end + + tile_coord.pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_pfaf); + + % ------------------------------------------------- + % + % EASE grid tile file: + % + % column 8: Pfafstetter index (same as column 2) + % column 9: Pfafstetter ID + % + % non-EASE grid tile file: + % + % column 9: i_indg2 [for ocean tiles] *OR* pfaf [for non-ocean tiles] + % column 10: j_indg2 + % column 11: frac_cell2 [for ocean tiles] *OR* frac_pfaf [for non-ocean tiles] + % column 12: dummy_index2 + + if ~isEASE + + tile_coord.frac_pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_frac_pfaf); + + tile_coord.i_indg2( ind_ocean) = tmpdata(ind_ocean, 9); + tile_coord.j_indg2( ind_ocean) = tmpdata(ind_ocean, 10); + tile_coord.frac_cell2(ind_ocean) = tmpdata(ind_ocean, 11); + + end + + % ------------------------------------------------- + % + % get additional info from (land-only) "catchment.def" file + + tctmp = read_catchmentdef( fname_catchmentdef ); + + % verify consistency of catchment.def file and tile file + + if tctmp.N_tile~=length(ind_land) + + error('mismatch between tile file and catchment.def file: N_tile') + + end + + if ( any(tile_coord.tile_id(ind_land) - tctmp.tile_id) | ... + any(tile_coord.pfaf( ind_land) - tctmp.pfaf ) ... + ) + + error('mismatch between tile file and catchment.def file: tile_id or pfaf') + + end + + % fill tile_coord with relevant info from catchment.def + + tile_coord.min_lon(ind_land) = tctmp.min_lon; + tile_coord.max_lon(ind_land) = tctmp.max_lon; + tile_coord.min_lat(ind_land) = tctmp.min_lat; + tile_coord.max_lat(ind_land) = tctmp.max_lat; + + tile_coord.elev( ind_land) = tctmp.elev; +end % if is_nc4 % ========== EOF =================================================== From 0a60096de1c07e9a27c1c585a5133efc736b7580 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 27 Jan 2025 07:26:09 -0500 Subject: [PATCH 7/9] updated CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 390b789..9e38816 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Updated subroutine read_obs_sm_ASCAT_EUMET() to work with both original and revised file name templates. - Updated subroutines read_obs_sm_ASCAT_EUMET(), read_obs_SMAP_halforbit_Tb(), read_obs_SMOS() and read_obs_MODIS_SCF() with hardcoded time ranges for when observations are available and should be read. - Revised variable names (SHORT_NAME) and descriptions (LONG_NAME) to match M21C file specs. +- Renamed tilecoord%pfaf to %pfaf_index; added matlab tile file reader. ### Fixed From fa232ef9c7ddac203f218943df81b28ad5fbccf7 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 27 Jan 2025 09:37:20 -0500 Subject: [PATCH 8/9] added matlab reader for catchment.def file; fixed name of tile coord field "pfaf_index" in read_tile_file.m --- .../util/shared/matlab/read_catchmentdef.m | 46 +++++++++++++++++++ .../util/shared/matlab/read_tile_file.m | 25 +++++----- 2 files changed, 58 insertions(+), 13 deletions(-) create mode 100644 GEOSldas_App/util/shared/matlab/read_catchmentdef.m diff --git a/GEOSldas_App/util/shared/matlab/read_catchmentdef.m b/GEOSldas_App/util/shared/matlab/read_catchmentdef.m new file mode 100644 index 0000000..3a12c93 --- /dev/null +++ b/GEOSldas_App/util/shared/matlab/read_catchmentdef.m @@ -0,0 +1,46 @@ +function [tile_coord ] = read_catchmentdef( fname ) + +% read (land) tile properties from "catchment.def" file (EASEv2 and cube-sphere) +% +% reichle, 24 Jan 2025 +% +% ------------------------------------------------------------- + +% read file + +disp(['reading from ', fname]) + +ifp = fopen( fname, 'r' ); + +% read header line + +tmpdata = fscanf( ifp, '%f', 1 ); + +tile_coord.N_tile = tmpdata(1); + +% read rest of data + +tmpdata = fscanf( ifp, '%f' ); + +fclose(ifp); + +disp('done reading "catchment.def" file') + +% -------------------------------------------------- + +% process data + +tmpdata = reshape(tmpdata, [7, tile_coord.N_tile])'; + +tile_coord.tile_id = tmpdata(:,1); +tile_coord.pfaf_index = tmpdata(:,2); + +tile_coord.min_lon = tmpdata(:,3); +tile_coord.max_lon = tmpdata(:,4); +tile_coord.min_lat = tmpdata(:,5); +tile_coord.max_lat = tmpdata(:,6); + +tile_coord.elev = tmpdata(:,7); + +% ========== EOF =================================================== + diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file.m b/GEOSldas_App/util/shared/matlab/read_tile_file.m index 2312229..e4fff6b 100644 --- a/GEOSldas_App/util/shared/matlab/read_tile_file.m +++ b/GEOSldas_App/util/shared/matlab/read_tile_file.m @@ -81,7 +81,6 @@ % rename some fields - if strcmp(this_varname,'pfaf_index' ), this_varname = 'pfaf'; end if strcmp(this_varname,'i_indg1' ), this_varname = 'i_indg'; end if strcmp(this_varname,'j_indg1' ), this_varname = 'j_indg'; end if strcmp(this_varname,'frac_cell1' ), this_varname = 'frac_cell'; end @@ -346,7 +345,7 @@ tile_coord.elev = tmpNaN; % elevation of tile tile_coord.area = tmpNaN; % area of "atm" grid cell - tile_coord.pfaf = tmpNaN; % index of (hydrological) Pfafstetter catchment + tile_coord.pfaf_index = tmpNaN; % index of (hydrological) Pfafstetter catchment tile_coord.frac_pfaf = tmpNaN; % area fraction of Pfafstetter catchment if ~isEASE @@ -372,14 +371,14 @@ if isEASE col_area = NaN; - col_pfaf = 2; % land tiles only - col_frac_pfaf = NaN; % land tiles only + col_pfafindex = 2; % land tiles only + col_fracpfaf = NaN; % land tiles only else col_area = 2; - col_pfaf = 9; % land tiles only - col_frac_pfaf = 11; + col_pfafindex = 9; % land tiles only + col_fracpfaf = 11; end @@ -390,7 +389,7 @@ end - tile_coord.pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_pfaf); + tile_coord.pfaf_index( ind_NOTocean) = tmpdata(ind_NOTocean, col_pfafindex); % ------------------------------------------------- % @@ -401,14 +400,14 @@ % % non-EASE grid tile file: % - % column 9: i_indg2 [for ocean tiles] *OR* pfaf [for non-ocean tiles] + % column 9: i_indg2 [for ocean tiles] *OR* pfaf_index [for non-ocean tiles] % column 10: j_indg2 - % column 11: frac_cell2 [for ocean tiles] *OR* frac_pfaf [for non-ocean tiles] + % column 11: frac_cell2 [for ocean tiles] *OR* frac_pfaf [for non-ocean tiles] % column 12: dummy_index2 if ~isEASE - tile_coord.frac_pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_frac_pfaf); + tile_coord.frac_pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_fracpfaf); tile_coord.i_indg2( ind_ocean) = tmpdata(ind_ocean, 9); tile_coord.j_indg2( ind_ocean) = tmpdata(ind_ocean, 10); @@ -430,11 +429,11 @@ end - if ( any(tile_coord.tile_id(ind_land) - tctmp.tile_id) | ... - any(tile_coord.pfaf( ind_land) - tctmp.pfaf ) ... + if ( any(tile_coord.tile_id( ind_land) - tctmp.tile_id ) | ... + any(tile_coord.pfaf_index(ind_land) - tctmp.pfaf_index) ... ) - error('mismatch between tile file and catchment.def file: tile_id or pfaf') + error('mismatch between tile file and catchment.def file: tile_id or pfaf_index') end From 5496e35d64df5e01820d3b195b3d638712675fe6 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 29 Jan 2025 13:06:59 -0500 Subject: [PATCH 9/9] in matlab tile file reader, synchronize field names with latest make_bcs commits (read_tile_file.m) --- .../util/shared/matlab/read_tile_file.m | 72 +++++++++---------- 1 file changed, 33 insertions(+), 39 deletions(-) diff --git a/GEOSldas_App/util/shared/matlab/read_tile_file.m b/GEOSldas_App/util/shared/matlab/read_tile_file.m index e4fff6b..874d689 100644 --- a/GEOSldas_App/util/shared/matlab/read_tile_file.m +++ b/GEOSldas_App/util/shared/matlab/read_tile_file.m @@ -79,13 +79,6 @@ tmpdata = ncread( fname_til, this_varname ); - % rename some fields - - if strcmp(this_varname,'i_indg1' ), this_varname = 'i_indg'; end - if strcmp(this_varname,'j_indg1' ), this_varname = 'j_indg'; end - if strcmp(this_varname,'frac_cell1' ), this_varname = 'frac_cell'; end - if strcmp(this_varname,'dummy_index1'), this_varname = 'dummy_index'; end - cmd = ['tile_coord.', this_varname, ' = tmpdata;']; %disp(cmd) @@ -232,16 +225,16 @@ % % header line 2: - tile_coord.N_grids = fscanf( ifp, '%f', 1 ); + tile_coord.N_Grids = fscanf( ifp, '%f', 1 ); - % verify N_grids (should be 1 for EASE and 2 for non-EASE) + % verify N_Grids (should be 1 for EASE and 2 for non-EASE) - if isEASE & tile_coord.N_grids~=1, error(['unexpected N_grids for EASE tile space: ', num2str(tile_coord.N_grids)]), end + if isEASE & tile_coord.N_Grids~=1, error(['unexpected N_Grids for EASE tile space: ', num2str(tile_coord.N_Grids)]), end - if ~isEASE & tile_coord.N_grids~=2, error(['unexpected N_grids for non-EASE tile space: ', num2str(tile_coord.N_grids)]), end + if ~isEASE & tile_coord.N_Grids~=2, error(['unexpected N_Grids for non-EASE tile space: ', num2str(tile_coord.N_Grids)]), end % Deal with older EASE bcs versions having (useless) header lines for grid 2, - % despite having the correct value of N_grids=1 in header line 2. + % despite having the correct value of N_Grids=1 in header line 2. if isEASE & ( ~isempty(findstr('NL3',fname_til)) | ~isempty(findstr('NL4',fname_til)) ) @@ -249,7 +242,7 @@ else - tmp_n_grids = tile_coord.N_grids; + tmp_n_grids = tile_coord.N_Grids; end @@ -267,11 +260,12 @@ if tmp_n_grids==2 - % NOTE: EASE NL3 and NL4 contain additional header lines for second grid, despite (correct) N_grids=1 + % NOTE: EASE NL3 and NL4 contain additional header lines for second grid, despite (correct) N_Grids=1 + % For these old versions, read (and then ignore) additional header lines. - tile_coord.Grid2_Name = fscanf( ifp, '%s', 1 ); - tile_coord.IM2 = fscanf( ifp, '%f', 1 ); - tile_coord.JM2 = fscanf( ifp, '%f', 1 ); + tile_coord.Grid_ocn_Name = fscanf( ifp, '%s', 1 ); + tile_coord.IM_ocn = fscanf( ifp, '%f', 1 ); + tile_coord.JM_ocn = fscanf( ifp, '%f', 1 ); end @@ -318,41 +312,41 @@ % % copy data into tile_coord structure - tile_coord.N_tile = size(tmpdata,1); % number of tiles (assign here in case of subsetting above) + tile_coord.N_tile = size(tmpdata,1); % number of tiles (assign here in case of subsetting above) % the following are universal (for EASE and non-EASE, all tile types) - tile_coord.tile_id = tmp_tileid; % tile ID + tile_coord.tile_id = tmp_tileid; % tile ID - tile_coord.typ = tmpdata(:, 1); % tile type + tile_coord.typ = tmpdata(:, 1); % tile type - tile_coord.com_lon = tmpdata(:, 3); % center-of-mass longitude of tile - tile_coord.com_lat = tmpdata(:, 4); % center-of-mass latitude of tile + tile_coord.com_lon = tmpdata(:, 3); % center-of-mass longitude of tile + tile_coord.com_lat = tmpdata(:, 4); % center-of-mass latitude of tile - tile_coord.i_indg = tmpdata(:, 5); % i index of tile on global "atm" (or EASE) grid - tile_coord.j_indg = tmpdata(:, 6); % j index of tile on global "atm" (or EASE) grid - tile_coord.frac_cell = tmpdata(:, 7); % area fraction of "atm" (or EASE) grid cell + tile_coord.i_indg = tmpdata(:, 5); % i index of tile on global "atm" (or EASE) grid + tile_coord.j_indg = tmpdata(:, 6); % j index of tile on global "atm" (or EASE) grid + tile_coord.frac_cell = tmpdata(:, 7); % area fraction of "atm" (or EASE) grid cell % initialize remaining fields (to be filled below) tmpNaN = NaN*ones(tile_coord.N_tile,1); - tile_coord.min_lon = tmpNaN; % min longitude of tile - tile_coord.max_lon = tmpNaN; % max longitude of tile - tile_coord.min_lat = tmpNaN; % min latitude of tile - tile_coord.max_lat = tmpNaN; % max latitude of tile + tile_coord.min_lon = tmpNaN; % min longitude of tile + tile_coord.max_lon = tmpNaN; % max longitude of tile + tile_coord.min_lat = tmpNaN; % min latitude of tile + tile_coord.max_lat = tmpNaN; % max latitude of tile - tile_coord.elev = tmpNaN; % elevation of tile + tile_coord.elev = tmpNaN; % elevation of tile - tile_coord.area = tmpNaN; % area of "atm" grid cell - tile_coord.pfaf_index = tmpNaN; % index of (hydrological) Pfafstetter catchment - tile_coord.frac_pfaf = tmpNaN; % area fraction of Pfafstetter catchment + tile_coord.area = tmpNaN; % area of "atm" grid cell + tile_coord.pfaf_index = tmpNaN; % index of (hydrological) Pfafstetter catchment + tile_coord.frac_pfaf = tmpNaN; % area fraction of Pfafstetter catchment if ~isEASE - tile_coord.i_indg2 = tmpNaN; % i index of tile on global "ocean" grid (grid 2) - tile_coord.j_indg2 = tmpNaN; % j index of tile on global "ocean" grid (grid 2) - tile_coord.frac_cell2 = tmpNaN; % area fraction of "ocean" grid cell (grid 2) + tile_coord.i_indg_ocn = tmpNaN; % i index of tile on global "ocean" grid (grid 2) + tile_coord.j_indg_ocn = tmpNaN; % j index of tile on global "ocean" grid (grid 2) + tile_coord.frac_cell_ocn = tmpNaN; % area fraction of "ocean" grid cell (grid 2) end @@ -409,9 +403,9 @@ tile_coord.frac_pfaf( ind_NOTocean) = tmpdata(ind_NOTocean, col_fracpfaf); - tile_coord.i_indg2( ind_ocean) = tmpdata(ind_ocean, 9); - tile_coord.j_indg2( ind_ocean) = tmpdata(ind_ocean, 10); - tile_coord.frac_cell2(ind_ocean) = tmpdata(ind_ocean, 11); + tile_coord.i_indg_ocn( ind_ocean) = tmpdata(ind_ocean, 9); + tile_coord.j_indg_ocn( ind_ocean) = tmpdata(ind_ocean, 10); + tile_coord.frac_cell_ocn(ind_ocean) = tmpdata(ind_ocean, 11); end