forked from hmchristensen/MUMIP-cg
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathfunc_geostrophic.ncl
181 lines (150 loc) · 5.93 KB
/
func_geostrophic.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
; func_geostrophic.ncl
function func_geostrophic(fout:file)
local fout,tv,tv_top,tv_bot,grav,R_dry,rpi,rday,rsiyea,rsiday,\
romega,f_cori_vec,f_cori,smooth_me,\
logpres,gradp,gradp_x,gradp_y,press_gradx_out,press_grady_out,vg_out,ug_out,\
dsizes,dsizes2,tmparray,lsm_tmp,lsm,lsm_perm,d_size,pgx_out_lsm,pgy_out_lsm,\
is_cyclic,guess_type,nscan,eps,relc, opt, rho_sfc
begin
;==================================================================================
; "geostrophic winds" : calculate forcing term from pressure gradient
;
; NB this file has been simplified cf. Cascade file as we are only considering ocean points
; => our partial derivatives along model levels are at constant height
;
; Take care if we want to move to include land points. Because ICON is not a true sigma co-ordinate system,
; we lose some of the nice properties of derivatives along model levels that are used in
; ECMWF documentation
print("procedure geostrophic winds")
;=====================
; load required variables in from file
pres_data_out = fout->pressure_t
lsm = fout->lsm
t_data_out = fout->temp_t
qv_data_out = fout->qv_t
ql_data_out = fout->ql_t
qi_data_out = fout->qi_t
; ps_data_out = fout->ps
; z_sfc = fout->z_sfc
lev_out = t_data_out&lev
time_out = t_data_out&time
n_lev_out = dimsizes(lev_out)
;=====================
; do calculation
; set up constants - following IFS convention
; grav = 9.80665
R_dry = 287.0
rpi = 4*atan(1.0);
rday = 86400;
rsiyea = 365.25*rday*2.0*rpi/6.283076;
rsiday=rday/(1.0+rday/rsiyea);
romega=2.0*rpi/rsiday;
f_cori_vec=2*romega*sin(pres_data_out&lat*rpi/180);
f_cori=conform_dims(dimsizes(pres_data_out),f_cori_vec,2)
delete([/f_cori_vec/])
;=====================
; virtual temperature
; p. 17, Rogers and Yau
qt_data = qv_data_out+ql_data_out+qi_data_out
; deduce specific mass of dry air
qd_data = 1 - qt_data
; convert to m.r.
rv_data_out = qv_data_out/qd_data
epsilon = 0.622
tv_top = t_data_out*(1 + rv_data_out/epsilon)
tv_bot = 1 + rv_data_out
tv = tv_top/tv_bot
copy_VarCoords(t_data_out,tv)
;----
; (1/rho) * partial d p by dx at constant eta
logpres = log(pres_data_out)
; HMC remove smoothing
; logpres = smth9(logpres,0.5,0.25,False) ; heavy smoothing
copy_VarCoords(pres_data_out,logpres)
copy_VarAtts(pres_data_out,logpres)
smooth_me = False
gradp = hor_derivative(logpres,smooth_me)
gradp_x = gradp[0]
gradp_y = gradp[1]
delete([/gradp,logpres/])
press_gradx_out = R_dry*gradp_x*tv
press_grady_out = R_dry*gradp_y*tv
;-- save estimate of rho_sfc for ustar computation
rho_sfc = pres_data_out(:,n_lev_out-1,:,:)/(R_dry*tv(:,n_lev_out-1,:,:))
rho_sfc!0="time"
rho_sfc&time=pres_data_out&time
rho_sfc!1="lat"
rho_sfc&lat=pres_data_out&lat
rho_sfc!2="lon"
rho_sfc&lon=pres_data_out&lon
delete([/tv,tv_top,tv_bot/])
delete([/R_dry,rpi,rday,rsiyea,rsiday,romega/])
delete([/gradp_x,gradp_y/])
;; === Now we use ICON consistent Z-sfc, no need for this masking. ===
; ;=====================================
; ; over land mask out geostrophic winds and
; ; interpolate pressure gradients used to calc ug_out and vg_out
;
; dsizes = dimsizes(lsm)
; dsizes2 = (/dsizes(0),dsizes(1),dimsizes(time_out),dimsizes(lev_out)/)
; tmparray = new(dsizes2,typeof(lsm))
; delete([/dsizes,dsizes2/])
;
; lsm_tmp = conform(tmparray,lsm,(/0,1/))
; lsm_tmp!0 = "lat"
; lsm_tmp!1 = "lon"
; lsm_tmp!2 = "time"
; lsm_tmp!3 = "lev"
; lsm_perm = lsm_tmp(time | :, lev | :, lat | :, lon | :)
; delete([/lsm_tmp/])
;
; ; extend influence of land out over sea
; lsm_perm = smth9(lsm_perm,0.5,0.25,False) ; heavy smoothing
; lsm_perm = smth9(lsm_perm,0.5,0.25,False) ; heavy smoothing
; lsm_perm = smth9(lsm_perm,0.5,0.25,False) ; heavy smoothing
; d_size = dimsizes(lsm_perm)
; lsm_perm(:,:,:,0) = 1 ; mask edge:
; lsm_perm(:,:,0,:) = 1 ; set values to missing
; lsm_perm(:,:,:,d_size(3)-1) = 1 ;
; lsm_perm(:,:,d_size(2)-1,:) = 1 ;
; pgx_out_lsm = where(lsm_perm.ge.0.0003,press_gradx_out@_FillValue,press_gradx_out)
; pgy_out_lsm = where(lsm_perm.ge.0.0003,press_grady_out@_FillValue,press_grady_out)
; copy_VarCoords(press_gradx_out,pgx_out_lsm)
; copy_VarCoords(press_grady_out,pgy_out_lsm)
; copy_VarAtts(press_gradx_out,pgx_out_lsm)
; copy_VarAtts(press_grady_out,pgy_out_lsm)
; delete([/lsm_perm,press_gradx_out,press_grady_out/])
;
;
; ;======================
; ; use poisson grid fill to smoothly fill in missing values
;
; is_cyclic = False ; not cyclic data
; guess_type = 1 ; start with zonal means
; nscan = 200 ; no. iterations
; eps = 1.e-6 ; tolerance
; relc = 0.6 ; relaxation const
; opt = 0 ; dummy
;
; poisson_grid_fill(pgx_out_lsm,is_cyclic,guess_type,nscan,eps,relc,opt)
; poisson_grid_fill(pgy_out_lsm,is_cyclic,guess_type,nscan,eps,relc,opt)
;
; vg_out = pgx_out_lsm/f_cori
; ug_out = -pgy_out_lsm/f_cori
vg_out = press_gradx_out/f_cori
ug_out = -press_grady_out/f_cori
vg_out_flt = tofloat(vg_out)
ug_out_flt = tofloat(ug_out)
copy_VarCoords(pres_data_out,ug_out_flt)
copy_VarCoords(pres_data_out,vg_out_flt)
copy_VarAtts(pres_data_out,ug_out_flt)
copy_VarAtts(pres_data_out,vg_out_flt)
ug_out_flt@long_name = "Geostrophic zonal wind"
ug_out_flt@units = "m/s"
vg_out_flt@long_name = "Geostrophic meridional wind"
vg_out_flt@units = "m/s"
add_to_file(fout,ug_out_flt , "ug")
add_to_file(fout,vg_out_flt , "vg")
delete([/ug_out,vg_out,ug_out_flt,vg_out_flt,f_cori/])
return(rho_sfc)
end