| 1 |
function [mystruct_out]=gmaze_bulk_coare(mystruct_in,varargin); |
| 2 |
|
| 3 |
%varargin [hu ht hq] |
| 4 |
% [hu ht hq],albedo |
| 5 |
% [albedo HAS AN impact here] |
| 6 |
|
| 7 |
%version with shortened iteration; modified Rt and Rq |
| 8 |
%uses wave information wave period in s and wave ht in m |
| 9 |
%no wave, standard coare 2.6 charnock: jwave=0 |
| 10 |
%Oost et al. zo=50/2/pi L (u*/c)^4.5 if jwave=1 |
| 11 |
%taylor and yelland zo=1200 h*(L/h)^4.5 jwave=2 |
| 12 |
|
| 13 |
%inputs: |
| 14 |
atemp_in=mystruct_in.atemp; |
| 15 |
aqh_in=mystruct_in.aqh; |
| 16 |
uwind_in=mystruct_in.uwind; |
| 17 |
vwind_in=mystruct_in.vwind; |
| 18 |
swdn_in=mystruct_in.swdn; |
| 19 |
lwdn_in=mystruct_in.lwdn; |
| 20 |
slp_in=mystruct_in.slp; |
| 21 |
precip_in=mystruct_in.precip; |
| 22 |
sst_in=mystruct_in.sst; |
| 23 |
|
| 24 |
%compute wind speed: |
| 25 |
wspeed=uwind_in.*uwind_in+vwind_in.*vwind_in; |
| 26 |
wspeed=sqrt(wspeed); %wspeed=max(wspeed,0.5); |
| 27 |
|
| 28 |
%compute bulk water spec hum: |
| 29 |
nlon=size(wspeed,1); nlat=size(wspeed,2); |
| 30 |
Qs=reshape(qsee([reshape(sst_in,[nlon*nlat 1]) reshape(slp_in,[nlon*nlat 1])]),[nlon nlat]); |
| 31 |
% => gives qsat in g/kg |
| 32 |
%in ... exf_bulk_largeyeager04: |
| 33 |
%cvapor_fac=640380; |
| 34 |
%cvapor_exp=5107.4; |
| 35 |
%saltsat=0.980; |
| 36 |
%atmrho=1.2; |
| 37 |
%Qs = cvapor_fac*exp(-cvapor_exp./(sst_in+273.16))*saltsat/atmrho; |
| 38 |
% => gives qsat in kg/kg |
| 39 |
|
| 40 |
|
| 41 |
%change name of variable: |
| 42 |
u = wspeed; % wind speed (m/s) at height zu (m) |
| 43 |
us = 0*ones(size(u)); % surface current speed in the wind direction (m/s) |
| 44 |
ts = sst_in; % bulk water temperature (C) if jcool = 1, interface water T if jcool = 0 |
| 45 |
t = atemp_in-273.16; % bulk air temperature (C), height zt |
| 46 |
Qs = Qs/1000;%bulk water spec hum out of qsee -> /1000 => kg/kg |
| 47 |
Q = aqh_in; %bulk air spec hum is in kg/kg (from input) |
| 48 |
Rs = swdn_in; % downward solar flux (W/m^2) |
| 49 |
Rl = lwdn_in; % downard IR flux (W/m^2) |
| 50 |
rain = precip_in/3600; % rain rate (mm/hr) |
| 51 |
zi = 600*ones(size(u)); % PBL depth (m) |
| 52 |
P = slp_in; % Atmos surface pressure (mb) |
| 53 |
grav = 9.81*ones(size(u)); % grv(lat) |
| 54 |
if nargin==4 |
| 55 |
albed = varargin{2}*ones(size(u)); % albed |
| 56 |
else |
| 57 |
albed = 0.1*ones(size(u)); % albed |
| 58 |
end |
| 59 |
|
| 60 |
% Vectorization: |
| 61 |
[ny nx] = size(u); |
| 62 |
u = u(:); |
| 63 |
us = us(:); |
| 64 |
ts = ts(:); |
| 65 |
t = t(:); |
| 66 |
Qs = Qs(:); |
| 67 |
Q = Q(:); |
| 68 |
Rs = Rs(:); |
| 69 |
Rl = Rl(:); |
| 70 |
rain = rain(:); |
| 71 |
zi = zi(:); |
| 72 |
P = P(:); |
| 73 |
grav = grav(:); |
| 74 |
albed = albed(:); |
| 75 |
|
| 76 |
% 1D fields: |
| 77 |
if nargin==1 |
| 78 |
zu = 10; % wind speed measurement height (m) |
| 79 |
zt = 10; % air T measurement height (m) |
| 80 |
zq = 10; % air q measurement height (m) (specific humidity in g/kg) |
| 81 |
else |
| 82 |
zu=varargin{1}(1); |
| 83 |
zt=varargin{1}(2); |
| 84 |
zq=varargin{1}(3); |
| 85 |
end |
| 86 |
jcool = 1; % implement cool calculation skin switch, 0 = no, 1 = yes |
| 87 |
jwave = 0; % implement wave dependent roughness model |
| 88 |
twave = 0; % wave period (s) |
| 89 |
hwave = 0; % wave height (m) |
| 90 |
% add the adjusted height |
| 91 |
zsu = 10; %wind speed output height (m) |
| 92 |
zst = 10; %air T output height (m) |
| 93 |
zsq = 10; %air q output height (m) |
| 94 |
|
| 95 |
|
| 96 |
% BEGIN THE COMPUTATION: |
| 97 |
|
| 98 |
|
| 99 |
%******.*.*.*.*.* set constants .*.*.*.*.*.*.*.*.*.*.*.*.* |
| 100 |
Beta = 1.2; |
| 101 |
von = 0.4; |
| 102 |
fdg = 1.00; |
| 103 |
tdk = 273.16; |
| 104 |
Rgas = 287.1; |
| 105 |
cpa = 1004.67; |
| 106 |
be = 0.026; |
| 107 |
cpw = 4000; |
| 108 |
rhow = 1022; |
| 109 |
visw = 1e-6; |
| 110 |
tcw = 0.6; |
| 111 |
dter = 0.3; |
| 112 |
|
| 113 |
%.*.*.*.*.*.*.*.*.*.*.*.*.* air constants .*.*.*.*.*.*.*.*.*.*.*.* |
| 114 |
Rgas = 287.1; |
| 115 |
Le = (2.501-.00237.*ts).*1e6; |
| 116 |
cpa = 1004.67; |
| 117 |
cpv = cpa.*(1+0.84.*Q); |
| 118 |
rhoa = P.*100./(Rgas.*(t+tdk).*(1+0.61.*Q)); |
| 119 |
visa = 1.326e-5.*(1+6.542e-3.*t+8.301e-6.*t.*t-4.84e-9.*t.*t.*t); |
| 120 |
|
| 121 |
%.*.*.*.*.*.*.*.*.*.*.*.* cool skin constants .*.*.*.*.*.*.* |
| 122 |
Al = 2.1e-5.*(ts+3.2).^0.79; |
| 123 |
be = 0.026; |
| 124 |
cpw = 4000; |
| 125 |
rhow = 1022; |
| 126 |
visw = 1e-6; |
| 127 |
tcw = 0.6; |
| 128 |
bigc = 16.*grav.*cpw.*(rhow.*visw).^3./(tcw.*tcw.*rhoa.*rhoa); |
| 129 |
wetc = 0.622.*Le.*Qs./(Rgas.*(ts+tdk).^2); |
| 130 |
|
| 131 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.* wave parameters .*.*.*.*.*.*.*.*.* |
| 132 |
lwave=grav./2./pi.*twave.^2; |
| 133 |
cwave=grav./2./pi.*twave; |
| 134 |
|
| 135 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.* compute aux stuff .*.*.*.*.*.*.* |
| 136 |
% Net short wave |
| 137 |
Rns = Rs.*(1-albed); |
| 138 |
% Net long wave |
| 139 |
Rnl = 0.97.*(5.67e-8.*(ts-0.3.*jcool+tdk).^4-Rl); |
| 140 |
|
| 141 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.* Begin bulk loop .*.*.*.*.*.*.* |
| 142 |
|
| 143 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.* first guess .*.*.*.*.*.*.*.*.*.*.*.* |
| 144 |
du=u-us; |
| 145 |
dt=ts-t-.0098.*zt; |
| 146 |
dq=Qs-Q; |
| 147 |
ta=t+tdk; |
| 148 |
ug=.5; |
| 149 |
dter=0.3; |
| 150 |
dqer=wetc.*dter; |
| 151 |
ut=sqrt(du.*du+ug.*ug); |
| 152 |
% whoi using zsu (could be 10 or 2m)instad of 10m |
| 153 |
%u10=ut.*log(10./1e-4)./log(zu./1e-4); |
| 154 |
u10=ut.*log(zsu./1e-4)./log(zu./1e-4); |
| 155 |
usr=.035.*u10; |
| 156 |
zo10=0.011.*usr.*usr./grav+0.11.*visa./usr; |
| 157 |
Cd10=(von./log(zsu./zo10)).^2; |
| 158 |
Ch10=0.00115; |
| 159 |
Ct10=Ch10./sqrt(Cd10); |
| 160 |
zot10=zsu./exp(von./Ct10); |
| 161 |
Cd=(von./log(zu./zo10)).^2; |
| 162 |
Ct=von./log(zt./zot10); |
| 163 |
CC=von.*Ct./Cd; |
| 164 |
Ribcu=-zu./zi./.004./Beta.^3; |
| 165 |
Ribu=-grav.*zu./ta.*((dt-dter.*jcool)+.61.*ta.*dq)./ut.^2; |
| 166 |
nits=3; |
| 167 |
if Ribu<0; |
| 168 |
zetu=CC.*Ribu./(1+Ribu./Ribcu); |
| 169 |
else; |
| 170 |
zetu=CC.*Ribu.*(1+27./9.*Ribu./CC); |
| 171 |
end; |
| 172 |
L10=zu./zetu; |
| 173 |
if zetu>50; |
| 174 |
nits=1; |
| 175 |
end; |
| 176 |
usr=ut.*von./(log(zu./zo10)-psiu_30(zu./L10)); |
| 177 |
tsr=-(dt-dter.*jcool).*von.*fdg./(log(zt./zot10)-psit_30(zt./L10)); |
| 178 |
qsr=-(dq-wetc.*dter.*jcool).*von.*fdg./(log(zq./zot10)-psit_30(zq./L10)); |
| 179 |
|
| 180 |
tkt=.001; |
| 181 |
|
| 182 |
charn=0.011; |
| 183 |
if ut>10 |
| 184 |
charn=0.011+(ut-10)./(18-10).*(0.018-0.011); |
| 185 |
end; |
| 186 |
if ut>18 |
| 187 |
charn=0.018; |
| 188 |
end; |
| 189 |
|
| 190 |
%disp(usr) |
| 191 |
|
| 192 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.* bulk loop .*.*.*.*.*.*.*.*.*.*.*.* |
| 193 |
for i=1:nits; |
| 194 |
|
| 195 |
zet=von.*grav.*zu./ta.*(tsr.*(1+0.61.*Q)+.61.*ta.*qsr)./(usr.*usr)./(1+0.61.*Q); |
| 196 |
if jwave==0;zo=charn.*usr.*usr./grav+0.11.*visa./usr;end; |
| 197 |
if jwave==1;zo=50./2./pi.*lwave.*(usr./cwave).^4.5+0.11.*visa./usr;end;%Oost et al |
| 198 |
if jwave==2;zo=1200.*hwave.*(hwave./lwave).^4.5+0.11.*visa./usr;end;%Taylor and Yelland |
| 199 |
rr=zo.*usr./visa; % roughness Reynolds number |
| 200 |
L=zu./zet; |
| 201 |
zoq=min(1.15e-4,5.5e-5./rr.^.6); |
| 202 |
zot=zoq; |
| 203 |
usr=ut.*von./(log(zu./zo)-psiu_30(zu./L)); |
| 204 |
tsr=-(dt-dter.*jcool).*von.*fdg./(log(zt./zot)-psit_30(zt./L)); |
| 205 |
qsr=-(dq-wetc.*dter.*jcool).*von.*fdg./(log(zq./zoq)-psit_30(zq./L)); |
| 206 |
Bf=-grav./ta.*usr.*(tsr+.61.*ta.*qsr); |
| 207 |
if Bf>0 |
| 208 |
ug=Beta.*(Bf.*zi).^.333; |
| 209 |
else |
| 210 |
ug=.2; |
| 211 |
end; |
| 212 |
ut=sqrt(du.*du+ug.*ug); |
| 213 |
Rnl=0.97.*(5.67e-8.*(ts-dter.*jcool+tdk).^4-Rl); |
| 214 |
hsb=-rhoa.*cpa.*usr.*tsr; |
| 215 |
hlb=-rhoa.*Le.*usr.*qsr; |
| 216 |
qout=Rnl+hsb+hlb; |
| 217 |
dels=Rns.*(.065+11.*tkt-6.6e-5./tkt.*(1-exp(-tkt./8.0e-4))); % Eq.16 Shortwave |
| 218 |
qcol=qout-dels; |
| 219 |
alq=Al.*qcol+be.*hlb.*cpw./Le; % Eq. 7 Buoy flux water |
| 220 |
|
| 221 |
if alq>0; |
| 222 |
xlamx=6./(1+(bigc.*alq./usr.^4).^.75).^.333; % Eq 13 Saunders |
| 223 |
tkt=xlamx.*visw./(sqrt(rhoa./rhow).*usr); %Eq.11 Sub. thk |
| 224 |
else |
| 225 |
xlamx=6.0; |
| 226 |
tkt=min(.01,xlamx.*visw./(sqrt(rhoa./rhow).*usr)); %Eq.11 Sub. thk |
| 227 |
end; |
| 228 |
|
| 229 |
dter=qcol.*tkt./tcw;% Eq.12 Cool skin |
| 230 |
dqer=wetc.*dter; |
| 231 |
|
| 232 |
end;%bulk iter loop |
| 233 |
|
| 234 |
tau = rhoa.*usr.*usr.*du./ut; %stress |
| 235 |
hsb = -rhoa.*cpa.*usr.*tsr; |
| 236 |
hlb = -rhoa.*Le.*usr.*qsr; |
| 237 |
Rnl = 0.97.*(5.67e-8.*(ts-dter.*jcool+tdk).^4-Rl); |
| 238 |
|
| 239 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.* rain heat flux .*.*.*.*.*.*.*.* |
| 240 |
if 0 |
| 241 |
dwat=2.11e-5.*((t+tdk)./tdk).^1.94;%! water vapour diffusivity |
| 242 |
dtmp=(1.+3.309e-3.*t-1.44e-6.*t.*t).*0.02411./(rhoa.*cpa); %!heat diffusivity |
| 243 |
alfac= 1./(1+(wetc.*Le.*dwat)./(cpa.*dtmp)); %! wet bulb factor |
| 244 |
RF= rain.*alfac.*cpw.*((ts-t-dter.*jcool)+(Qs-Q-dqer.*jcool).*Le./cpa)./3600; |
| 245 |
%.*.*.*.*.*.* add the momentum flux due to rainfall by kelan .*.*.*.*.*.* |
| 246 |
|
| 247 |
tau_r =0.85 .*rain./3600.*du; %%du==wind spd relative to sea surface |
| 248 |
|
| 249 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.*.*.* Webb et al. correection .*.*.*.*.*.*.*.*.*.*.*.* |
| 250 |
wbar=1.61.*hlb./Le./(1+1.61.*Q)./rhoa+hsb./rhoa./cpa./ta;%formulation in hlb already includes webb |
| 251 |
%wbar=1.61.*hlb./Le./rhoa+(1+1.61.*Q).*hsb./rhoa./cpa./ta; |
| 252 |
hl_webb=rhoa.*wbar.*Q.*Le; |
| 253 |
%.*.*.*.*.*.*.*.*.*.*.*.*.*.* compute transfer coeffs relative to ut @meas. ht .*.*. |
| 254 |
Cd=tau./rhoa./ut./max(.1,du); |
| 255 |
Ch=-usr.*tsr./ut./(dt-dter.*jcool); |
| 256 |
Ce=-usr.*qsr./(dq-dqer.*jcool)./ut; |
| 257 |
%.*.*.*.*.*.*.*.*.*.*.*.* 10-m neutral coeff realtive to ut .*.*.*.*.*.*.*.* |
| 258 |
Cdn_10=von.*von./log(zsu./zo)./log(zsu./zo); %von=0.4; == Frank's CD |
| 259 |
Chn_10=von.*von.*fdg./log(zsu./zo)./log(zsu./zot); |
| 260 |
Cen_10=von.*von.*fdg./log(zsu./zo)./log(zsu./zoq); |
| 261 |
end % 0/1 |
| 262 |
Cdn_10=von.*von./log(zsu./zo)./log(zsu./zo); %von=0.4; == Frank's CD |
| 263 |
|
| 264 |
|
| 265 |
%.*.*.*Kelan adds azeta() to adjust met variables to new heights |
| 266 |
%.*.*.*copy from cor3_0af.for |
| 267 |
% zsu % adjusted wind speed height |
| 268 |
% zst % adjusted air temp height |
| 269 |
% zsq % adjusted specific humidity height |
| 270 |
% ts - dter.*jcool === sst of Bradley's; |
| 271 |
|
| 272 |
if 0 |
| 273 |
% kelan add this function call azeta.m in matlab version |
| 274 |
% Frank missed this call in his cor30a.for version |
| 275 |
zetU=azeta(von,grav,ta,Q,usr,tsr,qsr,zsu); |
| 276 |
zetT=azeta(von,grav,ta,Q,usr,tsr,qsr,zst); |
| 277 |
% same as above if zst=zsu |
| 278 |
zetQ=azeta(von,grav,ta,Q,usr,tsr,qsr,zsq); |
| 279 |
% same as above if zsq=zst=zsu |
| 280 |
|
| 281 |
u_zs=usr./von.*(log(zsu./zo)-psiu_30(zetU)); |
| 282 |
t_zs=(ts-dter.*jcool)+tsr./von.*(log(zst./zot)-psit_30(zetT)); |
| 283 |
q_zs=(Qs-dqer)+qsr./von.*(log(zsq./zoq)-psit_30(zetQ)); % !kg./kg |
| 284 |
qa_zs=1000..*q_zs; % !g./kg |
| 285 |
|
| 286 |
e_zs=q_zs.*P./(0.62197+0.378.*q_zs); % !mb |
| 287 |
|
| 288 |
%call humidity(t_zs,p,esat_zs); % !mb |
| 289 |
esat_zs = (1.0007+3.46e-6.*P).*6.1121.*exp(17.502.*t_zs./(240.97+t_zs)); |
| 290 |
rh_zs=e_zs./esat_zs; |
| 291 |
end % 0/1 |
| 292 |
|
| 293 |
% Output: |
| 294 |
y(1,:,:) = -reshape(hsb,[ny nx]); |
| 295 |
y(2,:,:) = -reshape(hlb,[ny nx]); |
| 296 |
y(3,:,:) = -reshape(Rnl,[ny nx]); |
| 297 |
y(4,:,:) = reshape(Rns,[ny nx]); |
| 298 |
y(5,:,:) = reshape(tau,[ny nx]); |
| 299 |
y(6,:,:) = reshape(Cdn_10,[ny nx]); |
| 300 |
y(7,:,:) = reshape(rhoa,[ny nx]); |
| 301 |
|
| 302 |
|
| 303 |
% compute wind stress in two directions: |
| 304 |
ustress = rhoa.*usr.*usr.*uwind_in(:)./ut; |
| 305 |
vstress = rhoa.*usr.*usr.*vwind_in(:)./ut; |
| 306 |
% compute evaporation: |
| 307 |
flamb=2500000; |
| 308 |
evap=hlb/flamb; |
| 309 |
% change sign and convert from kg/m^2/s to m/s via rhoConstFresh |
| 310 |
rhoConstFresh=999.8; |
| 311 |
evap = -evap/rhoConstFresh; |
| 312 |
% reshape: |
| 313 |
hlb=-reshape(hlb,[ny nx]); |
| 314 |
hsb=-reshape(hsb,[ny nx]); |
| 315 |
ustress=reshape(ustress,[ny nx]); |
| 316 |
vstress=reshape(vstress,[ny nx]); |
| 317 |
evap=-reshape(evap,[ny nx]); |
| 318 |
|
| 319 |
%mystruct_out=struct('hl',hlb,'hs',hsb,'evap',evap,'ustress',ustress,'vstress',vstress); |
| 320 |
mystruct_out=struct('hl',real(hlb),'hs',real(hsb),'evap',real(evap),'ustress',real(ustress),'vstress',real(vstress)); |
| 321 |
|
| 322 |
|
| 323 |
function psi=psit_30(zet) |
| 324 |
x=(1-15*zet).^.5; |
| 325 |
psik=2*log((1+x)/2); |
| 326 |
x=(1-34.15*zet).^.3333; |
| 327 |
psic=1.5*log((1+x+x.*x)/3)-sqrt(3)*atan((1+2*x)/sqrt(3))+4*atan(1)/sqrt(3); |
| 328 |
f=zet.*zet./(1+zet.*zet); |
| 329 |
psi=(1-f).*psik+f.*psic; |
| 330 |
|
| 331 |
ii=find(zet>0); |
| 332 |
if ~isempty(ii); |
| 333 |
%psi=-4.7*zet; |
| 334 |
c=min(50,.35*zet); |
| 335 |
psi(ii)=-((1+2/3*zet(ii)).^1.5+.6667*(zet(ii)-14.28)./exp(c(ii))+8.525); |
| 336 |
end; |
| 337 |
|
| 338 |
|
| 339 |
function psi=psiu_30(zet) |
| 340 |
|
| 341 |
x=(1-15*zet).^.25; |
| 342 |
psik=2*log((1+x)/2)+log((1+x.*x)/2)-2*atan(x)+2*atan(1); |
| 343 |
x=(1-10.15*zet).^.3333; |
| 344 |
psic=1.5*log((1+x+x.*x)/3)-sqrt(3)*atan((1+2*x)/sqrt(3))+4*atan(1)/sqrt(3); |
| 345 |
f=zet.*zet./(1+zet.*zet); |
| 346 |
psi=(1-f).*psik+f.*psic; |
| 347 |
ii=find(zet>0); |
| 348 |
if ~isempty(ii); |
| 349 |
|
| 350 |
%psi(ii)=-4.7*zet(ii); |
| 351 |
%c(ii)=min(50,.35*zet(ii)); |
| 352 |
c=min(50,.35*zet); |
| 353 |
psi(ii)=-((1+1.0*zet(ii)).^1.0+.667*(zet(ii)-14.28)./exp(c(ii))+8.525); |
| 354 |
end; |
| 355 |
|
| 356 |
function s=qsee(y) |
| 357 |
x=y(:,1); |
| 358 |
p=y(:,2); |
| 359 |
es=6.112.*exp(17.502.*x./(x+240.97))*.98.*(1.0007+3.46e-6*p); |
| 360 |
s=es*621.97./(p-.378*es); |
| 361 |
|
| 362 |
function y=qsat(y) |
| 363 |
x=y(:,1);%temp |
| 364 |
p=y(:,2);%pressure |
| 365 |
es=6.112.*exp(17.502.*x./(x+241.0)).*(1.0007+3.46e-6*p); |
| 366 |
y=es*622./(p-.378*es); |
| 367 |
|
| 368 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph