| 1 |
gforget |
1.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