Skip to content

Commit de453ad

Browse files
biljanaorescaninbiljanaorescanin
authored
Merge branch 'develop' into feature/wjiang/fix_remap
2 parents bdf3bc8 + 8185f71 commit de453ad

File tree

7 files changed

+216
-35
lines changed

7 files changed

+216
-35
lines changed

‎GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -55,6 +55,7 @@ file(GLOB MAKE_BCS_PYTHON CONFIGURE_DEPENDS "./make_bcs*.py")
5555
list(FILTER MAKE_BCS_PYTHON EXCLUDE REGEX "make_bcs_shared.py")
5656
install(PROGRAMS ${MAKE_BCS_PYTHON} DESTINATION bin)
5757
install(PROGRAMS TileFile_ASCII_to_nc4.py DESTINATION bin)
58+
install(PROGRAMS ExtractBCsFromOrig.py DESTINATION bin)
5859

5960
set(file ./make_bcs_shared.py)
6061
configure_file(${file} ${file} @ONLY)

‎GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/ExtractBCsFromOrig.py‎

Lines changed: 163 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,163 @@
1+
#!/usr/bin/env python3
2+
3+
from netCDF4 import Dataset
4+
from scipy.io import FortranFile
5+
import numpy as np
6+
import sys
7+
import os
8+
import glob
9+
10+
def get_i_j_pf(til_file):
11+
"""
12+
This is just for cubed-sphere text tile.
13+
"""
14+
ii = []
15+
jj = []
16+
pfaf = []
17+
with open(til_file, 'r') as f:
18+
for line in f:
19+
data = line.split()
20+
if data[0] == '100' :
21+
ii.append(int(data[4]))
22+
jj.append(int(data[5]))
23+
pfaf.append(int(data[8]))
24+
return ii, jj, pfaf
25+
26+
def create_reduced_nc4(input_file, output_file, indices):
27+
"""
28+
Create a new nc4 file with with indices of for the input file.
29+
30+
Parameters:
31+
input_file (str): Path to the input nc4 file
32+
output_file (str): Path to the output nc4 file
33+
indices(int array):
34+
"""
35+
36+
# Open the input file in read mode
37+
try:
38+
src = Dataset(input_file, 'r')
39+
print(f"Successfully opened input file: {input_file}")
40+
except Exception as e:
41+
print(f"Error opening input file: {e}")
42+
return
43+
44+
# Create the output file in write mode
45+
try:
46+
dst = Dataset(output_file, 'w', format='NETCDF4')
47+
print(f"Creating output file: {output_file}")
48+
except Exception as e:
49+
print(f"Error creating output file: {e}")
50+
src.close()
51+
return
52+
53+
# Copy global attributes
54+
dst.setncatts(src.__dict__)
55+
56+
# Copy dimensions with reduced size for unlimited/time dimensions
57+
for name, dimension in src.dimensions.items():
58+
if name == 'tile':
59+
dst.createDimension(name, len(indices))
60+
else:
61+
dst.createDimension(name, len(dimension))
62+
# Copy variables with reduced data
63+
for name, variable in src.variables.items():
64+
# Create variable with same attributes but possibly reduced dimensions
65+
dims = variable.dimensions
66+
new_dims = []
67+
for dim in dims:
68+
new_dims.append(dim)
69+
# Create the variable in the output file
70+
dst_var = dst.createVariable(
71+
name,
72+
variable.datatype,
73+
new_dims,
74+
)
75+
76+
# Copy variable attributes
77+
dst[name].setncatts(variable.__dict__)
78+
79+
# Get the data and reduce it
80+
if len(dims) == 2 :
81+
dst_var[:,:] = variable[:,indices]
82+
if len(dims) == 1 :
83+
dst_var[:] = variable[indices]
84+
85+
86+
# Close both files
87+
src.close()
88+
dst.close()
89+
print(f"Successfully created reduced file: {output_file}")
90+
91+
def create_reduced_bin(input_file, output_file, indices):
92+
"""
93+
Create a new ibinary file with with indices of for the input file.
94+
95+
Parameters:
96+
input_file (str): Path to the input binary file
97+
output_file (str): Path to the output binaryfile
98+
indices(int array):
99+
"""
100+
nland = len(indices)
101+
fout = FortranFile(output_file, 'w')
102+
with FortranFile(input_file, 'r') as f:
103+
while True :
104+
try:
105+
a = f.read_reals(dtype=np.float32)
106+
b = f.read_reals(dtype=np.float32)
107+
a[12] = np.float32(nland)
108+
fout.write_record(a)
109+
b = b[indices]
110+
fout.write_record(b)
111+
except :
112+
break
113+
fout.close()
114+
print(f"Successfully created reduced file: {output_file}")
115+
116+
# Example usage
117+
if __name__ == "__main__":
118+
# Replace these with your actual file paths
119+
# input_nc4 = "/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/v13/land/CF1440x6C/clsm/catchcn_params.nc4"
120+
# output_nc4 = "output_file_reduced.nc"
121+
bcs_dir = sys.argv[1]
122+
bcs_ver = sys.argv[2]
123+
air_res = sys.argv[3]
124+
ocn_res = sys.argv[4]
125+
agname = air_res + '_' + air_res
126+
orig_tile = bcs_dir + '/' + bcs_ver + '/geometry/' + agname + '/' + agname +'-Pfafstetter.til'
127+
if not os.path.exists(orig_tile) :
128+
print( "The original tile file must exist " + orig_tile)
129+
exit()
130+
131+
mom_tile = 'til/' + air_res + '_' + ocn_res + '-Pfafstetter.til'
132+
if not os.path.exists(mom_tile) :
133+
print( "The MOM tile file must exist " + mom_tile)
134+
exit()
135+
136+
ii1, jj1, pf1 = get_i_j_pf(orig_tile)
137+
ii2, jj2, pf2 = get_i_j_pf(mom_tile)
138+
i1 = 0
139+
indices =[]
140+
for i2 in range(len(ii2)):
141+
match = ii1[i1] == ii2[i2] and jj1[i1] == jj2[i2] and pf1[i1] == pf2[i2]
142+
while not match :
143+
i1 = i1 +1
144+
match = ii1[i1] == ii2[i2] and jj1[i1] == jj2[i2] and pf1[i1] == pf2[i2]
145+
indices.append(i1)
146+
i1 = i1 + 1
147+
148+
catch_params_file = bcs_dir + '/' + bcs_ver + '/land/' + air_res + '/clsm/catch_params.nc4'
149+
catchcn_params_file = bcs_dir + '/' + bcs_ver + '/land/' + air_res + '/clsm/catchcn_params.nc4'
150+
mom_catch_params_file = 'clsm/catch_params.nc4'
151+
mom_catchcn_params_file = 'clsm/catchcn_params.nc4'
152+
153+
create_reduced_nc4(catch_params_file, mom_catch_params_file, indices)
154+
create_reduced_nc4(catchcn_params_file, mom_catchcn_params_file, indices)
155+
156+
files = glob.glob(bcs_dir + '/' + bcs_ver + '/land/' + air_res + '/*.da*')
157+
for file in files:
158+
fname = os.path.basename(file)
159+
if 'vegdyn' in fname:
160+
create_reduced_nc4( file,'clsm/vegdyn.data', indices)
161+
else:
162+
shortname = 'clsm/' + fname.split('_')[0] + '.dat'
163+
create_reduced_bin(file, shortname, indices)

