forked from KuangLab-Harvard/SAM_SRCv6.11
-
Notifications
You must be signed in to change notification settings - Fork 0
/
oceflx.f90
239 lines (193 loc) · 8.36 KB
/
oceflx.f90
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
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
subroutine oceflx(rbot, ubot, vbot, tbot, qbot, thbot, zbot, ts, qs, &
shf, lhf, taux, tauy)
use params, only: salt_factor
implicit none
!
! Compute ocean to atmosphere surface fluxes of sensible, latent heat
! and stress components:
!
! Assume:
! 1) Neutral 10m drag coeff:
! cdn = .0027/U10N + .000142 + .0000764 U10N
! 2) Neutral 10m stanton number:
! ctn = .0327 sqrt(cdn), unstable
! ctn = .0180 sqrt(cdn), stable
! 3) Neutral 10m dalton number:
! cen = .0346 sqrt(cdn)
!
! Note:
! 1) here, tstar = <WT>/U*, and qstar = <WQ>/U*.
! 2) wind speeds should all be above a minimum speed (umin)
!
!--------------------------Code History---------------------------------
!
! Original: Bill Large/M.Vertenstein, Sep. 1995 for CCM3.5
! Standardized: L. Buja, Feb 1996
! Reviewed: B. Briegleb, March 1996
! Adopted for LES by Marat Khairoutdinov, July 1998
!
!------------------------------Arguments--------------------------------
!
! Input arguments
!
!
real rbot ! Density at bottom model level
real ubot ! Bottom level u wind
real vbot ! Bottom level v wind
real tbot ! Bottom level temperature
real qbot ! Bottom level specific humidity
real thbot ! Bottom level potential temperature ---> SEE NOTE BELOW.
real zbot ! Bottom level height above surface
real ts ! Surface temperature
real qs ! Surface saturation specific humidity
!
! Output arguments
!
real shf ! Initial sensible heat flux (Km/s)
real lhf ! Initial latent heat flux (m/s)
real taux ! X surface stress (N/m2)
real tauy ! Y surface stress (N/m2)
!
!---------------------------Local variables-----------------------------
!
real ustar ! ustar
real tstar ! tstar
real qstar ! qstar
real u10n ! neutral 10 m wind speed over ocean
real vmag ! Surface wind magnitude
real thvbot ! Bottom lev virtual potential temp
real delt ! potential T difference (K)
real delq ! specific humidity difference (kg/kg)
real rdn ! sqrt of neutral exchange coeff (momentum)
real rhn ! sqrt of neutral exchange coeff (heat)
real ren ! sqrt of neutral exchange coeff (tracers)
real rd ! sqrt of exchange coeff (momentum)
real rh ! sqrt of exchange coeff (heat)
real re ! sqrt of exchange coeff (tracers)
real hol ! Ref hgt (10m) / monin-obukhov length
real xsq ! Temporary variable
real xqq ! Temporary variable
real alz ! ln(zbot/z10)
real cp ! Specific heat of moist air
real tau ! Reference height stress
real psimh ! Stability funct at ref lev (momentum)
real psixh ! Stability funct at ref lev (heat & tracers)
real stable ! Stability factor
real bn ! exchange coef funct for interpolation
real bh ! exchange coef funct for interpolation
real fac ! interpolation factor
real ln0 ! log factor for interpolation
real ln3 ! log factor for interpolation
real ztref ! reference height for air temperature
real gravit ! =9.81 m/s2
real ltheat ! Latent heat for given srf conditions
real xkar ! Von Karman constant
real zref ! 10m reference height
real umin ! Minimum wind speed at bottom level
parameter (xkar = 0.4 , &
zref = 10.0 , &
umin = 1.0, &
ztref=2.0 )
real qsatw, qsati
external qsatw
!--------------------------Statement functions--------------------------
real psimhu ! Unstable part of psimh
real psixhu ! Unstable part of psixh
real cdn ! Neutral drag coeff at bottom model level
real xd ! Dummy argument
real Umps ! Wind velocity (m/sec)
cdn(Umps) = 0.0027 / Umps + .000142 + .0000764 * Umps
psimhu(xd) = log((1.+xd*(2.+xd))*(1.+xd*xd)/8.) - 2.*atan(xd) + 1.571
psixhu(xd) = 2. * log((1. + xd*xd)/2.)
!---------------------------------------------------------------
! Set up necessary variables
!---------------------------------------------------------------
if(zbot.le.zref) then
print*,'oceflx: zbot should be > 10m'
call task_abort
end if
vmag = max(umin, sqrt(ubot **2+vbot **2))
thvbot = thbot * (1.0 + 0.606*qbot )
delt = thbot - ts
delq = qbot - qs
alz = log(zbot /zref)
cp = 1005.
gravit = 9.81
ltheat = 2.51e6
!---------------------------------------------------------------
! First iteration to converge on Z/L and hence the fluxes
!---------------------------------------------------------------
! Initial guess for roots of neutral exchange coefficients,
! assume z/L=0. and u10n is approximated by vmag.
! Stable if (thbot > ts ).
stable = 0.5 + sign(0.5 , delt)
rdn = sqrt(cdn(vmag))
rhn = (1.-stable) * 0.0327 + stable * 0.018
ren = 0.0346
! Initial guess of ustar,tstar and qstar
ustar = rdn*vmag
tstar = rhn*delt
qstar = ren*delq
! Compute stability and evaluate all stability functions
! Stable if (thbot > ts or hol > 0 )
hol = xkar*gravit*zbot*(tstar/thvbot+qstar/(1./0.606+qbot))/ustar**2
hol = sign( min(abs(hol),10.), hol )
stable = 0.5 + sign(0.5 , hol)
xsq = max(sqrt(abs(1. - 16.*hol)) , 1.)
xqq = sqrt(xsq)
psimh = -5. * hol * stable + (1.-stable)*psimhu(xqq)
psixh = -5. * hol * stable + (1.-stable)*psixhu(xqq)
! Shift 10m neutral wind speed using old rdn coefficient
rd = rdn / (1.+rdn/xkar*(alz-psimh))
u10n = vmag * rd / rdn
! Update the neutral transfer coefficients at 10m and neutral stability
rdn = sqrt(cdn(u10n))
ren = 0.0346
rhn = (1.-stable) * 0.0327 + stable * 0.018
! Shift all coeffs to measurement height and stability
rd = rdn / (1.+rdn/xkar*(alz-psimh))
rh = rhn / (1.+rhn/xkar*(alz-psixh))
re = ren / (1.+ren/xkar*(alz-psixh))
! Update ustar, tstar, qstar using updated, shifted coeffs
ustar = rd * vmag
tstar = rh * delt
qstar = re * delq
!---------------------------------------------------------------
! Second iteration to converge on Z/L and hence the fluxes
!---------------------------------------------------------------
! Recompute stability & evaluate all stability functions
! Stable if (thbot > ts or hol > 0 )
hol = xkar*gravit*zbot*(tstar/thvbot+qstar/(1./0.606+qbot))/ustar**2
hol = sign( min(abs(hol),10.), hol )
stable = 0.5 + sign(0.5 , hol)
xsq = max(sqrt(abs(1. - 16.*hol)) , 1.)
xqq = sqrt(xsq)
psimh = -5. * hol * stable + (1.-stable)*psimhu(xqq)
psixh = -5. * hol * stable + (1.-stable)*psixhu(xqq)
! Shift 10m neutral wind speed using old rdn coefficient
rd = rdn / (1.+rdn/xkar*(alz-psimh))
u10n = vmag * rd / rdn
! Update the neutral transfer coefficients at 10m and neutral stability
rdn = sqrt(cdn(u10n))
ren = 0.0346
rhn = (1.-stable) * 0.0327 + stable * 0.018
! Shift all coeffs to measurement height and stability
rd = rdn / (1.+rdn/xkar*(alz-psimh))
rh = rhn / (1.+rhn/xkar*(alz-psixh))
re = ren / (1.+ren/xkar*(alz-psixh))
!---------------------------------------------------------------
! Compute the fluxes
!---------------------------------------------------------------
! Update ustar, tstar, qstar using updated, shifted coeffs
ustar = rd * vmag
tstar = rh * delt
qstar = re * delq
! Compute surface stress components
tau = rbot * ustar * ustar
taux = -tau * ubot / vmag
tauy = -tau * vbot / vmag
! Compute heat flux components at current surface temperature
! (Define positive latent and sensible heat as upwards into the atm)
shf = -tau * tstar / ustar / rbot
lhf = -tau * qstar / ustar / rbot
end subroutine oceflx