| 1 |
function [mystruct_out]=exf_bulk_largeyeager04(mystruct_in,varargin); |
| 2 |
|
| 3 |
%varargin [hu ht hq] |
| 4 |
% [hu ht hq],dragCoeff |
| 5 |
|
| 6 |
tmp1=fieldnames(mystruct_in); |
| 7 |
if length(tmp1)>5; doNetFluxes=1; end; |
| 8 |
|
| 9 |
% formulae in short: |
| 10 |
% wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v) |
| 11 |
% Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir |
| 12 |
% Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap |
| 13 |
% = -Evap * Lvap |
| 14 |
% with Ws = wind speed = sqrt(del.u^2 +del.v^2) ; |
| 15 |
% del.T = Tair - Tsurf ; del.Q = Qair - Qsurf ; |
| 16 |
% Cd,Ch,Ce = drag coefficient, Stanton number and Dalton number |
| 17 |
% respectively [no-units], function of height & stability |
| 18 |
|
| 19 |
|
| 20 |
%parameters from exf_readparms.F |
| 21 |
if nargin==1 |
| 22 |
hu=10; |
| 23 |
ht=10; %differ from model version because all CORE data are given at 10m |
| 24 |
hq=10; |
| 25 |
else |
| 26 |
hu=varargin{1}(1); |
| 27 |
ht=varargin{1}(2); |
| 28 |
hq=varargin{1}(3); |
| 29 |
end |
| 30 |
zref=10; |
| 31 |
umin=0.5; |
| 32 |
karman=0.4; |
| 33 |
gravity_mks=9.81; |
| 34 |
cen2kel=273.150; |
| 35 |
humid_fac=0.606; |
| 36 |
cvapor_fac=640380; |
| 37 |
cvapor_fac_ice=11637800; |
| 38 |
cvapor_exp=5107.4; |
| 39 |
cvapor_exp_ice=5897.8; |
| 40 |
saltsat=0.980; |
| 41 |
atmrho=1.2; |
| 42 |
gamma_blk=0.01; |
| 43 |
cdrag_1=0.0027000; |
| 44 |
cdrag_2=0.0001420; |
| 45 |
cdrag_3=0.0000764; |
| 46 |
cstanton_1 = 0.0327; |
| 47 |
cstanton_2 = 0.0180; |
| 48 |
cdalton = 0.0346; |
| 49 |
niter_bulk=2; |
| 50 |
psim_fac=5; |
| 51 |
atmcp=1005; |
| 52 |
flamb=2500000; |
| 53 |
rhoConstFresh=999.8; |
| 54 |
stefanBoltzmann = 5.670e-8; |
| 55 |
ocean_emissivity=5.50e-8 / 5.670e-8; |
| 56 |
albedo=0.1; |
| 57 |
|
| 58 |
%-- Set surface parameters : |
| 59 |
zwln = log(hu/zref); |
| 60 |
ztln = log(ht/zref); |
| 61 |
czol = hu*karman*gravity_mks; |
| 62 |
|
| 63 |
%inputs: |
| 64 |
atemp_in=mystruct_in.atemp; |
| 65 |
aqh_in=mystruct_in.aqh; |
| 66 |
uwind_in=mystruct_in.uwind; |
| 67 |
vwind_in=mystruct_in.vwind; |
| 68 |
sst_in=mystruct_in.sst; |
| 69 |
|
| 70 |
|
| 71 |
%compute wind speed: |
| 72 |
wspeed=uwind_in.*uwind_in+vwind_in.*vwind_in; |
| 73 |
wspeed=sqrt(wspeed); wspeed=max(wspeed,umin); |
| 74 |
|
| 75 |
%- Surface Temp. |
| 76 |
Tsf = sst_in + cen2kel; |
| 77 |
|
| 78 |
%--- Compute turbulent surface fluxes |
| 79 |
%- Pot. Temp and saturated specific humidity |
| 80 |
|
| 81 |
t0 = atemp_in.*(1 + humid_fac*aqh_in); |
| 82 |
tmpbulk = cvapor_fac*exp(-cvapor_exp./Tsf); |
| 83 |
ssq = saltsat*tmpbulk/atmrho; |
| 84 |
deltap = atemp_in + gamma_blk*ht - Tsf; |
| 85 |
delq = aqh_in - ssq; |
| 86 |
|
| 87 |
%-- initial guess for exchange coefficients: |
| 88 |
% take U_N = del.U ; stability from del.Theta ; |
| 89 |
stable = 0.5 + 0.5*sign(deltap); |
| 90 |
if nargin<3 |
| 91 |
tmpbulk = cdrag_1./wspeed + cdrag_2 + cdrag_3.*wspeed; |
| 92 |
else |
| 93 |
%use prescribed drag coefficient |
| 94 |
tmpbulk=varargin{2}; |
| 95 |
end |
| 96 |
rdn = sqrt(tmpbulk); |
| 97 |
rhn = (1-stable)*cstanton_1 + stable*cstanton_2; |
| 98 |
ren = cdalton; |
| 99 |
%-- calculate turbulent scales |
| 100 |
ustar=rdn.*wspeed; |
| 101 |
tstar=rhn.*deltap; |
| 102 |
qstar=ren*delq; |
| 103 |
|
| 104 |
%--- iterate with psi-functions to find transfer coefficients |
| 105 |
for iter=1:niter_bulk |
| 106 |
|
| 107 |
huol = ( tstar./t0 + qstar./(1/humid_fac+aqh_in) )*czol./(ustar.*ustar); |
| 108 |
|
| 109 |
%cph The following is different in Large&Pond1981 code: |
| 110 |
%cph huol = max(huol,zolmin) with zolmin = -100 |
| 111 |
tmpbulk = min(abs(huol),10); |
| 112 |
huol = tmpbulk.*sign(huol); |
| 113 |
htol = huol*ht/hu; |
| 114 |
hqol = huol*hq/hu; |
| 115 |
stable = 0.5 + 0.5*sign(huol); |
| 116 |
|
| 117 |
% Evaluate all stability functions assuming hq = ht. |
| 118 |
% The following is different in Large&Pond1981 code: |
| 119 |
% xsq = max(sqrt(abs(1 - 16.*huol)),1) |
| 120 |
xsq = sqrt( abs(1 - huol*16) ); |
| 121 |
x = sqrt(xsq); |
| 122 |
psimh = -psim_fac*huol.*stable + (1-stable).* ... |
| 123 |
( log( (1 + 2*x + xsq).*(1+xsq)*.125 ) - 2*atan(x) + 0.5*pi ); |
| 124 |
|
| 125 |
|
| 126 |
% The following is different in Large&Pond1981 code: |
| 127 |
% xsq = max(sqrt(abs(1 - 16.*htol)),1) |
| 128 |
xsq = sqrt( abs(1 - htol*16) ); |
| 129 |
psixh = -psim_fac*htol.*stable + (1-stable).*( 2*log(0.5*(1+xsq)) ); |
| 130 |
|
| 131 |
%- Shift wind speed using old coefficient |
| 132 |
usn = wspeed./(1 + rdn/karman.*(zwln-psimh) ); |
| 133 |
usm = max(usn, umin); |
| 134 |
|
| 135 |
%- Update the 10m, neutral stability transfer coefficients |
| 136 |
if nargin<3 |
| 137 |
tmpbulk = cdrag_1./usm + cdrag_2 + cdrag_3.*usm; |
| 138 |
else |
| 139 |
%use prescribed drag coefficient |
| 140 |
tmpbulk=varargin{2}; |
| 141 |
end |
| 142 |
rdn = sqrt(tmpbulk); |
| 143 |
rhn = (1-stable)*cstanton_1 + stable*cstanton_2; |
| 144 |
ren = cdalton; |
| 145 |
|
| 146 |
%- Shift all coefficients to the measurement height and stability. |
| 147 |
rd = rdn./(1 + rdn.*(zwln-psimh)/karman); |
| 148 |
rh = rhn./(1 + rhn.*(ztln-psixh)/karman); |
| 149 |
re = ren./(1 + ren.*(ztln-psixh)/karman); |
| 150 |
|
| 151 |
%-- Update ustar, tstar, qstar using updated, shifted coefficients. |
| 152 |
ustar = rd.*wspeed; |
| 153 |
qstar = re.*delq; |
| 154 |
tstar = rh.*deltap; |
| 155 |
|
| 156 |
% end of iteration loop |
| 157 |
end%for iter=1:niter_bulk |
| 158 |
|
| 159 |
%- Coeff: |
| 160 |
tau = atmrho*rd.*wspeed; |
| 161 |
|
| 162 |
%- Turbulent Fluxes |
| 163 |
hs = atmcp*tau.*tstar; |
| 164 |
hl = flamb*tau.*qstar; |
| 165 |
% change sign and convert from kg/m^2/s to m/s via rhoConstFresh |
| 166 |
evap = -tau.*qstar/rhoConstFresh; |
| 167 |
|
| 168 |
% ustress(i,j,bi,bj) = tau*rd*ws*cw(i,j,bi,bj) |
| 169 |
% vstress(i,j,bi,bj) = tau*rd*ws*sw(i,j,bi,bj) |
| 170 |
%- jmc: below is how it should be written ; different from above when |
| 171 |
% both wind-speed & 2 compon. of the wind are specified, and in |
| 172 |
%- this case, formula below is better. |
| 173 |
ustress = tau.*rd.*uwind_in; |
| 174 |
vstress = tau.*rd.*vwind_in; |
| 175 |
|
| 176 |
mystruct_out=struct('hl',hl,'hs',hs,'evap',evap,'ustress',ustress,'vstress',vstress); |
| 177 |
|
| 178 |
|
| 179 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph