Skip to content

Commit

Permalink
fix rs
Browse files Browse the repository at this point in the history
  • Loading branch information
@Kang-Sun-UB
Kang-Sun-UB committed Dec 16, 2024
1 parent f77625f commit 30a2268
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 17 deletions.
25 changes: 20 additions & 5 deletions CAREER/tempo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -928,15 +928,30 @@ def latlon2m(lat1,lon1,lat2,lon2,m_per_lat,m_per_lon):
lat2=self['latc'][:-2,:-2],lon2=self['lonc'][:-2,:-2],
m_per_lat=self['m_per_lat'],m_per_lon=self['m_per_lon'][1:-1,1:-1]
)
self[f'd{field}dx'] = dfdx

self[f'd{field}dx'] = dfdx
self[f'd{field}dy'] = dfdy
self[f'd{field}dr'] = dfdr
self[f'd{field}ds'] = dfds

### proven to be equivalent to the simple formula below
# # gradient (vec{g}) = gx*vec{x}+gy*vec{y} = gr*vec{r}+gy*vec{r}
# gx = (dfdx-dfdy*self['xdoty'])/(1-np.square(self['xdoty']))
# gy = (dfdy-dfdx*self['xdoty'])/(1-np.square(self['xdoty']))

# gr = (dfdr-dfds*self['rdots'])/(1-np.square(self['rdots']))
# gs = (dfds-dfdr*self['rdots'])/(1-np.square(self['rdots']))

# # directional derivatives
# self[field+'_DD_xy'] = windx*gx + windy*gy + \
# self['xdoty']*(windx*gy + windy*gx)
# self[field+'_DD_rs'] = windr*gr + winds*gs + \
# self['rdots']*(windr*gs + winds*gr)
### proven to be equivalent to the simple formula below

# directional derivatives
self[field+'_DD_xy'] = windx*dfdx + windy*dfdy + \
self['xdoty']*(windx*dfdy + windy*dfdx)
self[field+'_DD_rs'] = windr*dfdr + winds*dfds + \
self['rdots']*(windr*dfds + winds*dfdr)
self[field+'_DD_xy'] = windx*dfdx + windy*dfdy
self[field+'_DD_rs'] = windr*dfdr + winds*dfds
if keep_single_xyrs:
with warnings.catch_warnings():
warnings.filterwarnings(action='ignore', message='Mean of empty slice')
Expand Down
19 changes: 7 additions & 12 deletions popy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1983,8 +1983,8 @@ def F_grads(c,dy,dx_vec,dd_vec,finite_difference_order):
### grad(vcd) dot wind
dcdx,dcdy,dcdr,dcds = F_grads(vcd,dy,dx_vec,dd_vec,finite_difference_order)
wind_column_xy = dcdx*self['wind_e'] + dcdy*self['wind_n']
wind_column_rs = dcdr*self['wind_ne'] + dcds*self['wind_nw'] \
+ r_dot_s*(dcdr*self['wind_nw'] + dcds*self['wind_ne'])
wind_column_rs = dcdr*self['wind_ne'] + dcds*self['wind_nw']

with warnings.catch_warnings():
warnings.filterwarnings(action='ignore', message='Mean of empty slice')
wind_column = np.nanmean(np.array([wind_column_xy,wind_column_rs]),axis=0)
Expand All @@ -2006,8 +2006,7 @@ def F_grads(c,dy,dx_vec,dd_vec,finite_difference_order):
if z0 is not None:
dcdx,dcdy,dcdr,dcds = F_grads(z0,dy,dx_vec,dd_vec,finite_difference_order)
wind_topo_xy = vcd*(dcdx*self['wind_e'] + dcdy*self['wind_n'])
wind_topo_rs = vcd*(dcdr*self['wind_ne'] + dcds*self['wind_nw'] \
+ r_dot_s*(dcdr*self['wind_nw'] + dcds*self['wind_ne']))
wind_topo_rs = vcd*(dcdr*self['wind_ne'] + dcds*self['wind_nw'])
with warnings.catch_warnings():
warnings.filterwarnings(action='ignore', message='Mean of empty slice')
wind_topo = np.nanmean(np.array([wind_topo_xy,wind_topo_rs]),axis=0)
Expand Down Expand Up @@ -2043,10 +2042,8 @@ def F_grads(c,dy,dx_vec,dd_vec,finite_difference_order):
np.power(a0,order)*(dpdx*self['wind_e'] + dpdy*self['wind_n'])

wind_albedo_rs = order*pa*np.power(a0,order-1)*\
(dcdr*self['wind_ne'] + dcds*self['wind_nw']+\
r_dot_s*(dcdr*self['wind_nw'] + dcds*self['wind_ne']))+\
np.power(a0,order)*(dpdr*self['wind_ne'] + dpds*self['wind_nw']+\
r_dot_s*(dpdr*self['wind_nw'] + dpds*self['wind_ne']))
(dcdr*self['wind_ne'] + dcds*self['wind_nw'])+\
np.power(a0,order)*(dpdr*self['wind_ne'] + dpds*self['wind_nw'])

wind_albedo = np.nanmean(np.array([wind_albedo_xy,wind_albedo_rs]),axis=0)
wind_albedo[np.isnan(wind_albedo_xy) & np.isnan(wind_albedo_rs)] = np.nan
Expand All @@ -2067,10 +2064,8 @@ def F_grads(c,dy,dx_vec,dd_vec,finite_difference_order):
np.power(a0,order)*(dpdx*self['wind_e'] + dpdy*self['wind_n'])

wind_albedo_rs = order*pa*np.power(a0,order-1)*\
(dcdr*self['wind_ne'] + dcds*self['wind_nw']+\
r_dot_s*(dcdr*self['wind_nw'] + dcds*self['wind_ne']))+\
np.power(a0,order)*(dpdr*self['wind_ne'] + dpds*self['wind_nw']+\
r_dot_s*(dpdr*self['wind_nw'] + dpds*self['wind_ne']))
(dcdr*self['wind_ne'] + dcds*self['wind_nw'])+\
np.power(a0,order)*(dpdr*self['wind_ne'] + dpds*self['wind_nw'])

wind_albedo = np.nanmean(np.array([wind_albedo_xy,wind_albedo_rs]),axis=0)
wind_albedo[np.isnan(wind_albedo_xy) & np.isnan(wind_albedo_rs)] = np.nan
Expand Down

0 comments on commit 30a2268

Please sign in to comment.