‎GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/TileFile_ASCII_to_nc4.F90‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ program TileFile_ASCII_to_nc4
4747

4848
open (newunit=unit, file=trim(tile_file), form='formatted', action='read')
4949

50-
read (unit,*) tmpline ! header line 1: N_tile [maxcat] nx ny (see below)
50+
read (unit,"(A)") tmpline ! header line 1: N_tile [maxcat] nx ny (see below)
5151
read (unit,*) N_grid ! header line 2: N_grid [=1 for EASE, =2 otherwise]
5252
read (unit,*) gName1 ! header line 3: name of atm grid
5353
read (unit,*) n_lon1 ! header line 4: N_lon of atm grid
@@ -76,11 +76,11 @@ program TileFile_ASCII_to_nc4
7676
! in some legacy bcs, dummy ocean grid info is included in header (despite N_grid=1);
7777
! read next line and decide if it is dummy header or info for first tile
7878

79-
read (unit,*) tmpline
79+
read (unit,"(A)") tmpline
8080
if (index(tmpline,'OCEAN')/=0) then
8181
read (unit,*)
8282
read (unit,*)
83-
read (unit,*) tmpline
83+
read (unit,"(A)") tmpline
8484
endif
8585

8686
else
@@ -90,7 +90,7 @@ program TileFile_ASCII_to_nc4
9090
read (unit,*) gName2
9191
read (unit,*) n_lon2
9292
read (unit,*) n_lat2
93-
read (unit,*) tmpline ! read info for first tile (to accommodate legacy EASE grid issues above)
93+
read (unit,"(A)") tmpline ! read info for first tile (to accommodate legacy EASE grid issues above)
9494

