Skip to content

Commit 0bee935

Browse files
kshedstromkshedstrom
committed
Added split_diatom scripts and wrf_to_roms.py
1 parent f02324a commit 0bee935

File tree

5 files changed

+374
-3
lines changed

5 files changed

+374
-3
lines changed

‎Troubles

+9
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
For the future:
2+
3+
/import/home/kshedstrom/python/lib/python3.6/site-packages/pyroms_toolbox/BGrid_GFDL/get_nc_BGrid_GFDL.py:141: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.
4+
Check the NumPy 1.11 release notes for more information.
5+
h[:,0] = h[:,-2]
6+
/import/home/kshedstrom/python/lib/python3.6/site-packages/pyroms_toolbox/BGrid_GFDL/get_nc_BGrid_GFDL.py:145: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.
7+
Check the NumPy 1.11 release notes for more information.
8+
f[:,0] = f[:,-2]
9+

‎examples/cobalt-preproc/Boundary_bio/split_bdry_diatom.py

+147
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
import numpy as np
2+
import netCDF4
3+
import sys
4+
5+
# This program splits the large phytoplankton into large and medium.
6+
# True for nitrogen, iron, and silicon.
7+
# Copy the original file before running this, then operate on the copy.
8+
# Could probably do this with ncap2 instead.
9+
#
10+
ncfile = sys.argv[1]
11+
nc = netCDF4.Dataset(ncfile, 'a', format='NETCDF3_CLASSIC')
12+
spval = -10000000000.
13+
14+
nlg = nc.variables['nlg_north'][:]
15+
nlg *= 0.5
16+
nc.variables['nlg_north'][:] = nlg
17+
18+
nc.createVariable('nmd_north', 'f8', ('ocean_time', 's_rho', 'xi_rho'), fill_value=spval)
19+
nc.variables['nmd_north'].long_name = 'Medium Phytoplankton Nitrogen north boundary condition'
20+
nc.variables['nmd_north'].units = 'mol/kg'
21+
nc.variables['nmd_north'].time = 'ocean_time'
22+
nc.variables['nmd_north'].field = 'nmd_north, scalar, series'
23+
nc.variables['nmd_north'][:] = nlg
24+
25+
felg = nc.variables['felg_north'][:]
26+
felg *= 0.5
27+
nc.variables['felg_north'][:] = felg
28+
29+
nc.createVariable('femd_north', 'f8', ('ocean_time', 's_rho', 'xi_rho'), fill_value=spval)
30+
nc.variables['femd_north'].long_name = 'Medium Phytoplankton Iron north boundary condition'
31+
nc.variables['femd_north'].units = 'mol/kg'
32+
nc.variables['femd_north'].time = 'ocean_time'
33+
nc.variables['femd_north'].field = 'femd_north, scalar, series'
34+
nc.variables['femd_north'][:] = felg
35+
36+
silg = nc.variables['silg_north'][:]
37+
silg *= 0.5
38+
nc.variables['silg_north'][:] = silg
39+
40+
nc.createVariable('simd_north', 'f8', ('ocean_time', 's_rho', 'xi_rho'), fill_value=spval)
41+
nc.variables['simd_north'].long_name = 'Medium Phytoplankton Silicon north boundary condition'
42+
nc.variables['simd_north'].units = 'mol/kg'
43+
nc.variables['simd_north'].time = 'ocean_time'
44+
nc.variables['simd_north'].field = 'simd_north, scalar, series'
45+
nc.variables['simd_north'][:] = silg
46+
47+
48+
nlg = nc.variables['nlg_south'][:]
49+
nlg *= 0.5
50+
nc.variables['nlg_south'][:] = nlg
51+
52+
nc.createVariable('nmd_south', 'f8', ('ocean_time', 's_rho', 'xi_rho'), fill_value=spval)
53+
nc.variables['nmd_south'].long_name = 'Medium Phytoplankton Nitrogen south boundary condition'
54+
nc.variables['nmd_south'].units = 'mol/kg'
55+
nc.variables['nmd_south'].time = 'ocean_time'
56+
nc.variables['nmd_south'].field = 'nmd_south, scalar, series'
57+
nc.variables['nmd_south'][:] = nlg
58+
59+
felg = nc.variables['felg_south'][:]
60+
felg *= 0.5
61+
nc.variables['felg_south'][:] = felg
62+
63+
nc.createVariable('femd_south', 'f8', ('ocean_time', 's_rho', 'xi_rho'), fill_value=spval)
64+
nc.variables['femd_south'].long_name = 'Medium Phytoplankton Iron south boundary condition'
65+
nc.variables['femd_south'].units = 'mol/kg'
66+
nc.variables['femd_south'].time = 'ocean_time'
67+
nc.variables['femd_south'].field = 'femd_south, scalar, series'
68+
nc.variables['femd_south'][:] = felg
69+
70+
silg = nc.variables['silg_south'][:]
71+
silg *= 0.5
72+
nc.variables['silg_south'][:] = silg
73+
74+
nc.createVariable('simd_south', 'f8', ('ocean_time', 's_rho', 'xi_rho'), fill_value=spval)
75+
nc.variables['simd_south'].long_name = 'Medium Phytoplankton Silicon south boundary condition'
76+
nc.variables['simd_south'].units = 'mol/kg'
77+
nc.variables['simd_south'].time = 'ocean_time'
78+
nc.variables['simd_south'].field = 'simd_south, scalar, series'
79+
nc.variables['simd_south'][:] = silg
80+
81+
nlg = nc.variables['nlg_west'][:]
82+
nlg *= 0.5
83+
nc.variables['nlg_west'][:] = nlg
84+
85+
nc.createVariable('nmd_west', 'f8', ('ocean_time', 's_rho', 'eta_rho'), fill_value=spval)
86+
nc.variables['nmd_west'].long_name = 'Medium Phytoplankton Nitrogen west boundary condition'
87+
nc.variables['nmd_west'].units = 'mol/kg'
88+
nc.variables['nmd_west'].time = 'ocean_time'
89+
nc.variables['nmd_west'].field = 'nmd_west, scalar, series'
90+
nc.variables['nmd_west'][:] = nlg
91+
92+
felg = nc.variables['felg_west'][:]
93+
felg *= 0.5
94+
nc.variables['felg_west'][:] = felg
95+
96+
nc.createVariable('femd_west', 'f8', ('ocean_time', 's_rho', 'eta_rho'), fill_value=spval)
97+
nc.variables['femd_west'].long_name = 'Medium Phytoplankton Iron west boundary condition'
98+
nc.variables['femd_west'].units = 'mol/kg'
99+
nc.variables['femd_west'].time = 'ocean_time'
100+
nc.variables['femd_west'].field = 'femd_west, scalar, series'
101+
nc.variables['femd_west'][:] = felg
102+
103+
silg = nc.variables['silg_west'][:]
104+
silg *= 0.5
105+
nc.variables['silg_west'][:] = silg
106+
107+
nc.createVariable('simd_west', 'f8', ('ocean_time', 's_rho', 'eta_rho'), fill_value=spval)
108+
nc.variables['simd_west'].long_name = 'Medium Phytoplankton Silicon west boundary condition'
109+
nc.variables['simd_west'].units = 'mol/kg'
110+
nc.variables['simd_west'].time = 'ocean_time'
111+
nc.variables['simd_west'].field = 'simd_west, scalar, series'
112+
nc.variables['simd_west'][:] = silg
113+
114+
nlg = nc.variables['nlg_east'][:]
115+
nlg *= 0.5
116+
nc.variables['nlg_east'][:] = nlg
117+
118+
nc.createVariable('nmd_east', 'f8', ('ocean_time', 's_rho', 'eta_rho'), fill_value=spval)
119+
nc.variables['nmd_east'].long_name = 'Medium Phytoplankton Nitrogen east boundary condition'
120+
nc.variables['nmd_east'].units = 'mol/kg'
121+
nc.variables['nmd_east'].time = 'ocean_time'
122+
nc.variables['nmd_east'].field = 'nmd_east, scalar, series'
123+
nc.variables['nmd_east'][:] = nlg
124+
125+
felg = nc.variables['felg_east'][:]
126+
felg *= 0.5
127+
nc.variables['felg_east'][:] = felg
128+
129+
nc.createVariable('femd_east', 'f8', ('ocean_time', 's_rho', 'eta_rho'), fill_value=spval)
130+
nc.variables['femd_east'].long_name = 'Medium Phytoplankton Iron east boundary condition'
131+
nc.variables['femd_east'].units = 'mol/kg'
132+
nc.variables['femd_east'].time = 'ocean_time'
133+
nc.variables['femd_east'].field = 'femd_east, scalar, series'
134+
nc.variables['femd_east'][:] = felg
135+
136+
silg = nc.variables['silg_east'][:]
137+
silg *= 0.5
138+
nc.variables['silg_east'][:] = silg
139+
140+
nc.createVariable('simd_east', 'f8', ('ocean_time', 's_rho', 'eta_rho'), fill_value=spval)
141+
nc.variables['simd_east'].long_name = 'Medium Phytoplankton Silicon east boundary condition'
142+
nc.variables['simd_east'].units = 'mol/kg'
143+
nc.variables['simd_east'].time = 'ocean_time'
144+
nc.variables['simd_east'].field = 'simd_east, scalar, series'
145+
nc.variables['simd_east'][:] = silg
146+
147+
nc.close()

