diff --git a/CHANGELOG.md b/CHANGELOG.md index 6ee77c0..38bbef5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,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 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_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 new file mode 100644 index 0000000..874d689 --- /dev/null +++ b/GEOSldas_App/util/shared/matlab/read_tile_file.m @@ -0,0 +1,446 @@ +function [ tile_coord ] = read_til_file( fname_til, fname_catchmentdef ) + +% 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 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, 27 Jan 2025 +% +% ------------------------------------------------------------- + +disp(['reading from ', fname_til]) + +% detect nc4 vs ASCII + +is_nc4 = 1; + +if ~strcmp(fname_til(end-2:end),'nc4'), is_nc4=0; end + +% ------------------------------------------------------------- + +if is_nc4 + + % read nc4 file + + s = ncinfo( fname_til ); + + % 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_til, '/', 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_til, this_varname ); + + cmd = ['tile_coord.', this_varname, ' = tmpdata;']; + + %disp(cmd) + eval(cmd) + + end + +else + + % ------------------------------------------------------------- + % + % read ASCII file + + % 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 + + % ------------------------------------------------------------- + % + % 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 + % For these old versions, read (and then ignore) additional header lines. + + 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 + + % -------------- + % + % 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_index = tmpNaN; % index of (hydrological) Pfafstetter catchment + tile_coord.frac_pfaf = tmpNaN; % area fraction of Pfafstetter catchment + + if ~isEASE + + 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 + + % -------------------------------------------------- + % + % 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_pfafindex = 2; % land tiles only + col_fracpfaf = NaN; % land tiles only + + else + + col_area = 2; + col_pfafindex = 9; % land tiles only + col_fracpfaf = 11; + + end + + + if ~isnan(col_area) + + tile_coord.area = tmpdata(:, col_area); % tile area + + end + + tile_coord.pfaf_index( ind_NOTocean) = tmpdata(ind_NOTocean, col_pfafindex); + + % ------------------------------------------------- + % + % 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_index [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_fracpfaf); + + 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 + + % ------------------------------------------------- + % + % 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_index(ind_land) - tctmp.pfaf_index) ... + ) + + error('mismatch between tile file and catchment.def file: tile_id or pfaf_index') + + 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 =================================================== + 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)