9595
endif
9696

‎GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_cube.py‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -41,9 +41,12 @@
4141
bin/CombineRasters.x -f 0 -t {NT} {OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter >/dev/null
4242
bin/CombineRasters.x -t {NT} CF{NC}x6C{SGNAME} {OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter
4343
bin/mk_runofftbl.x -g CF{NC}x6C{SGNAME}_{OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv}
44-
setenv OMP_NUM_THREADS {NCPUS}
45-
if ({SKIPLAND} != True) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C{SGNAME}_{OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv}
46-
setenv OMP_NUM_THREADS 1
44+
45+
if ({SKIPLAND} != True) then
46+
bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C{SGNAME}_{OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} -p no
47+
bin/ExtractBCsFromOrig.py {BCS_DIR} {lbcsv} CF{NC}x6C{SGNAME} {OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}
48+
endif
49+
4750
chmod 755 bin/create_README.csh
4851
bin/create_README.csh
4952
endif
@@ -144,6 +147,7 @@ def make_bcs_cube(config):
144147
SCRATCH_DIR = scratch_dir, \
145148
bin_dir = bin_dir, \
146149
MAKE_BCS_INPUT_DIR = config['inputdir'], \
150+
BCS_DIR = config['bcs_dir'], \
147151
DATENAME = DATENAME, \
148152
POLENAME = POLENAME, \
149153
OCEAN_VERSION = OCEAN_VERSION, \

‎GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,8 @@ def get_configs_from_answers(answers):
4343
else:
4444
make_bcs_input_dir = '/nobackup/gmao_SIteam/ModelData/make_bcs_inputs/'
4545

46+
bcs_dir ='/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/'
47+
4648
user = get_user()
4749
expdir = answers["out_path"]
4850
now = datetime.now()
@@ -148,6 +150,7 @@ def get_configs_from_answers(answers):
148150
config ['expdir'] = expdir
149151
config ['outdir'] = outdir
150152
config ['inputdir'] = make_bcs_input_dir
153+
config ['bcs_dir'] = bcs_dir
151154
config ['NCPUS'] = 16
152155

153156
for x in answers.get('Stretched_CS',[]):

‎GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py‎

Lines changed: 24 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -112,31 +112,30 @@ def get_script_mv(grid_type):
112112
ln -s vegdyn_{RC}.dat vegdyn_{RC}.nc4
113113
cd ../../
114114
115-
/bin/mv clsm/ar.new \\
116-
clsm/bf.dat \\
117-
clsm/ts.dat \\
118-
clsm/catchment.def \\
119-
clsm/cti_stats.dat \\
120-
clsm/tau_param.dat \\
121-
clsm/soil_param.dat \\
122-
clsm/mosaic_veg_typs_fracs \\
123-
clsm/soil_param.first \\
124-
clsm/bad_sat_param.tiles \\
125-
clsm/README \\
126-
clsm/lai.* \\
127-
clsm/AlbMap* \\
128-
clsm/g5fmt \\
129-
clsm/vegetation.hst2 \\
130-
clsm/pfaf_fractions.dat \\
131-
clsm/plots \\
132-
clsm/CLM_veg_typs_fracs \\
133-
clsm/Grid2Catch_TransferData.nc \\
134-
clsm/CLM_NDep_SoilAlb_T2m \\
135-
clsm/CLM4.5_abm_peatf_gdp_hdm_fc \\
136-
clsm/catch_params.nc4 \\
137-
clsm/catchcn_params.nc4 \\
138-
clsm/country_and_state_code.data \\
139-
land/{GRIDNAME}/clsm/
115+
/bin/mv clsm/ar.new land/{GRIDNAME}/clsm/
116+
/bin/mv clsm/bf.dat land/{GRIDNAME}/clsm/
117+
/bin/mv clsm/ts.dat land/{GRIDNAME}/clsm/
118+
/bin/mv clsm/catchment.def land/{GRIDNAME}/clsm/
119+
/bin/mv clsm/cti_stats.dat land/{GRIDNAME}/clsm/
120+
/bin/mv clsm/tau_param.dat land/{GRIDNAME}/clsm/
121+
/bin/mv clsm/soil_param.dat land/{GRIDNAME}/clsm/
122+
/bin/mv clsm/mosaic_veg_typs_fracs land/{GRIDNAME}/clsm/
123+
/bin/mv clsm/soil_param.first land/{GRIDNAME}/clsm/
124+
/bin/mv clsm/bad_sat_param.tiles land/{GRIDNAME}/clsm/
125+
/bin/mv clsm/README land/{GRIDNAME}/clsm/
126+
/bin/mv clsm/lai.* land/{GRIDNAME}/clsm/
127+
/bin/mv clsm/AlbMap* land/{GRIDNAME}/clsm/
128+
/bin/mv clsm/g5fmt land/{GRIDNAME}/clsm/
129+
/bin/mv clsm/vegetation.hst2 land/{GRIDNAME}/clsm/
130+
/bin/mv clsm/pfaf_fractions.dat land/{GRIDNAME}/clsm/
131+
/bin/mv clsm/plots land/{GRIDNAME}/clsm/
132+
/bin/mv clsm/CLM_veg_typs_fracs land/{GRIDNAME}/clsm/
133+
/bin/mv clsm/Grid2Catch_TransferData.nc land/{GRIDNAME}/clsm/
134+
/bin/mv clsm/CLM_NDep_SoilAlb_T2m land/{GRIDNAME}/clsm/
135+
/bin/mv clsm/CLM4.5_abm_peatf_gdp_hdm_fc land/{GRIDNAME}/clsm/
136+
/bin/mv clsm/catch_params.nc4 land/{GRIDNAME}/clsm/
137+
/bin/mv clsm/catchcn_params.nc4 land/{GRIDNAME}/clsm/
138+
/bin/mv clsm/country_and_state_code.data land/{GRIDNAME}/clsm/
140139
141140
"""
142141
mv_template = mv_template + get_change_til_file(grid_type)

‎GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ PROGRAM mkCatchParam
3737
integer :: NC = i_raster, NR = j_raster
3838
character*5 :: LBCSV = 'UNDEF'
3939
character*128 :: Gridname = ''
40+
character*128 :: withbcs = ''
4041
character*128 :: ARG, MaskFile
4142
character*256 :: CMD
4243
character*1 :: opt
@@ -47,7 +48,7 @@ PROGRAM mkCatchParam
4748
integer :: I, J, command_argument_count, nxt
4849
real*8 :: dx, dy, lon0
4950
logical :: regrid
50-
character(len=400), dimension (6) :: Usage
51+
character(len=128), dimension (7) :: Usage
5152
character*128 :: Grid2
5253
character*2 :: poles
5354
character*128 :: fnameRst = '' ! a.k.a. "gfile[r]" in mod_process_hres_data.F90
@@ -115,7 +116,8 @@ PROGRAM mkCatchParam
115116
USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 "
116117
USAGE(4) =" -g: Gridname (name of the .til or .rst file *without* file extension) "
117118
USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC "
118-
USAGE(6) =" -v LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08 v09 ) "
119+
USAGE(6) =" -v: Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08 v09 ) "
120+
USAGE(7) =" -p: if no, it creates catchment_def and nc4 tile files then exits "
119121

120122
! Process Arguments
121123
!------------------
@@ -130,7 +132,7 @@ PROGRAM mkCatchParam
130132
write (log_file,'(a)')trim(cmd)
131133
write (log_file,'(a)')' '
132134
endif
133-
135+
withbcs = 'yes'
134136
I = command_argument_count()
135137
if(I < 1 .or. I > 10) then
136138
write (log_file,'(a)') "Wrong Number of arguments: ", i
@@ -165,6 +167,8 @@ PROGRAM mkCatchParam
165167
call init_bcs_config (trim(LBCSV)) ! get bcs details from version string
166168
case ('b')
167169
DL = trim(arg)
170+
case ('p')
171+
withbcs = trim(arg)
168172
case default
169173
do j = 1,size(usage)
170174
print "(sp,a100)", Usage(j)
@@ -274,7 +278,14 @@ PROGRAM mkCatchParam
274278
endif
275279
write (log_file,'(a)')' '
276280

281+
if (trim(withbcs) == 'no') then
282+
write (log_file,'(a)')'Skipping MOM BCs. BCs will be extracted from the corresponding BCs. '
283+
close (log_file,status='keep')
284+
call exit(0)
285+
endif
286+
277287
call MAPL_ReadTilingNC4( trim(fnameTil)//".nc4", iTable = iTable)
288+
278289
N_land = count(iTable(:,0) == 100) ! n_land = number of land tiles
279290
allocate(tile_j_dum, source = iTable(1:n_land,7)) ! possible used in cti_stats.dat
280291
deallocate (iTable)

0 commit comments

Comments
 (0)