‎examples/cobalt-preproc/Initial_bio/README

+7-3
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,8 @@
1-
# first run make_ic_file_bio.py to create initial conditions from COBALT results
2-
3-
# second run make_ic_file_bio_addons.py to create the WOA and GLODAP fields and merge with original file
1+
# First, run make_remap_weights_file.py to create the weights files.
2+
#
3+
# Second, run make_ic_file_bio.py to create initial conditions from COBALT results.
4+
#
5+
# Third, run make_ic_file_bio_addons.py to create the WOA and GLODAP fields and merge with original file.
6+
#
7+
# Fourth, if running with coastal diatoms, run the split_diatom scripts.
48

‎examples/cobalt-preproc/Initial_bio/split_diatom.py

+74
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
import numpy as np
2+
import netCDF4
3+
import sys
4+
5+
# This program splits the large phytoplankton into large and medium.
6+
# True for nitrogen, iron, and silicon.
7+
# Copy the original file before running this, then operate on the copy.
8+
# Could probably do this with ncap2 instead.
9+
#
10+
ncfile = sys.argv[1]
11+
nc = netCDF4.Dataset(ncfile, 'a', format='NETCDF3_CLASSIC')
12+
spval = 1.e+37
13+
14+
nlg = nc.variables['nlg'][:]
15+
nlg *= 0.5
16+
nc.variables['nlg'][:] = nlg
17+
18+
nc.createVariable('nmd', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
19+
#nc.createVariable('nmd', 'f8', ('ocean_time', 'two', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
20+
nc.variables['nmd'].long_name = 'Medium Phytoplankton Nitrogen'
21+
nc.variables['nmd'].units = 'mol/kg'
22+
nc.variables['nmd'].time = 'ocean_time'
23+
nc.variables['nmd'].field = 'nmd, scalar, series'
24+
nc.variables['nmd'].grid = 'grid'
25+
nc.variables['nmd'].location = 'face'
26+
nc.variables['nmd'].coordinates = 'lon_rho lat_rho s_rho ocean_time'
27+
nc.variables['nmd'][:] = nlg
28+
29+
felg = nc.variables['felg'][:]
30+
felg *= 0.5
31+
nc.variables['felg'][:] = felg
32+
33+
nc.createVariable('femd', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
34+
#nc.createVariable('femd', 'f8', ('ocean_time', 'two', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
35+
nc.variables['femd'].long_name = 'Medium Phytoplankton Iron'
36+
nc.variables['femd'].units = 'mol/kg'
37+
nc.variables['femd'].time = 'ocean_time'
38+
nc.variables['femd'].field = 'femd, scalar, series'
39+
nc.variables['femd'].grid = 'grid'
40+
nc.variables['femd'].location = 'face'
41+
nc.variables['femd'].coordinates = 'lon_rho lat_rho s_rho ocean_time'
42+
nc.variables['femd'][:] = felg
43+
44+
silg = nc.variables['silg'][:]
45+
silg *= 0.5
46+
nc.variables['silg'][:] = silg
47+
48+
nc.createVariable('simd', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
49+
#nc.createVariable('simd', 'f8', ('ocean_time', 'two', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
50+
nc.variables['simd'].long_name = 'Medium Phytoplankton Silicon'
51+
nc.variables['simd'].units = 'mol/kg'
52+
nc.variables['simd'].time = 'ocean_time'
53+
nc.variables['simd'].field = 'simd, scalar, series'
54+
nc.variables['simd'].grid = 'grid'
55+
nc.variables['simd'].location = 'face'
56+
nc.variables['simd'].coordinates = 'lon_rho lat_rho s_rho ocean_time'
57+
nc.variables['simd'][:] = silg
58+
59+
mu_mem_lg = nc.variables['mu_mem_lg'][:]
60+
mu_mem_lg *= 0.5
61+
nc.variables['mu_mem_lg'][:] = mu_mem_lg
62+
63+
nc.createVariable('mu_mem_md', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
64+
#nc.createVariable('mu_mem_md', 'f8', ('ocean_time', 'two', 's_rho', 'eta_rho', 'xi_rho'), fill_value=spval)
65+
nc.variables['mu_mem_md'].long_name = 'Medium Phytoplankton Silicon'
66+
nc.variables['mu_mem_md'].units = 'mol/kg'
67+
nc.variables['mu_mem_md'].time = 'ocean_time'
68+
nc.variables['mu_mem_md'].field = 'mu_mem_md, scalar, series'
69+
nc.variables['mu_mem_md'].grid = 'grid'
70+
nc.variables['mu_mem_md'].location = 'face'
71+
nc.variables['mu_mem_md'].coordinates = 'lon_rho lat_rho s_rho ocean_time'
72+
nc.variables['mu_mem_md'][:] = mu_mem_lg
73+
74+
nc.close()

‎examples/stuff/wrf_to_roms.py

+137
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
import numpy as np
2+
import netCDF4
3+
import sys
4+
import os
5+
from mpl_toolkits.basemap import Basemap
6+
import matplotlib.colors as colors
7+
from scipy.signal import medfilt2d
8+
from scipy.interpolate import interp2d
9+
#from scipy.interpolate import griddata
10+
11+
import pyroms
12+
import pyroms_toolbox
13+
from bathy_smoother import *
14+
import creep
15+
16+
fname = "geo_em.d01.nc"
17+
bathyfile = '/archive/u1/uaf/kate/bathy/ARDEM_1.05.nc'
18+
f2name = "Bering_WRF_grid.nc"
19+
ncid = netCDF4.Dataset(fname, "r")
20+
21+
rlat = ncid.variables['XLAT_M'][0,:,:]
22+
rlon = ncid.variables['XLONG_M'][0,:,:]
23+
latv = ncid.variables['XLAT_V'][0,1:-1,:]
24+
lonv = ncid.variables['XLONG_V'][0,1:-1,:]
25+
latu = ncid.variables['XLAT_U'][0,:,1:-1]
26+
lonu = ncid.variables['XLONG_U'][0,:,1:-1]
27+
f = ncid.variables['F'][:]
28+
cosang = ncid.variables['COSALPHA'][:]
29+
sinang = ncid.variables['SINALPHA'][:]
30+
mask = ncid.variables['LANDMASK'][:]
31+
f = ncid.variables['F'][:]
32+
corner_lats = ncid.corner_lats
33+
corner_lons = ncid.corner_lons
34+
35+
map_proj = ncid.MAP_PROJ
36+
if (map_proj == 1):
37+
my_proj = 'lcc'
38+
lat_1 = ncid.TRUELAT1
39+
lat_2 = ncid.TRUELAT2
40+
lon_0 = ncid.STAND_LON
41+
llcrnrlon = corner_lons[12]
42+
llcrnrlat = corner_lats[12]
43+
urcrnrlon = corner_lons[14]
44+
urcrnrlat = corner_lats[14]
45+
46+
ncid.close()
47+
48+
map = Basemap(projection=my_proj, lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, \
49+
llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat, urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat, \
50+
resolution='h')
51+
52+
# define the 4 corners of the grid
53+
# first point is the top left corner then counter clock wise rotation
54+
lon0=corner_lons[13] ; lat0=corner_lats[13]
55+
lon1=corner_lons[12] ; lat1=corner_lats[12]
56+
lon2=corner_lons[15] ; lat2=corner_lats[15]
57+
lon3=corner_lons[14] ; lat3=corner_lats[14]
58+
59+
#generate the new grid
60+
lonp=np.array([lon0, lon1, lon2, lon3])
61+
latp=np.array([lat0, lat1, lat2, lat3])
62+
63+
# shift data so lons go from 0 to 360 instead of -180 to 180.
64+
lonp = np.where(lonp < 0, lonp+360, lonp)
65+
66+
beta = np.array([1, 1, 1, 1])
67+
68+
Mp, Lp = rlon.shape
69+
hgrd = pyroms.grid.Gridgen(lonp, latp, beta, (Mp+1,Lp+1), proj=map)
70+
71+
lonv, latv = map(hgrd.x_vert, hgrd.y_vert, inverse=True)
72+
hgrd = pyroms.grid.CGrid_geo(lonv, latv, map)
73+
hgrd.lon_rho = np.where(hgrd.lon_rho < 0, hgrd.lon_rho+360, hgrd.lon_rho)
74+
75+
hgrd.mask_rho = (1.0 - mask)
76+
77+
ncid = netCDF4.Dataset(bathyfile, "r")
78+
79+
lons = ncid.variables['lon'][:]
80+
lats = ncid.variables['lat'][:]
81+
topo = ncid.variables['z'][:]
82+
ncid.close()
83+
84+
# depth positive
85+
topo = -topo
86+
87+
# fix minimum depth
88+
hmin = 5
89+
topo = pyroms_toolbox.change(topo, '<', hmin, hmin)
90+
91+
# interpolate new bathymetry
92+
lon, lat = meshgrid(lons, lats)
93+
h = griddata(lon.flat,lat.flat,topo.flat,hgrd.lon_rho,hgrd.lat_rho)
94+
#h = griddata((lon,lat),topo.flat,(hgrd.lon_rho,hgrd.lat_rho),method='nearest')
95+
#interpf = interp2d(lons,lats,topo)
96+
#lon_rho = reshape(hgrd.lon_rho,Lp*Mp)
97+
#lat_rho = reshape(hgrd.lat_rho,Lp*Mp)
98+
#h = interpf(lon_rho,lat_rho)
99+
100+
# insure that depth is always deeper than hmin
101+
h = pyroms_toolbox.change(h, '<', hmin, hmin)
102+
103+
# check bathymetry roughness
104+
hgrd.mask_rho = reshape(hgrd.mask_rho, (Mp,Lp))
105+
RoughMat = bathy_tools.RoughnessMatrix(h, hgrd.mask_rho)
106+
print 'Max Roughness value is: ', RoughMat.max()
107+
#h = creep.cslf(h, nan, -200., 200.)
108+
h = np.where(isnan(h), 5500.0, h)
109+
hraw = h.copy()
110+
111+
# smooth the raw bathy using the direct iterative method from Martinho and Batteen
112+
# (2006)
113+
rx0_max = 0.35
114+
h = bathy_smoothing.smoothing_Positive_rx0(hgrd.mask_rho, h, rx0_max)
115+
116+
# check bathymetry roughness again
117+
RoughMat = bathy_tools.RoughnessMatrix(h, hgrd.mask_rho)
118+
print 'Max Roughness value is: ', RoughMat.max()
119+
h = pyroms_toolbox.shapiro_filter.shapiro2(h, 2)
120+
121+
hgrd.h = h
122+
# define vertical grd
123+
Vtrans = 2
124+
theta_s = 7.0
125+
theta_b = 0.1
126+
Tcline = 250
127+
N = 50
128+
vgrd = pyroms.vgrid.s_coordinate_4(h, theta_b, theta_s, Tcline, N, hraw=hraw)
129+
vgrd.h = h # what the hell??
130+
131+
#ROMS grid
132+
grd_name = 'Bering_WRF'
133+
grd = pyroms.grid.ROMS_Grid(grd_name, hgrd, vgrd)
134+
135+
#write grid to netcdf file
136+
pyroms.grid.write_ROMS_grid(grd, filename=f2name)
137+

0 commit comments

Comments
 (0)