1 |
C $Header: /u/gcmpack/MITgcm/pkg/bulk_force/bulkf_formula_lay.F,v 1.1 2006/05/25 17:32:29 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "BULK_FORCE_OPTIONS.h" |
5 |
|
6 |
CBOP |
7 |
C !ROUTINE: BULKF_FORMULA_LAY |
8 |
C !INTERFACE: |
9 |
SUBROUTINE BULKF_FORMULA_LAY( |
10 |
I uw, vw, ws, Ta, Qa, tsfCel, |
11 |
O flwupa, flha, fsha, df0dT, |
12 |
O ust, vst, evp, ssq, dEvdT, |
13 |
I iceornot, i,j,bi,bj,myThid ) |
14 |
|
15 |
C !DESCRIPTION: \bv |
16 |
C *==========================================================* |
17 |
C | SUBROUTINE BULKF_FORMULA_LAY |
18 |
C | o Calculate bulk formula fluxes over open ocean or seaice |
19 |
C | Large and Yeager, 2004, NCAR/TN-460+STR. |
20 |
C *==========================================================* |
21 |
C \ev |
22 |
C |
23 |
C === Turbulent Fluxes : |
24 |
C * use the approach "B": shift coeff to height & stability of the |
25 |
C atmosphere state (instead of "C": shift temp & humid to the height |
26 |
C of wind, then shift the coeff to this height & stability of the atmos). |
27 |
C * similar to EXF (except over sea-ice) ; default parameter values |
28 |
C taken from Large & Yeager. |
29 |
C * assume that Qair & Tair inputs are from the same height (zq=zt) |
30 |
C * formulae in short: |
31 |
C wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v) |
32 |
C Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir |
33 |
C Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap |
34 |
C = -Evap * Lvap |
35 |
C with Ws = wind speed = sqrt(del.u^2 +del.v^2) ; |
36 |
C del.T = Tair - Tsurf ; del.Q = Qair - Qsurf ; |
37 |
C Cd,Ch,Ce = drag coefficient, Stanton number and Dalton number |
38 |
C respectively [no-units], function of height & stability |
39 |
|
40 |
C !USES: |
41 |
IMPLICIT NONE |
42 |
C === Global variables === |
43 |
#include "EEPARAMS.h" |
44 |
#include "SIZE.h" |
45 |
#include "PARAMS.h" |
46 |
#include "BULKF_PARAMS.h" |
47 |
|
48 |
C !INPUT/OUTPUT PARAMETERS: |
49 |
C input: |
50 |
_RL uw ! zonal wind speed (at grid center) [m/s] |
51 |
_RL vw ! meridional wind speed (at grid center) [m/s] |
52 |
_RL ws ! wind speed [m/s] at height zwd |
53 |
_RL Ta ! air temperature [K] at height zth |
54 |
_RL Qa ! specific humidity [kg/kg] at heigth zth |
55 |
_RL tsfCel ! sea-ice or sea surface temperature [oC] |
56 |
INTEGER iceornot ! 0=open water, 1=sea-ice, 2=sea-ice with snow |
57 |
INTEGER i,j, bi,bj !current grid point indices |
58 |
INTEGER myThid ! my Thread Id number |
59 |
C output: |
60 |
_RL flwupa ! upward long wave radiation (>0 upward) [W/m2] |
61 |
_RL flha ! latent heat flux (>0 downward) [W/m2] |
62 |
_RL fsha ! sensible heat flux (>0 downward) [W/m2] |
63 |
_RL df0dT ! derivative of heat flux with respect to Tsf [W/m2/K] |
64 |
_RL ust ! zonal wind stress (at grid center) [N/m2] |
65 |
_RL vst ! meridional wind stress (at grid center)[N/m2] |
66 |
_RL evp ! evaporation rate (over open water) [kg/m2/s] |
67 |
_RL ssq ! surface specific humidity [kg/kg] |
68 |
_RL dEvdT ! derivative of evap. with respect to Tsf [kg/m2/s/K] |
69 |
CEOP |
70 |
|
71 |
#ifdef ALLOW_BULK_FORCE |
72 |
|
73 |
C == Local variables == |
74 |
_RL dflhdT ! derivative of latent heat with respect to T |
75 |
_RL dfshdT ! derivative of sensible heat with respect to T |
76 |
_RL dflwupdT ! derivative of long wave with respect to T |
77 |
|
78 |
_RL Tsf ! surface temperature [K] |
79 |
_RL Ts2 ! surface temperature square [K^2] |
80 |
c _RL ht ! height for air temperature [m] |
81 |
c _RL hq ! height for humidity [m] |
82 |
c _RL hu ! height for wind speed [m] |
83 |
c _RL zref ! reference height [m] |
84 |
_RL wsm ! limited wind speed [m/s] (> umin) |
85 |
_RL usn ! neutral, zref (=10m) wind speed [m/s] |
86 |
_RL usm ! usn but limited [m/s] (> umin) |
87 |
c _RL umin ! minimum wind speed used for drag-coeff [m/s] |
88 |
_RL lath ! latent heat of vaporization or sublimation [J/kg] |
89 |
_RL t0 ! virtual temperature [K] |
90 |
_RL delth ! potential temperature diff [K] |
91 |
_RL delq ! specific humidity difference [kg/kg] |
92 |
_RL ustar ! friction velocity [m/s] |
93 |
_RL tstar ! temperature scale [K] |
94 |
_RL qstar ! humidity scale [kg/kg] |
95 |
_RL rd ! = sqrt(Cd) [-] |
96 |
_RL re ! = Ce / sqrt(Cd) [-] |
97 |
_RL rh ! = Ch / sqrt(Cd) [-] |
98 |
_RL rdn, ren, rhn ! neutral, zref (=10m) values of rd, re, rh |
99 |
_RL stable ! = 1 if stable ; = 0 if unstable |
100 |
_RL huol ! stability parameter at zwd [-] (=z/Monin-Obuklov length) |
101 |
_RL htol ! stability parameter at zth [-] |
102 |
_RL x ! stability function [-] |
103 |
_RL xsq ! = x^2 [-] |
104 |
_RL psimh ! momentum stability function |
105 |
_RL psixh ! latent & sensib. stability function |
106 |
_RL czol ! = zref*Karman_cst*gravity |
107 |
_RL zwln ! = log(zwd/zref) |
108 |
_RL ztln ! = log(zth/zref) |
109 |
c _RL cdalton ! coeff to evaluate Dalton Number |
110 |
c _RL mixratio |
111 |
c _RL ea |
112 |
c _RL psim_fac |
113 |
_RL tau ! surface stress coef = rhoA * Ws * Cd |
114 |
_RL csha ! sensib.heat flx coef = rhoA * Ws * Ch * CpAir |
115 |
_RL clha ! latent heat flx coef = rhoA * Ws * Ce * Lvap |
116 |
c _RL zice |
117 |
c _RL ssq0, ssq1, ssq2 ! constant used in saturated specific humidity |
118 |
c _RL p0 ! reference sea-level atmospheric pressure [mb] |
119 |
_RL qs1w, qs2w ! above freezing saturated specific humidity |
120 |
_RL qs1i, qs2i ! below freezing saturated specific humidity |
121 |
_RL tmpBlk |
122 |
_RL half, one, two |
123 |
INTEGER iter |
124 |
|
125 |
C == external Functions |
126 |
|
127 |
C-- Constant |
128 |
DATA half, one, two |
129 |
& / 0.5 _d 0 , 1. _d 0 , 2. _d 0 / |
130 |
c DATA ssq0, ssq1, ssq2 |
131 |
c & / 3.797915 _d 0 , 7.93252 _d -6 , 2.166847 _d -3 / |
132 |
c DATA p0 / 1013. _d 0 / |
133 |
DATA qs1w, qs2w |
134 |
& / 640.38 _d 3 , 5107.0 _d -0 / |
135 |
DATA qs1i, qs2i |
136 |
& / 11637.80 _d 3 , 5897.8 _d -0 / |
137 |
|
138 |
C-- Set surface parameters : |
139 |
c zice = 0.0005 _d 0 |
140 |
zwln = LOG(zwd/zref) |
141 |
ztln = LOG(zth/zref) |
142 |
czol = zref*xkar*gravity |
143 |
|
144 |
C- Surface Temp. |
145 |
Tsf = tsfCel+Tf0kel |
146 |
Ts2 = Tsf*Tsf |
147 |
C- Wind speed |
148 |
IF (ws.EQ.0. _d 0) THEN |
149 |
ws = SQRT(uw*uw + vw*vw) |
150 |
ENDIF |
151 |
wsm = MAX(ws,umin) |
152 |
|
153 |
C--- Compute turbulent surface fluxes |
154 |
C- Pot. Temp and saturated specific humidity |
155 |
t0 = Ta*(one + humid_fac*Qa) |
156 |
IF ( iceornot.EQ.0 ) THEN |
157 |
lath=Lvap |
158 |
ssq = saltQsFac*qs1w*EXP( -qs2w/Tsf ) / rhoA |
159 |
dEvdT = qs2w |
160 |
ELSE |
161 |
lath = Lvap+Lfresh |
162 |
ssq = qs1i*EXP( -qs2i/Tsf ) / rhoA |
163 |
dEvdT = qs2i |
164 |
ENDIF |
165 |
c ssq = ssq0*EXP( lath*(ssq1-ssq2/Tsf) ) / p0 |
166 |
c dEvdT = lath*ssq2 |
167 |
|
168 |
delth = Ta + gamma_blk*zth - Tsf |
169 |
delq = Qa - ssq |
170 |
|
171 |
C-- initial guess for exchange coefficients: |
172 |
C take U_N = del.U ; stability from del.Theta ; |
173 |
stable = half + SIGN(half, delth) |
174 |
tmpBlk = cdrag_1/wsm + cdrag_2 + cdrag_3*wsm |
175 |
rdn = SQRT(tmpBlk) |
176 |
rhn = stable*cStantonS + (one-stable)*cStantonU |
177 |
ren = cDalton |
178 |
c rdn=xkar/(LOG(zref/zice)) |
179 |
c rhn=rdn |
180 |
c ren=rdn |
181 |
C-- calculate turbulent scales |
182 |
ustar=rdn*wsm |
183 |
tstar=rhn*delth |
184 |
qstar=ren*delq |
185 |
|
186 |
C--- iterate with psi-functions to find transfer coefficients |
187 |
DO iter=1,blk_nIter |
188 |
|
189 |
huol = ( tstar/t0 |
190 |
& +qstar/(Qa + one/humid_fac) |
191 |
& )*czol/(ustar*ustar) |
192 |
huol = SIGN( MIN(abs(huol),10. _d 0), huol) |
193 |
stable = half + SIGN(half, huol) |
194 |
xsq = SQRT( ABS(one - huol*16. _d 0) ) |
195 |
x = SQRT(xsq) |
196 |
psimh = -5. _d 0*huol*stable |
197 |
& + (one-stable)* |
198 |
& ( LOG( (one + two*x + xsq)*(one+xsq)*.125 ) |
199 |
& -two*ATAN(x) + half*pi ) |
200 |
htol = huol*zth/zwd |
201 |
xsq = SQRT( ABS(one - htol*16. _d 0) ) |
202 |
psixh = -5. _d 0*htol*stable |
203 |
& + (one-stable)*( two*LOG(half*(one+xsq)) ) |
204 |
|
205 |
C- Shift wind speed using old coefficient |
206 |
usn = ws/(one + rdn/xkar*(zwln-psimh) ) |
207 |
usm = MAX(usn, umin) |
208 |
|
209 |
C- Update the 10m, neutral stability transfer coefficients |
210 |
tmpBlk = cdrag_1/usm + cdrag_2 + cdrag_3*usm |
211 |
rdn = SQRT(tmpBlk) |
212 |
rhn = stable*cStantonS + (one-stable)*cStantonU |
213 |
ren = cDalton |
214 |
|
215 |
C- Shift all coefficients to the measurement height and stability. |
216 |
rd = rdn/(1. _d 0 + rdn*(zwln-psimh)/xkar) |
217 |
rh = rhn/(1. _d 0 + rhn*(ztln-psixh)/xkar) |
218 |
re = ren/(1. _d 0 + ren*(ztln-psixh)/xkar) |
219 |
|
220 |
C-- Update ustar, tstar, qstar using updated, shifted coefficients. |
221 |
ustar = rd*wsm |
222 |
qstar = re*delq |
223 |
tstar = rh*delth |
224 |
|
225 |
ENDDO |
226 |
|
227 |
C- Coeff: |
228 |
tau = rhoA*rd*ws |
229 |
csha = cpAir*tau*rh |
230 |
clha = lath*tau*re |
231 |
|
232 |
C- Turbulent Fluxes |
233 |
fsha = csha*delth |
234 |
flha = clha*delq |
235 |
evp = -flha/lath |
236 |
ust = tau*rd*uw |
237 |
vst = tau*rd*vw |
238 |
|
239 |
C- surf.Temp derivative of turbulent Fluxes |
240 |
dEvdT = tau*re*ssq*dEvdT/Ts2 |
241 |
dflhdT = -lath*dEvdT |
242 |
dfshdT = -csha |
243 |
|
244 |
C--- Upward long wave radiation |
245 |
IF ( iceornot.EQ.0 ) THEN |
246 |
flwupa = ocean_emissivity*stefan*Ts2*Ts2 |
247 |
dflwupdT= ocean_emissivity*stefan*Ts2*Tsf*4. _d 0 |
248 |
ELSEIF (iceornot.EQ.2) THEN |
249 |
flwupa = snow_emissivity*stefan*Ts2*Ts2 |
250 |
dflwupdT = snow_emissivity*stefan*Ts2*Tsf*4. _d 0 |
251 |
ELSE |
252 |
flwupa = ice_emissivity*stefan*Ts2*Ts2 |
253 |
dflwupdT = ice_emissivity*stefan*Ts2*Tsf*4. _d 0 |
254 |
ENDIF |
255 |
|
256 |
C- Total derivative with respect to surface temperature |
257 |
df0dT = -dflwupdT+dfshdT+dflhdT |
258 |
|
259 |
#endif /*ALLOW_BULK_FORCE*/ |
260 |
|
261 |
RETURN |
262 |
END |