-
Notifications
You must be signed in to change notification settings - Fork 0
/
surface.f90
349 lines (276 loc) · 11.4 KB
/
surface.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
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
subroutine surface()
use vars
use params
use microphysics, only: micro_field, index_water_vapor
implicit none
real qvs(0:nx,1-YES3D:ny),t_s, q_s, u_h0
real taux0, tauy0, xlmo
real diag_ustar, coef, coef1
real fluxq0_coef, ustar0, u10n
real ws0
integer i,j
real(8) buffer(2), buffer1(2)
real dummy1, dummy2
! LES mode:
call t_startf ('surface')
if(.not.SFC_FLX_FXD) then
if(sstxy(1,1).le.-100.) then
print*,'surface: sst is undefined. Quitting...'
call task_abort()
end if
if(OCEAN) then
if(LES) then
call oceflx(pres(1),u0(1)+ug, v0(1)+vg,t0(1)-gamaz(1), q0(1),t0(1),z(1),&
sstxy(1,1)+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
fluxq0_coef, ustar0, u10n)
if(SFC_TAU_FXD) then
u_h0 = max(1.,sqrt((u0(1)+ug)**2+(v0(1)+vg)**2))
taux0 = -(u0(1)+ug)/u_h0*tau0*rhow(1)
tauy0 = -(v0(1)+vg)/u_h0*tau0*rhow(1)
else
tau0=sqrt( taux0**2 + tauy0**2)/rhow(1)
end if
if(doFixedWindSpeedForSurfaceFluxes) then
call oceflx(pres(1),WindSpeedForFluxes,0.,t0(1)-gamaz(1), q0(1),t0(1),z(1),&
sstxy(1,1)+t00, fluxt0, fluxq0, dummy1, dummy2, q_s, &
fluxq0_coef, ustar0, u10n)
end if
fluxbt(:,:) = fluxt0
fluxbq(:,:) = fluxq0
fluxbu(:,:) = taux0/rhow(1)
fluxbv(:,:) = tauy0/rhow(1)
! extra parameters for surface fluxes of water isotopologues
fluxbq_coef(:,:) = fluxq0_coef
qsat_surf(:,:) = q_s
ustar(:,:) = ustar0
u10arr(:,:) = u10n !bloss(2018-02): For aerosol surface flux in M2005_PA
end if ! LES
if(CEM) then
qvs(0:nx,1-YES3D:ny) = micro_field(0:nx,1-YES3D:ny,1,index_water_vapor)
do j=1,ny
do i=1,nx
if(doFixedWindSpeedForSurfaceFluxes) then
! fixed wind speed
call oceflx(pres(1),WindSpeedForFluxes, 0., &
t(i,j,1)-gamaz(1),qv(i,j,1),t(i,j,1),z(1), &
sstxy(i,j)+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
fluxq0_coef, ustar0, u10n)
else
! normal call
call oceflx(pres(1),0.5*(u(i+1,j,1)+u(i,j,1))+ug, &
0.5*(v(i,j+YES3D,1)+v(i,j,1))+vg, &
t(i,j,1)-gamaz(1),qv(i,j,1),t(i,j,1),z(1), &
sstxy(i,j)+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
fluxq0_coef, ustar0, u10n)
end if
fluxbt(i,j) = fluxt0
fluxbq(i,j) = fluxq0
! extra parameters for surface fluxes of water isotopologues
fluxbq_coef(i,j) = fluxq0_coef
qsat_surf(i,j) = q_s
ustar(i,j) = ustar0
u10arr(i,j) = u10n
call oceflx(pres(1),u(i,j,1)+ug, &
0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg, &
0.5*(t(i-1,j,1)+t(i,j,1))-gamaz(1),0.5*(qvs(i-1,j)+qvs(i,j)), &
0.5*(t(i-1,j,1)+t(i,j,1)),z(1), &
0.5*(sstxy(i-1,j)+sstxy(i,j))+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
fluxq0_coef, ustar0, u10n)
if(SFC_TAU_FXD) then
u_h0 = max(1.,sqrt((u(i,j,1)+ug)**2+ &
(0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg)**2))
taux0 = -(u(i,j,1)+ug)/u_h0*tau0*rhow(1)
end if
fluxbu(i,j) = taux0/rhow(1)
call oceflx(pres(1),0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug, &
v(i,j,1)+vg, &
0.5*(t(i,j-YES3D,1)+t(i,j,1))-gamaz(1),0.5*(qvs(i,j-YES3D)+qvs(i,j)), &
0.5*(t(i,j-YES3D,1)+t(i,j,1)),z(1), &
0.5*(sstxy(i,j-YES3D)+sstxy(i,j))+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
fluxq0_coef, ustar0, u10n)
if(SFC_TAU_FXD) then
u_h0 = max(1.,sqrt( &
(0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug)**2+ &
(v(i,j,1)+vg)**2))
tauy0 = -(v(i,j,1)+vg)/u_h0*tau0*rhow(1)
end if
fluxbv(i,j) = tauy0/rhow(1)
end do
end do
end if ! CEM
end if ! OCEAN
if(LAND) then
if(doFixedWindSpeedForSurfaceFluxes) then
write(*,*) 'Fixed Wind Speed not yet implemented for land fluxes'
STOP 'in surface.f90'
end if
if(LES) then
coef = (1000./pres0)**(rgas/cp)
coef1 = (1000./pres(1))**(rgas/cp)
t_s = (sstxy(1,1)+t00)*coef
q_s = soil_wetness*qsatw(sstxy(1,1)+t00,pres(1))
call landflx(pres(1),(t0(1)-gamaz(1))*coef1, t_s, &
q0(1), q_s, u0(1)+ug, v0(1)+vg, z(1), z0, &
fluxt0, fluxq0, taux0, tauy0, xlmo)
if(SFC_TAU_FXD) then
u_h0 = max(1.,sqrt((u0(1)+ug)**2+(v0(1)+vg)**2))
taux0 = -(u0(1)+ug)/u_h0*tau0*rhow(1)
tauy0 = -(v0(1)+vg)/u_h0*tau0*rhow(1)
else
tau0=sqrt( taux0**2 + tauy0**2)/rhow(1)
end if
fluxbt(:,:) = fluxt0
fluxbq(:,:) = fluxq0
fluxbu(:,:) = taux0/rhow(1)
fluxbv(:,:) = tauy0/rhow(1)
end if ! LES
if(CEM) then
coef = (1000./pres0)**(rgas/cp)
coef1 = (1000./pres(1))**(rgas/cp)
qvs(0:nx,1-YES3D:ny) = micro_field(0:nx,1-YES3D:ny,1,index_water_vapor)
do j=1,ny
do i=1,nx
t_s = (sstxy(i,j)+t00)*coef
q_s = soil_wetness*qsatw(sstxy(i,j)+t00,pres(1))
call landflx(pres(1),(t(i,j,1)-gamaz(1))*coef1, t_s, &
qv(i,j,1), q_s, 0.5*(u(i+1,j,1)+u(i,j,1))+ug, &
0.5*(v(i,j+YES3D,1)+v(i,j,1))+vg, z(1), z0, &
fluxt0, fluxq0, taux0, tauy0, xlmo)
fluxbt(i,j) = fluxt0
fluxbq(i,j) = fluxq0
t_s = (0.5*(sstxy(i-1,j)+sstxy(i,j))+t00)*coef
q_s = soil_wetness*qsatw(0.5*(sstxy(i-1,j)+sstxy(i,j))+t00,pres(1))
call landflx(pres(1),(0.5*(t(i-1,j,1)+t(i,j,1))-gamaz(1))*coef1, t_s, &
0.5*(qvs(i-1,j)+qvs(i,j)), q_s, u(i,j,1)+ug, &
0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg, &
z(1), z0, fluxt0, fluxq0, taux0, tauy0, xlmo)
if(SFC_TAU_FXD) then
u_h0 = max(1.,sqrt((u(i,j,1)+ug)**2+ &
(0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg)**2))
taux0 = -(u(i,j,1)+ug)/u_h0*tau0*rhow(1)
end if
fluxbu(i,j) = taux0/rhow(1)
t_s = (0.5*(sstxy(i,j-YES3D)+sstxy(i,j))+t00)*coef
q_s = soil_wetness*qsatw(0.5*(sstxy(i,j-YES3D)+sstxy(i,j))+t00,pres(1))
call landflx(pres(1),(0.5*(t(i,j-YES3D,1)+t(i,j,1))-gamaz(1))*coef1, t_s, &
0.5*(qvs(i,j-YES3D)+qvs(i,j)), q_s, &
0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug, &
v(i,j,1)+vg, &
z(1), z0, fluxt0, fluxq0, taux0, tauy0, xlmo)
if(SFC_TAU_FXD) then
u_h0 = max(1.,sqrt( &
(0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug)**2+ &
(v(i,j,1)+vg)**2))
tauy0 = -(v(i,j,1)+vg)/u_h0*tau0*rhow(1)
end if
fluxbv(i,j) = tauy0/rhow(1)
end do
end do
end if ! CEM
end if ! LAND
end if! .not.SFC_FLX_FXD
if(SFC_FLX_FXD) then
u_h0 = max(1.,sqrt((u0(1)+ug)**2+(v0(1)+vg)**2))
if(.not.SFC_TAU_FXD) then
if(OCEAN) z0 = 0.0001 ! for LAND z0 should be set in namelist (default z0=0.035)
tau0 = diag_ustar(z(1), &
bet(1)*(fluxt0+epsv*(t0(1)-gamaz(1))*fluxq0),u_h0,z0)**2
end if ! .not.SFC_TAU_FXD
if(LES) then
taux0 = -(u0(1)+ug)/u_h0*tau0
tauy0 = -(v0(1)+vg)/u_h0*tau0
fluxbu(:,:) = taux0
fluxbv(:,:) = tauy0
else
fluxbu(:,:) = -(u(1:nx,1:ny,1)+ug)/u_h0*tau0
fluxbv(:,:) = -(v(1:nx,1:ny,1)+vg)/u_h0*tau0
end if
fluxbt(:,:) = fluxt0
fluxbq(:,:) = fluxq0
end if ! SFC_FLX_FXD
!
! Homogenize the surface scalar fluxes if needed for sensitivity studies
!
if(dosfchomo) then
fluxt0 = 0.
fluxq0 = 0.
do j=1,ny
do i=1,nx
fluxt0 = fluxt0 + fluxbt(i,j)
fluxq0 = fluxq0 + fluxbq(i,j)
end do
end do
fluxt0 = fluxt0 / float(nx*ny)
fluxq0 = fluxq0 / float(nx*ny)
if(dompi) then
buffer(1) = fluxt0
buffer(2) = fluxq0
call task_sum_real8(buffer,buffer1,2)
fluxt0 = buffer1(1) /float(nsubdomains)
fluxq0 = buffer1(2) /float(nsubdomains)
end if ! dompi
fluxbt(:,:) = fluxt0
fluxbq(:,:) = fluxq0
end if
shf_xy(:,:) = shf_xy(:,:) + fluxbt(:,:) * dtfactor
lhf_xy(:,:) = lhf_xy(:,:) + fluxbq(:,:) * dtfactor
call t_stopf ('surface')
end
! ----------------------------------------------------------------------
!
! DISCLAIMER : this code appears to be correct but has not been
! very thouroughly tested. If you do notice any
! anomalous behaviour then please contact Andy and/or
! Bjorn
!
! Function diag_ustar: returns value of ustar using the below
! similarity functions and a specified buoyancy flux (bflx) given in
! kinematic units
!
! phi_m (zeta > 0) = (1 + am * zeta)
! phi_m (zeta < 0) = (1 - bm * zeta)^(-1/4)
!
! where zeta = z/lmo and lmo = (theta_rev/g*vonk) * (ustar^2/tstar)
!
! Ref: Businger, 1973, Turbulent Transfer in the Atmospheric Surface
! Layer, in Workshop on Micormeteorology, pages 67-100.
!
! Code writen March, 1999 by Bjorn Stevens
!
! Code corrected 8th June 1999 (obukhov length was wrong way up,
! so now used as reciprocal of obukhov length)
real function diag_ustar(z,bflx,wnd,z0)
implicit none
real, parameter :: vonk = 0.4 ! von Karmans constant
real, parameter :: g = 9.81 ! gravitational acceleration
real, parameter :: am = 4.8 ! " " "
real, parameter :: bm = 19.3 ! " " "
real, parameter :: eps = 1.e-10 ! non-zero, small number
real, intent (in) :: z ! height where u locates
real, intent (in) :: bflx ! surface buoyancy flux (m^2/s^3)
real, intent (in) :: wnd ! wind speed at z
real, intent (in) :: z0 ! momentum roughness height
integer :: iterate
real :: lnz, klnz, c1, x, psi1, zeta, rlmo, ustar
lnz = log(z/z0)
klnz = vonk/lnz
c1 = 3.14159/2. - 3.*log(2.)
ustar = wnd*klnz
if (bflx /= 0.0) then
do iterate=1,4
rlmo = -bflx * vonk/(ustar**3 + eps) !reciprocal of
!obukhov length
zeta = z*rlmo
if (zeta > 0.) then
ustar = vonk*wnd /(lnz + am*zeta)
else
x = sqrt( sqrt( 1.0 - bm*zeta ) )
psi1 = 2.*log(1.0+x) + log(1.0+x*x) - 2.*atan(x) + c1
ustar = wnd*vonk/(lnz - psi1)
end if
end do
end if
diag_ustar = ustar
return
end function diag_ustar
! ----------------------------------------------------------------------