gabriel/Test
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Update calc_ref_state.py and mycalc_ocean.py to calculate z-dep ocean area

Browse Source
This commit is contained in:
Gabriel Wolf 2018-09-21 10:45:42 +01:00
parent 822aafa0aa
commit f2103892c4
2 changed files with 44 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
calc_ref_state.py
View File

@ -109,6 +109,10 @@ except:
grd.p = plev
grd.Up = 'dbar'
grd.Np = len(plev)
if grd.lon[-1]+grd.lon[0] >= 1:
grd.symlon = 1
else:
grd.symlon = 0
lat_xy, lon_xy = np.meshgrid(grd.lat,grd.lon)
p_yz, lat_yz = np.meshgrid(grd.z, grd.lat)
lat_xyz, lon_xyz, p_xyz = np.meshgrid(grd.lat, grd.lon, grd.p)

51
mycalc_ocean.py
View File

@ -2,7 +2,7 @@
# calc_rho_from_theta
# calc_theta_from_temp
# gsw_g
def calc_area_lat(grd,r_earth):
def calc_area_lat(grd,r_earth,calc_old=0):
# def calc_area_lat(grd,r_earth) #########################################################
# written by : Gabriel Wolf, g.a.wolf@reading.ac.uk
# adapted from get_area_latitude.m of Remi Tailleux (09/2018)
@ -19,16 +19,45 @@ def calc_area_lat(grd,r_earth):
# convert deg to rad for sin/cos
deg_to_rad = np.pi/180.0
lat_radians = grd.lat*deg_to_rad
# compute length of longitudinal and latitudinal elements
print 'area calculation not valid for non-equidistant input lon-grid'
print 'area calculation not valid for non-equidistant input lat-grid'
dx = r_earth*deg_to_rad*(grd.Nlon/360.0)
dy = r_earth*deg_to_rad*abs(grd.lat[1]-grd.lat[0])
# calculate area
area_xyz = np.zeros([grd.Nlon,grd.Nlat,grd.Nz])
for i_lon in range(0,grd.Nlon):
for i_z in range(0,grd.Nz):
area_xyz[i_lon,:,i_z] = (dx*dy)*np.cos(lat_radians)
if calc_old == 1:
print '*** old version of area calculation ***'
# compute length of longitudinal and latitudinal elements
print '*** area calculation not valid for non-equidistant input lon-grid ***'
print '*** area calculation not valid for non-equidistant input lat-grid ***'
dx = r_earth*deg_to_rad*(grd.Nlon/360.0)
dy = r_earth*deg_to_rad*abs(grd.lat[1]-grd.lat[0])
# calculate area
area_xyz = np.zeros([grd.Nlon,grd.Nlat,grd.Nz])
for i_lon in range(0,grd.Nlon):
for i_z in range(0,grd.Nz):
area_xyz[i_lon,:,i_z] = (dx*dy)*np.cos(lat_radians)
else:
lat_d, lon_d, z_d = np.meshgrid(grd.lat,grd.lon,grd.z)
if grd.symlon == 1:
delta_lon = 360.0
index_d = np.append(np.arange(1,grd.Nlon),0)
lon_e = 0.5 * (lon_d + lon_d[index_d,:,:])
lon_e[-1,:,:] = lon_e[-1,:,:] + 0.5 * delta_lon
index_d = np.append(grd.Nlon-1,np.arange(0,grd.Nlon-1))
lon_w = 0.5 * (lon_d + lon_d[index_d,:,:])
lon_w[0,:,:] = lon_w[0,:,:] - 0.5 * delta_lon
else:
lon_d = np.append(2*grd.lon[0] - grd.lon[1],grd.lon)
lon_d = np.append(lon_d,2*grd.lon[-1] - grd.lon[-2])
dummy1, lon_d, dummy2 = np.meshgrid(grd.lat,lon_d,grd.z)
index_d = np.arange(1,grd.Nlon+1)
lon_e = 0.5 * (lon_d[index_d,:,:] + delta_lon + lon_d[index_d+1,:,:])
lon_w = 0.5 * (lon_d[index_d,:,:] + delta_lon + lon_d[index_d-1,:,:])
lat_d = np.append(max(2*grd.lat[0] - grd.lat[1],-180.0-grd.lat[0]),grd.lat)
lat_d = np.append(lat_d,min(2*grd.lat[-1] - grd.lat[-2],180.0-grd.lat[-1]))
lat_d, dummy1, dummy2 = np.meshgrid(lat_d,grd.lon,grd.z)
index_d = np.arange(1,grd.Nlat+1)
lat_n = 0.5 * (lat_d[:,index_d,:] + lat_d[:,index_d+1,:])
lat_s = 0.5 * (lat_d[:,index_d,:] + lat_d[:,index_d-1,:])
r_z = r_earth - abs(grd.z)
#area_xyz = (deg_to_rad*r_z)**2*(lon_e-lon_w)*(lat_n-lat_s)*np.cos(lat_d[:,index_d,:]*deg_to_rad)
area_xyz = deg_to_rad*r_z**2*(lon_e-lon_w)*(np.sin(lat_n*deg_to_rad)-np.sin(lat_s*deg_to_rad))
del delta_lon, lon_d, lat_d, z_d, dummy1, dummy2, index_d, r_z
# return horizontal area for input grid
return area_xyz