ROMS
Loading...
Searching...
No Matches
tl_bulk_flux_mod Module Reference

Functions/Subroutines

subroutine, public tl_bulk_flux (ng, tile)
 
subroutine tl_bulk_flux_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, rmask, umask, vmask, rmask_wet, umask_wet, vmask_wet, alpha, tl_alpha, beta, tl_beta, rho, tl_rho, t, tl_t, hair, pair, tair, uwind, vwind, awave, pwave, rain, lhflx, tl_lhflx, lrflx, tl_lrflx, shflx, tl_shflx, srflx, stflux, tl_stflux, evap, tl_evap, sustr, tl_sustr, svstr, tl_svstr)
 
real(r8) function, public tl_bulk_psiu (tl_zol, zol, pi)
 
real(r8) function, public tl_bulk_psit (tl_zol, zol, pi)
 

Function/Subroutine Documentation

◆ tl_bulk_flux()

subroutine, public tl_bulk_flux_mod::tl_bulk_flux ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 65 of file tl_bulk_flux.F.

66!***********************************************************************
67!
68 USE mod_param
69 USE mod_forces
70 USE mod_grid
71 USE mod_mixing
72 USE mod_ocean
73 USE mod_stepping
74!
75! Imported variable declarations.
76!
77 integer, intent(in) :: ng, tile
78!
79! Local variable declarations.
80!
81 character (len=*), parameter :: MyFile = &
82 & __FILE__
83!
84# include "tile.h"
85!
86# ifdef PROFILE
87 CALL wclock_on (ng, itlm, 17, __line__, myfile)
88# endif
89 CALL tl_bulk_flux_tile (ng, tile, &
90 & lbi, ubi, lbj, ubj, &
91 & imins, imaxs, jmins, jmaxs, &
92 & nrhs(ng), &
93# ifdef MASKING
94 & grid(ng) % rmask, &
95 & grid(ng) % umask, &
96 & grid(ng) % vmask, &
97# endif
98# ifdef WET_DRY
99 & grid(ng) % rmask_wet, &
100 & grid(ng) % umask_wet, &
101 & grid(ng) % vmask_wet, &
102# endif
103 & mixing(ng) % alpha, &
104 & mixing(ng) % tl_alpha, &
105 & mixing(ng) % beta, &
106 & mixing(ng) % tl_beta, &
107 & ocean(ng) % rho, &
108 & ocean(ng) % tl_rho, &
109 & ocean(ng) % t, &
110 & ocean(ng) % tl_t, &
111 & forces(ng) % Hair, &
112 & forces(ng) % Pair, &
113 & forces(ng) % Tair, &
114 & forces(ng) % Uwind, &
115 & forces(ng) % Vwind, &
116# ifdef CLOUDS
117 & forces(ng) % cloud, &
118# endif
119# ifdef BBL_MODEL_NOT_YET
120 & forces(ng) % Awave, &
121 & forces(ng) % Pwave, &
122# endif
123 & forces(ng) % rain, &
124 & forces(ng) % lhflx, &
125 & forces(ng) % tl_lhflx, &
126 & forces(ng) % lrflx, &
127 & forces(ng) % tl_lrflx, &
128 & forces(ng) % shflx, &
129 & forces(ng) % tl_shflx, &
130 & forces(ng) % srflx, &
131 & forces(ng) % stflux, &
132 & forces(ng) % tl_stflux, &
133# ifdef EMINUSP
134 & forces(ng) % evap, &
135 & forces(ng) % tl_evap, &
136# endif
137 & forces(ng) % sustr, &
138 & forces(ng) % tl_sustr, &
139 & forces(ng) % svstr, &
140 & forces(ng) % tl_svstr)
141# ifdef PROFILE
142 CALL wclock_off (ng, itlm, 17, __line__, myfile)
143# endif
144!
145 RETURN
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
integer, parameter itlm
Definition mod_param.F:663
integer, dimension(:), allocatable nrhs
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3

References mod_forces::forces, mod_grid::grid, mod_param::itlm, mod_mixing::mixing, mod_stepping::nrhs, mod_ocean::ocean, tl_bulk_flux_tile(), wclock_off(), and wclock_on().

Here is the call graph for this function:

◆ tl_bulk_flux_tile()

subroutine tl_bulk_flux_mod::tl_bulk_flux_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:,lbj:), intent(in) rmask,
real(r8), dimension(lbi:,lbj:), intent(in) umask,
real(r8), dimension(lbi:,lbj:), intent(in) vmask,
real(r8), dimension(lbi:,lbj:), intent(in) rmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) umask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) vmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) alpha,
real(r8), dimension(lbi:,lbj:), intent(in) tl_alpha,
real(r8), dimension(lbi:,lbj:), intent(in) beta,
real(r8), dimension(lbi:,lbj:), intent(in) tl_beta,
real(r8), dimension(lbi:,lbj:,:), intent(in) rho,
real(r8), dimension(lbi:,lbj:,:), intent(in) tl_rho,
real(r8), dimension(lbi:,lbj:,:,:,:), intent(in) t,
real(r8), dimension(lbi:,lbj:,:,:,:), intent(in) tl_t,
real(r8), dimension(lbi:,lbj:), intent(in) hair,
real(r8), dimension(lbi:,lbj:), intent(in) pair,
real(r8), dimension(lbi:,lbj:), intent(in) tair,
real(r8), dimension(lbi:,lbj:), intent(in) uwind,
real(r8), dimension(lbi:,lbj:), intent(in) vwind,
real(r8), dimension(lbi:,lbj:), intent(in) awave,
real(r8), dimension(lbi:,lbj:), intent(in) pwave,
real(r8), dimension(lbi:,lbj:), intent(in) rain,
real(r8), dimension(lbi:,lbj:), intent(inout) lhflx,
real(r8), dimension(lbi:,lbj:), intent(inout) tl_lhflx,
real(r8), dimension(lbi:,lbj:), intent(inout) lrflx,
real(r8), dimension(lbi:,lbj:), intent(inout) tl_lrflx,
real(r8), dimension(lbi:,lbj:), intent(inout) shflx,
real(r8), dimension(lbi:,lbj:), intent(inout) tl_shflx,
real(r8), dimension(lbi:,lbj:), intent(inout) srflx,
real(r8), dimension(lbi:,lbj:,:), intent(in) stflux,
real(r8), dimension(lbi:,lbj:,:), intent(inout) tl_stflux,
real(r8), dimension(lbi:,lbj:), intent(out) evap,
real(r8), dimension(lbi:,lbj:), intent(out) tl_evap,
real(r8), dimension(lbi:,lbj:), intent(out) sustr,
real(r8), dimension(lbi:,lbj:), intent(out) tl_sustr,
real(r8), dimension(lbi:,lbj:), intent(out) svstr,
real(r8), dimension(lbi:,lbj:), intent(out) tl_svstr )
private

Definition at line 149 of file tl_bulk_flux.F.

181!***********************************************************************
182!
183 USE mod_param
184 USE mod_scalars
185!
188# ifdef DISTRIBUTE
190# endif
191!
192! Imported variable declarations.
193!
194 integer, intent(in) :: ng, tile
195 integer, intent(in) :: LBi, UBi, LBj, UBj
196 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
197 integer, intent(in) :: nrhs
198!
199# ifdef ASSUMED_SHAPE
200# ifdef MASKING
201 real(r8), intent(in) :: rmask(LBi:,LBj:)
202 real(r8), intent(in) :: umask(LBi:,LBj:)
203 real(r8), intent(in) :: vmask(LBi:,LBj:)
204# endif
205# ifdef WET_DRY
206 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
207 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
208 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
209# endif
210 real(r8), intent(in) :: alpha(LBi:,LBj:)
211 real(r8), intent(in) :: beta(LBi:,LBj:)
212 real(r8), intent(in) :: rho(LBi:,LBj:,:)
213 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
214 real(r8), intent(in) :: Hair(LBi:,LBj:)
215 real(r8), intent(in) :: Pair(LBi:,LBj:)
216 real(r8), intent(in) :: Tair(LBi:,LBj:)
217 real(r8), intent(in) :: Uwind(LBi:,LBj:)
218 real(r8), intent(in) :: Vwind(LBi:,LBj:)
219 real(r8), intent(in) :: tl_alpha(LBi:,LBj:)
220 real(r8), intent(in) :: tl_beta(LBi:,LBj:)
221 real(r8), intent(in) :: tl_rho(LBi:,LBj:,:)
222 real(r8), intent(in) :: tl_t(LBi:,LBj:,:,:,:)
223# ifdef CLOUDS
224 real(r8), intent(in) :: cloud(LBi:,LBj:)
225# endif
226# ifdef BBL_MODEL_NOT_YET
227 real(r8), intent(in) :: Awave(LBi:,LBj:)
228 real(r8), intent(in) :: Pwave(LBi:,LBj:)
229# endif
230 real(r8), intent(in) :: rain(LBi:,LBj:)
231 real(r8), intent(in) :: stflux(LBi:,LBj:,:)
232
233 real(r8), intent(inout) :: lhflx(LBi:,LBj:)
234 real(r8), intent(inout) :: lrflx(LBi:,LBj:)
235 real(r8), intent(inout) :: shflx(LBi:,LBj:)
236 real(r8), intent(inout) :: srflx(LBi:,LBj:)
237 real(r8), intent(inout) :: tl_lhflx(LBi:,LBj:)
238 real(r8), intent(inout) :: tl_lrflx(LBi:,LBj:)
239 real(r8), intent(inout) :: tl_shflx(LBi:,LBj:)
240 real(r8), intent(inout) :: tl_stflux(LBi:,LBj:,:)
241
242# ifdef EMINUSP
243 real(r8), intent(out) :: evap(LBi:,LBj:)
244 real(r8), intent(out) :: tl_evap(LBi:,LBj:)
245# endif
246 real(r8), intent(out) :: sustr(LBi:,LBj:)
247 real(r8), intent(out) :: svstr(LBi:,LBj:)
248 real(r8), intent(out) :: tl_sustr(LBi:,LBj:)
249 real(r8), intent(out) :: tl_svstr(LBi:,LBj:)
250# else
251# ifdef MASKING
252 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
253 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
254 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
255# endif
256# ifdef WET_DRY
257 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
258 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
259 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
260# endif
261 real(r8), intent(in) :: alpha(LBi:UBi,LBj:UBj)
262 real(r8), intent(in) :: beta(LBi:UBi,LBj:UBj)
263 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
264 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
265 real(r8), intent(in) :: Hair(LBi:UBi,LBj:UBj)
266 real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
267 real(r8), intent(in) :: Tair(LBi:UBi,LBj:UBj)
268 real(r8), intent(in) :: Uwind(LBi:UBi,LBj:UBj)
269 real(r8), intent(in) :: Vwind(LBi:UBi,LBj:UBj)
270 real(r8), intent(in) :: alpha(LBi:UBi,LBj:UBj)
271 real(r8), intent(in) :: beta(LBi:UBi,LBj:UBj)
272 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
273 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
274 real(r8), intent(in) :: tl_alpha(LBi:UBi,LBj:UBj)
275 real(r8), intent(in) :: tl_beta(LBi:UBi,LBj:UBj)
276 real(r8), intent(in) :: tl_rho(LBi:UBi,LBj:UBj,N(ng))
277 real(r8), intent(in) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
278# ifdef CLOUDS
279 real(r8), intent(in) :: cloud(LBi:UBi,LBj:UBj)
280# endif
281# ifdef BBL_MODEL_NOT_YET
282 real(r8), intent(in) :: Awave(LBi:UBi,LBj:UBj)
283 real(r8), intent(in) :: Pwave(LBi:UBi,LBj:UBj)
284# endif
285 real(r8), intent(in) :: rain(LBi:UBi,LBj:UBj)
286 real(r8), intent(in) :: stflx(LBi:UBi,LBj:UBj,NT(ng))
287
288 real(r8), intent(in) :: lhflx(LBi:UBi,LBj:UBj)
289 real(r8), intent(in) :: lrflx(LBi:UBi,LBj:UBj)
290 real(r8), intent(in) :: shflx(LBi:UBi,LBj:UBj)
291 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
292 real(r8), intent(in) :: stflx(LBi:UBi,LBj:UBj,NT(ng))
293 real(r8), intent(inout) :: tl_lhflx(LBi:UBi,LBj:UBj)
294 real(r8), intent(inout) :: tl_lrflx(LBi:UBi,LBj:UBj)
295 real(r8), intent(inout) :: tl_shflx(LBi:UBi,LBj:UBj)
296 real(r8), intent(inout) :: tl_stflux(LBi:UBi,LBj:UBj,NT(ng))
297
298# ifdef EMINUSP
299 real(r8), intent(out) :: evap(LBi:UBi,LBj:UBj)
300 real(r8), intent(out) :: tl_evap(LBi:UBi,LBj:UBj)
301# endif
302 real(r8), intent(out) :: sustr(LBi:UBi,LBj:UBj)
303 real(r8), intent(out) :: svstr(LBi:UBi,LBj:UBj)
304 real(r8), intent(out) :: tl_sustr(LBi:UBi,LBj:UBj)
305 real(r8), intent(out) :: tl_svstr(LBi:UBi,LBj:UBj)
306# endif
307!
308! Local variable declarations.
309!
310 integer :: Iter, i, j, k
311 integer, parameter :: IterMax = 3
312
313 real(r8), parameter :: eps = 1.0e-20_r8
314 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
315
316 real(r8) :: Bf, Cd, Hl, Hlw, Hscale, Hs, Hsr, IER
317 real(r8) :: tl_Bf, tl_Cd, tl_Hl, tl_Hlw, tl_Hsr, tl_Hs
318 real(r8) :: PairM, RH, Taur
319 real(r8) :: Wspeed, ZQoL, ZToL
320
321 real(r8) :: cff, cff1, cff2, diffh, diffw, oL, upvel
322 real(r8) :: twopi_inv, wet_bulb
323 real(r8) :: tl_wet_bulb, tl_Wspeed, tl_upvel
324 real(r8) :: fac, tl_fac, fac1, tl_fac1, fac2, tl_fac2
325 real(r8) :: fac3, tl_fac3, tl_cff, ad_upvel
326 real(r8) :: adfac, adfac1, adfac2
327# ifdef LONGWAVE
328 real(r8) :: e_sat, vap_p
329# endif
330# ifdef COOL_SKIN
331 real(r8) :: Clam, Fc, Hcool, Hsb, Hlb, Qbouy, Qcool, lambd
332 real(r8) :: tl_Clam, tl_Hcool, tl_Hsb, tl_Hlb
333 real(r8) :: tl_Qbouy, tl_Qcool, tl_lambd
334# endif
335
336 real(r8), dimension(IminS:ImaxS) :: CC
337 real(r8), dimension(IminS:ImaxS) :: Cd10
338 real(r8), dimension(IminS:ImaxS) :: Ch10
339 real(r8), dimension(IminS:ImaxS) :: Ct10
340 real(r8), dimension(IminS:ImaxS) :: charn
341 real(r8), dimension(IminS:ImaxS) :: Ct
342 real(r8), dimension(IminS:ImaxS) :: Cwave
343 real(r8), dimension(IminS:ImaxS) :: Cwet
344 real(r8), dimension(IminS:ImaxS) :: delQ
345 real(r8), dimension(IminS:ImaxS) :: delQc
346 real(r8), dimension(IminS:ImaxS) :: delT
347 real(r8), dimension(IminS:ImaxS) :: delTc
348 real(r8), dimension(IminS:ImaxS) :: delW
349 real(r8), dimension(IminS:ImaxS) :: L
350 real(r8), dimension(IminS:ImaxS) :: L10
351 real(r8), dimension(IminS:ImaxS) :: Lwave
352 real(r8), dimension(IminS:ImaxS) :: Q
353 real(r8), dimension(IminS:ImaxS) :: Qair
354 real(r8), dimension(IminS:ImaxS) :: Qpsi
355 real(r8), dimension(IminS:ImaxS) :: Qsea
356 real(r8), dimension(IminS:ImaxS) :: Qstar
357 real(r8), dimension(IminS:ImaxS) :: rhoAir
358 real(r8), dimension(IminS:ImaxS) :: rhoSea
359 real(r8), dimension(IminS:ImaxS) :: Ri
360 real(r8), dimension(IminS:ImaxS) :: Ribcu
361 real(r8), dimension(IminS:ImaxS) :: Rr
362 real(r8), dimension(IminS:ImaxS) :: Scff
363 real(r8), dimension(IminS:ImaxS) :: TairC
364 real(r8), dimension(IminS:ImaxS) :: TairK
365 real(r8), dimension(IminS:ImaxS) :: Tcff
366 real(r8), dimension(IminS:ImaxS) :: Tpsi
367 real(r8), dimension(IminS:ImaxS) :: TseaC
368 real(r8), dimension(IminS:ImaxS) :: TseaK
369 real(r8), dimension(IminS:ImaxS) :: Tstar
370 real(r8), dimension(IminS:ImaxS) :: u10
371 real(r8), dimension(IminS:ImaxS) :: VisAir
372 real(r8), dimension(IminS:ImaxS) :: Wgus
373 real(r8), dimension(IminS:ImaxS) :: Wmag
374 real(r8), dimension(IminS:ImaxS) :: Wpsi
375 real(r8), dimension(IminS:ImaxS) :: Wstar
376 real(r8), dimension(IminS:ImaxS) :: Zetu
377 real(r8), dimension(IminS:ImaxS) :: Zo10
378 real(r8), dimension(IminS:ImaxS) :: ZoT10
379 real(r8), dimension(IminS:ImaxS) :: ZoL
380 real(r8), dimension(IminS:ImaxS) :: ZoQ
381 real(r8), dimension(IminS:ImaxS) :: ZoT
382 real(r8), dimension(IminS:ImaxS) :: ZoW
383 real(r8), dimension(IminS:ImaxS) :: ZWoL
384
385 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hlv
386 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LHeat
387 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LRad
388 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: SHeat
389 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: SRad
390 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taux
391 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauy
392
393 real(r8), dimension(IminS:ImaxS) :: tl_CC
394 real(r8), dimension(IminS:ImaxS) :: tl_charn
395 real(r8), dimension(IminS:ImaxS) :: tl_Ct
396 real(r8), dimension(IminS:ImaxS) :: tl_Cd10
397 real(r8), dimension(IminS:ImaxS) :: tl_Ct10
398 real(r8), dimension(IminS:ImaxS) :: tl_Cwet
399 real(r8), dimension(IminS:ImaxS) :: tl_delQ
400 real(r8), dimension(IminS:ImaxS) :: tl_delQc
401 real(r8), dimension(IminS:ImaxS) :: tl_delT
402 real(r8), dimension(IminS:ImaxS) :: tl_delTc
403 real(r8), dimension(IminS:ImaxS) :: tl_delW
404 real(r8), dimension(IminS:ImaxS) :: tl_L
405 real(r8), dimension(IminS:ImaxS) :: tl_L10
406 real(r8), dimension(IminS:ImaxS) :: tl_Q
407 real(r8), dimension(IminS:ImaxS) :: tl_Qpsi
408 real(r8), dimension(IminS:ImaxS) :: tl_Qsea
409 real(r8), dimension(IminS:ImaxS) :: tl_Qstar
410 real(r8), dimension(IminS:ImaxS) :: tl_rhoSea
411 real(r8), dimension(IminS:ImaxS) :: tl_Ri
412 real(r8), dimension(IminS:ImaxS) :: tl_Rr
413 real(r8), dimension(IminS:ImaxS) :: tl_Scff
414 real(r8), dimension(IminS:ImaxS) :: tl_Tcff
415 real(r8), dimension(IminS:ImaxS) :: tl_Tpsi
416 real(r8), dimension(IminS:ImaxS) :: tl_TseaC
417 real(r8), dimension(IminS:ImaxS) :: tl_TseaK
418 real(r8), dimension(IminS:ImaxS) :: tl_Tstar
419 real(r8), dimension(IminS:ImaxS) :: tl_u10
420 real(r8), dimension(IminS:ImaxS) :: tl_Wgus
421 real(r8), dimension(IminS:ImaxS) :: tl_Wpsi
422 real(r8), dimension(IminS:ImaxS) :: tl_Wstar
423 real(r8), dimension(IminS:ImaxS) :: tl_Zetu
424 real(r8), dimension(IminS:ImaxS) :: tl_Zo10
425 real(r8), dimension(IminS:ImaxS) :: tl_ZoT10
426 real(r8), dimension(IminS:ImaxS) :: tl_ZoL
427 real(r8), dimension(IminS:ImaxS) :: tl_ZoQ
428 real(r8), dimension(IminS:ImaxS) :: tl_ZoT
429 real(r8), dimension(IminS:ImaxS) :: tl_ZoW
430 real(r8), dimension(IminS:ImaxS) :: tl_ZWoL
431
432 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Hlv
433 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_LHeat
434 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_LRad
435 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_SHeat
436 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Taux
437 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Tauy
438
439# include "set_bounds.h"
440!
441!=======================================================================
442! Atmosphere-Ocean bulk fluxes parameterization.
443!=======================================================================
444!
445 hscale=rho0*cp
446 twopi_inv=0.5_r8/pi
447 DO j=jstr-1,jendr
448 DO i=istr-1,iendr
449!
450! Input bulk parameterization fields.
451!
452 wmag(i)=sqrt(uwind(i,j)*uwind(i,j)+vwind(i,j)*vwind(i,j))
453 pairm=pair(i,j)
454 tairc(i)=tair(i,j)
455 tairk(i)=tairc(i)+273.16_r8
456 tseac(i)=t(i,j,n(ng),nrhs,itemp)
457 tl_tseac(i)=tl_t(i,j,n(ng),nrhs,itemp)
458 tseak(i)=tseac(i)+273.16_r8
459 tl_tseak(i)=tl_tseac(i)
460 rhosea(i)=rho(i,j,n(ng))+1000.0_r8
461 tl_rhosea(i)=tl_rho(i,j,n(ng))
462 rh=hair(i,j)
463 srad(i,j)=srflx(i,j)*hscale
464 tcff(i)=alpha(i,j)
465 tl_tcff(i)=tl_alpha(i,j)
466 scff(i)=beta(i,j)
467 tl_scff(i)=tl_beta(i,j)
468!
469! Initialize.
470!
471 deltc(i)=0.0_r8
472 tl_deltc(i)=0.0_r8
473 delqc(i)=0.0_r8
474 tl_delqc(i)=0.0_r8
475 lheat(i,j)=lhflx(i,j)*hscale
476 tl_lheat(i,j)=tl_lhflx(i,j)*hscale
477 sheat(i,j)=shflx(i,j)*hscale
478 tl_sheat(i,j)=tl_shflx(i,j)*hscale
479 taur=0.0_r8
480 taux(i,j)=0.0_r8
481 tl_taux(i,j)=0.0_r8
482 tauy(i,j)=0.0_r8
483 tl_tauy(i,j)=0.0_r8
484!
485!-----------------------------------------------------------------------
486! Compute net longwave radiation (W/m2), LRad.
487!-----------------------------------------------------------------------
488
489# if defined LONGWAVE
490!
491! Use Berliand (1952) formula to calculate net longwave radiation.
492! The equation for saturation vapor pressure is from Gill (Atmosphere-
493! Ocean Dynamics, pp 606). Here the coefficient in the cloud term
494! is assumed constant, but it is a function of latitude varying from
495! 1.0 at poles to 0.5 at the Equator).
496!
497 cff=(0.7859_r8+0.03477_r8*tairc(i))/ &
498 & (1.0_r8+0.00412_r8*tairc(i))
499 e_sat=10.0_r8**cff ! saturation vapor pressure (hPa or mbar)
500 vap_p=e_sat*rh ! water vapor pressure (hPa or mbar)
501 cff2=tairk(i)*tairk(i)*tairk(i)
502 cff1=cff2*tairk(i)
503 lrad(i,j)=-emmiss*stefbo* &
504 & (cff1*(0.39_r8-0.05_r8*sqrt(vap_p))* &
505 & (1.0_r8-0.6823_r8*cloud(i,j)*cloud(i,j))+ &
506 & cff2*4.0_r8*(tseak(i)-tairk(i)))
507 tl_lrad(i,j)=-emmiss*stefbo*cff2*4.0_r8*tl_tseak(i)
508
509# elif defined LONGWAVE_OUT
510!
511! Treat input longwave data as downwelling radiation only and add
512! outgoing IR from model sea surface temperature.
513!
514 lrad(i,j)=lrflx(i,j)*hscale- &
515 & emmiss*stefbo*tseak(i)*tseak(i)*tseak(i)*tseak(i)
516 tl_lrad(i,j)=tl_lrflx(i,j)*hscale- &
517 & 4.0_r8*emmiss*stefbo* &
518 & tl_tseak(i)*tseak(i)*tseak(i)*tseak(i)
519# else
520 lrad(i,j)=lrflx(i,j)*hscale
521 tl_lrad(i,j)=tl_lrflx(i,j)*hscale
522# endif
523# ifdef MASKING
524 lrad(i,j)=lrad(i,j)*rmask(i,j)
525 tl_lrad(i,j)=tl_lrad(i,j)*rmask(i,j)
526# endif
527# ifdef WET_DRY
528 lrad(i,j)=lrad(i,j)*rmask_wet(i,j)
529 tl_lrad(i,j)=tl_lrad(i,j)*rmask_wet(i,j)
530# endif
531!
532!-----------------------------------------------------------------------
533! Compute specific humidities (kg/kg).
534!
535! note that Qair(i) is the saturation specific humidity at Tair
536! Q(i) is the actual specific humidity
537! Qsea(i) is the saturation specific humidity at Tsea
538!
539! Saturation vapor pressure in mb is first computed and then
540! converted to specific humidity in kg/kg
541!
542! The saturation vapor pressure is computed from Teten formula
543! using the approach of Buck (1981):
544!
545! Esat(mb) = (1.0007_r8+3.46E-6_r8*PairM(mb))*6.1121_r8*
546! EXP(17.502_r8*TairC(C)/(240.97_r8+TairC(C)))
547!
548! The ambient vapor is found from the definition of the
549! Relative humidity:
550!
551! RH = W/Ws*100 ~ E/Esat*100 E = RH/100*Esat if RH is in %
552! E = RH*Esat if RH fractional
553!
554! The specific humidity is then found using the relationship:
555!
556! Q = 0.622 E/(P + (0.622-1)e)
557!
558! Q(kg/kg) = 0.62197_r8*(E(mb)/(PairM(mb)-0.378_r8*E(mb)))
559!
560!-----------------------------------------------------------------------
561!
562! Compute air saturation vapor pressure (mb), using Teten formula.
563!
564 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
565 & exp(17.502_r8*tairc(i)/(240.97_r8+tairc(i)))
566!
567! Compute specific humidity at Saturation, Qair (kg/kg).
568!
569 qair(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
570!
571! Compute specific humidity, Q (kg/kg).
572!
573 IF (rh.lt.2.0_r8) THEN !RH fraction
574 cff=cff*rh !Vapor pres (mb)
575 q(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff)) !Spec hum (kg/kg)
576 ELSE !RH input was actually specific humidity in g/kg
577 q(i)=rh/1000.0_r8 !Spec Hum (kg/kg)
578 END IF
579!
580! Compute water saturation vapor pressure (mb), using Teten formula.
581!
582 fac=17.502_r8*tseac(i)/(240.97_r8+tseac(i))
583 tl_fac=17.502_r8*tl_tseac(i)/(240.97_r8+tseac(i))- &
584 & tl_tseac(i)*fac/(240.97_r8+tseac(i))
585!^ cff=(1.0007_r8+3.46E-6_r8*PairM)*6.1121_r8*EXP(fac)
586!^
587 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
588 & exp(17.502_r8*tseac(i)/(240.97_r8+tseac(i)))
589 tl_cff=tl_fac*cff
590!
591! Vapor Pressure reduced for salinity (Kraus and Businger, 1994, pp42).
592!
593 cff=cff*0.98_r8
594 tl_cff=tl_cff*0.98_r8
595!
596! Compute Qsea (kg/kg) from vapor pressure.
597!
598 qsea(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
599 tl_fac=0.62197_r8*(tl_cff/(pairm-0.378_r8*cff))
600 tl_qsea(i)=tl_fac*(1.0_r8+0.378_r8*cff/((pairm-0.378_r8*cff)))
601!
602!-----------------------------------------------------------------------
603! Compute Monin-Obukhov similarity parameters for wind (Wstar),
604! heat (Tstar), and moisture (Qstar), Liu et al. (1979).
605!-----------------------------------------------------------------------
606!
607! Moist air density (kg/m3).
608!
609 rhoair(i)=pairm*100.0_r8/(blk_rgas*tairk(i)* &
610 & (1.0_r8+0.61_r8*q(i)))
611!
612! Kinematic viscosity of dry air (m2/s), Andreas (1989).
613!
614 visair(i)=1.326e-5_r8* &
615 & (1.0_r8+tairc(i)*(6.542e-3_r8+tairc(i)* &
616 & (8.301e-6_r8-4.84e-9_r8*tairc(i))))
617!
618! Compute latent heat of vaporization (J/kg) at sea surface, Hlv.
619!
620 hlv(i,j)=(2.501_r8-0.00237_r8*tseac(i))*1.0e+6_r8
621 tl_hlv(i,j)=-0.00237_r8*tl_tseac(i)*1.0e+6_r8
622!
623! Assume that wind is measured relative to sea surface and include
624! gustiness.
625!
626 wgus(i)=0.5_r8
627!
628! Note: tl_ZoW is zero at this point.
629!
630 tl_wgus(i)=0.0_r8
631 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
632 tl_delw(i)=0.0_r8
633 delq(i)=qsea(i)-q(i)
634 tl_delq(i)=tl_qsea(i)
635 delt(i)=tseac(i)-tairc(i)
636 tl_delt(i)=tl_tseac(i)
637!
638! Neutral coefficients.
639!
640 zow(i)=0.0001_r8
641 tl_zow(i)=0.0_r8
642 u10(i)=delw(i)*log(10.0_r8/zow(i))/log(blk_zw(ng)/zow(i))
643!
644! Note tl_delW and tl_ZoW are zero at this point.
645!
646 tl_u10(i)=0.0_r8
647 wstar(i)=0.035_r8*u10(i)
648 tl_wstar(i)=0.035_r8*tl_u10(i)
649 zo10(i)=0.011_r8*wstar(i)*wstar(i)/g+ &
650 & 0.11_r8*visair(i)/wstar(i)
651 tl_zo10(i)=0.022_r8*tl_wstar(i)*wstar(i)/g- &
652 & tl_wstar(i)*0.11_r8*visair(i)/(wstar(i)*wstar(i))
653 fac=log(10.0_r8/zo10(i))
654 tl_fac=-tl_zo10(i)/zo10(i)
655!^ Cd10(i)=(vonKar/fac)**2
656!^
657 cd10(i)=(vonkar/log(10.0_r8/zo10(i)))**2
658 tl_cd10(i)=-2.0_r8*tl_fac*cd10(i)/fac
659 ch10(i)=0.00115_r8
660 ct10(i)=ch10(i)/sqrt(cd10(i))
661 tl_ct10(i)=-0.5_r8*tl_cd10(i)*ct10(i)/cd10(i)
662 fac=vonkar/ct10(i)
663 tl_fac=-tl_ct10(i)*fac/ct10(i)
664!^ ZoT10(i)=10.0_r8/EXP(fac)
665!^
666 zot10(i)=10.0_r8/exp(vonkar/ct10(i))
667 tl_zot10(i)=-tl_fac*zot10(i)
668 fac=log(blk_zw(ng)/zo10(i))
669 tl_fac=-tl_zo10(i)/zo10(i)
670!^ Cd=(vonKar/fac)**2
671!^
672 cd=(vonkar/log(blk_zw(ng)/zo10(i)))**2
673 tl_cd=-2.0_r8*tl_fac*cd/fac
674!
675! Compute Richardson number.
676!
677 fac=log(blk_zt(ng)/zot10(i))
678 tl_fac=-tl_zot10(i)/zot10(i)
679!^ Ct(i)=vonKar/fac ! T transfer coefficient
680!^
681 ct(i)=vonkar/log(blk_zt(ng)/zot10(i)) ! T transfer coefficient
682 tl_ct(i)=-tl_fac*ct(i)/fac ! T transfer coefficient
683 cc(i)=vonkar*ct(i)/cd
684 tl_cc(i)=vonkar*tl_ct(i)/cd-tl_cd*cc(i)/cd
685 deltc(i)=0.0_r8
686 tl_deltc(i)=0.0_r8
687# ifdef COOL_SKIN
688 deltc(i)=blk_dter
689 tl_deltc(i)=0.0_r8
690# endif
691 ribcu(i)=-blk_zw(ng)/(blk_zabl*0.004_r8*blk_beta**3)
692 fac=1.0/(tairk(i)*delw(i)*delw(i))
693!
694! Note tl_delW is zero at this point.
695!
696!^ Ri(i)=-fac*g*blk_ZW(ng)*((delT(i)-delTc(i))+ &
697!^ & 0.61_r8*TairK(i)*delQ(i))
698!^
699 ri(i)=-g*blk_zw(ng)*((delt(i)-deltc(i))+ &
700 & 0.61_r8*tairk(i)*delq(i))/ &
701 & (tairk(i)*delw(i)*delw(i))
702 tl_ri(i)=-fac*g*blk_zw(ng)*((tl_delt(i)-tl_deltc(i))+ &
703 & 0.61_r8*tairk(i)*tl_delq(i))
704 IF (ri(i).lt.0.0_r8) THEN
705 zetu(i)=cc(i)*ri(i)/(1.0_r8+ri(i)/ribcu(i)) ! Unstable
706 tl_zetu(i)=(tl_cc(i)*ri(i)+cc(i)*tl_ri(i))/ &
707 & (1.0_r8+ri(i)/ribcu(i))- &
708 & (tl_ri(i)/ribcu(i))* &
709 & zetu(i)/(1.0_r8+ri(i)/ribcu(i))
710 ELSE
711 fac=3.0_r8*ri(i)/cc(i)
712 tl_fac=3.0_r8*tl_ri(i)/cc(i)-tl_cc(i)*fac/cc(i)
713!^ Zetu(i)=CC(i)*Ri(i)/(1.0_r8+fac) ! Stable
714!^
715 zetu(i)=cc(i)*ri(i)/(1.0_r8+3.0_r8*ri(i)/cc(i)) ! Stable
716 tl_zetu(i)=(tl_cc(i)*ri(i)+cc(i)*tl_ri(i))/(1.0_r8+fac)- &
717 & tl_fac*zetu(i)/(1.0_r8+fac)
718 END IF
719 l10(i)=blk_zw(ng)/zetu(i)
720 tl_l10(i)=-l10(i)*l10(i)*tl_zetu(i)/blk_zw(ng)
721!
722! First guesses for Monon-Obukhov similarity scales.
723!
724 fac1=blk_zw(ng)/l10(i)
725 tl_fac1=-tl_l10(i)*fac1/l10(i)
726 fac2=log(blk_zw(ng)/zo10(i))
727 tl_fac2=-tl_zo10(i)/zo10(i)
728!
729! Note tl_delW is zero at this point.
730!
731!^ Wstar(i)=delW(i)*vonKar/(fac2-bulk_psiu(fac1,pi))
732!^
733 wstar(i)=delw(i)*vonkar/(log(blk_zw(ng)/zo10(i))- &
734 & bulk_psiu(blk_zw(ng)/l10(i),pi))
735 tl_wstar(i)=-(tl_fac2- &
736 & tl_bulk_psiu(tl_fac1,fac1,pi))*wstar(i)/ &
737 & (fac2-bulk_psiu(fac1,pi))
738 fac1=blk_zt(ng)/l10(i)
739 tl_fac1=-tl_l10(i)*fac1/l10(i)
740 fac2=log(blk_zt(ng)/zot10(i))
741 tl_fac2=-tl_zot10(i)/zot10(i)
742!^ Tstar(i)=-(delT(i)-delTc(i))*vonKar/ &
743!^ & (fac2-bulk_psit(fac1,pi))
744!^
745 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
746 & (log(blk_zt(ng)/zot10(i))- &
747 & bulk_psit(blk_zt(ng)/l10(i),pi))
748 tl_tstar(i)=-(tl_delt(i)-tl_deltc(i))*vonkar/ &
749 & (fac2-bulk_psit(fac1,pi))- &
750 & (tl_fac2-tl_bulk_psit(tl_fac1,fac1,pi))*tstar(i)/ &
751 & (fac2-bulk_psit(fac1,pi))
752 fac1=blk_zq(ng)/l10(i)
753 tl_fac1=-tl_l10(i)*fac1/l10(i)
754 fac2=log(blk_zq(ng)/zot10(i))
755 tl_fac2=-tl_zot10(i)/zot10(i)
756!^ Qstar(i)=-(delQ(i)-delQc(i))*vonKar/ &
757!^ & (fac2-bulk_psit(fac1,pi))
758 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
759 & (log(blk_zq(ng)/zot10(i))- &
760 & bulk_psit(blk_zq(ng)/l10(i),pi))
761 tl_qstar(i)=-(tl_delq(i)-tl_delqc(i))*vonkar/ &
762 & (fac2-bulk_psit(fac1,pi))- &
763 & (tl_fac2-tl_bulk_psit(tl_fac1,fac1,pi))*qstar(i)/ &
764 & (fac2-bulk_psit(fac1,pi))
765!
766! Modify Charnock for high wind speeds. The 0.125 factor below is for
767! 1.0/(18.0-10.0).
768!
769 IF (delw(i).gt.18.0_r8) THEN
770 charn(i)=0.018_r8
771 tl_charn(i)=0.0_r8
772 ELSE IF ((10.0_r8.lt.delw(i)).and.(delw(i).le.18.0_r8)) THEN
773 charn(i)=0.011_r8+ &
774 & 0.125_r8*(0.018_r8-0.011_r8)*(delw(i)-10._r8)
775!
776! Note tl_delW is zero at this point.
777!
778 tl_charn(i)=0.0_r8
779 ELSE
780 charn(i)=0.011_r8
781 tl_charn(i)=0.0_r8
782 END IF
783# ifdef BBL_MODEL_NOT_YET
784 cwave(i)=g*pwave(i,j)*twopi_inv
785 lwave(i)=cwave(i)*pwave(i,j)
786# endif
787 END DO
788!
789! Iterate until convergence. It usually converges within 3 iterations.
790!
791 DO iter=1,itermax
792 DO i=istr-1,iendr
793# ifdef BBL_MODEL_NOT_YET
794!
795! Use wave info if we have it, two different options.
796!
797# ifdef WIND_WAVES
798 zow(i)=(25.0_r8/pi)*lwave(i)*(wstar(i)/cwave(i))**4.5+ &
799 & 0.11_r8*visair(i)/(wstar(i)+eps)
800 tl_zow(i)=(25.0_r8/pi)*lwave(i)*4.5_r8*tl_wstar(i)* &
801 & (wstar(i)/cwave(i))**3.5/cwave(i)- &
802 & tl_wstar(i)*0.11_r8*visair(i)/ &
803 & ((wstar(i)+eps)*(wstar(i)+eps))
804# else
805 zow(i)=1200._r8*awave(i,j)*(awave(i,j)/lwave(i))**4.5+ &
806 & 0.11_r8*visair(i)/(wstar(i)+eps)
807 tl_zow(i)=-tl_wstar(i)*0.11_r8*visair(i)/ &
808 & ((wstar(i)+eps)*(wstar(i)+eps))
809# endif
810# else
811 zow(i)=charn(i)*wstar(i)*wstar(i)/g+ &
812 & 0.11_r8*visair(i)/(wstar(i)+eps)
813 tl_zow(i)=2.0_r8*charn(i)*tl_wstar(i)*wstar(i)/g- &
814 & tl_wstar(i)*0.11_r8*visair(i)/ &
815 & ((wstar(i)+eps)*(wstar(i)+eps))
816# endif
817 rr(i)=zow(i)*wstar(i)/visair(i)
818 tl_rr(i)=(tl_zow(i)*wstar(i)+zow(i)*tl_wstar(i))/visair(i)
819!
820! Compute Monin-Obukhov stability parameter, Z/L.
821!
822 zoq(i)=min(1.15e-4_r8,5.5e-5_r8/rr(i)**0.6)
823 tl_zoq(i)= &
824 & -(0.5_r8-sign(0.5_r8,5.5e-5_r8/rr(i)**0.6-1.15e-4_r8)) &
825 & *0.6_r8*5.5e-5_r8*tl_rr(i)/rr(i)**1.6
826 zot(i)=zoq(i)
827 tl_zot(i)=tl_zoq(i)
828 zol(i)=vonkar*g*blk_zw(ng)* &
829 & (tstar(i)*(1.0_r8+0.61_r8*q(i))+ &
830 & 0.61_r8*tairk(i)*qstar(i))/ &
831 & (tairk(i)*wstar(i)*wstar(i)* &
832 & (1.0_r8+0.61_r8*q(i))+eps)
833 tl_zol(i)=vonkar*g*blk_zw(ng)* &
834 & (tl_tstar(i)*(1.0_r8+0.61_r8*q(i))+ &
835 & 0.61_r8*tairk(i)*tl_qstar(i))/ &
836 & (tairk(i)*wstar(i)*wstar(i)* &
837 & (1.0_r8+0.61_r8*q(i))+eps)- &
838 & 2.0_r8*tairk(i)*tl_wstar(i)*wstar(i)* &
839 & (1.0_r8+0.61_r8*q(i))*zol(i)/ &
840 & (tairk(i)*wstar(i)*wstar(i)* &
841 & (1.0_r8+0.61_r8*q(i))+eps)
842 l(i)=blk_zw(ng)/(zol(i)+eps)
843 tl_l(i)=-tl_zol(i)*l(i)/(zol(i)+eps)
844!
845! Evaluate stability functions at Z/L.
846!
847 wpsi(i)=bulk_psiu(zol(i),pi)
848 tl_wpsi(i)=tl_bulk_psiu(tl_zol(i),zol(i),pi)
849 fac=blk_zt(ng)/l(i)
850 tl_fac=-tl_l(i)*fac/l(i)
851!^ Tpsi(i)=bulk_psit(fac,pi)
852!^
853 tpsi(i)=bulk_psit(blk_zt(ng)/l(i),pi)
854 tl_tpsi(i)=tl_bulk_psit(tl_fac,fac,pi)
855 fac=blk_zq(ng)/l(i)
856 tl_fac=-tl_l(i)*fac/l(i)
857!^ Qpsi(i)=bulk_psit(fac,pi)
858!^
859 qpsi(i)=bulk_psit(blk_zq(ng)/l(i),pi)
860 tl_qpsi(i)=tl_bulk_psit(tl_fac,fac,pi)
861# ifdef COOL_SKIN
862 cwet(i)=0.622_r8*hlv(i,j)*qsea(i)/ &
863 & (blk_rgas*tseak(i)*tseak(i))
864 tl_cwet(i)=0.622_r8*(tl_hlv(i,j)*qsea(i)+ &
865 & hlv(i,j)*tl_qsea(i))/ &
866 & (blk_rgas*tseak(i)*tseak(i))- &
867 & 2.0_r8*blk_rgas*tl_tseak(i)*tseak(i)*cwet(i)/ &
868 & (blk_rgas*tseak(i)*tseak(i))
869 delqc(i)=cwet(i)*deltc(i)
870 tl_delqc(i)=tl_cwet(i)*deltc(i)+cwet(i)*tl_deltc(i)
871# endif
872!
873! Compute wind scaling parameters, Wstar.
874!
875 fac1=blk_zw(ng)/zow(i)
876 tl_fac1=-tl_zow(i)*fac1/zow(i)
877 fac2=log(fac1)
878 tl_fac2=tl_fac1/fac1
879 fac3=delw(i)*vonkar/(fac2-wpsi(i))
880 tl_fac3=tl_delw(i)*vonkar/(fac2-wpsi(i))- &
881 & (tl_fac2-tl_wpsi(i))*fac3/(fac2-wpsi(i))
882!^ Wstar(i)=MAX(eps,fac3)
883!^
884 wstar(i)=max(eps,delw(i)*vonkar/ &
885 & (log(blk_zw(ng)/zow(i))-wpsi(i)))
886 tl_wstar(i)=(0.5_r8-sign(0.5_r8,eps-fac3))*tl_fac3
887 fac1=blk_zt(ng)/zot(i)
888 tl_fac1=-tl_zot(i)*fac1/zot(i)
889 fac2=log(fac1)
890 tl_fac2=tl_fac1/fac1
891!^ Tstar(i)=-(delT(i)-delTc(i))*vonKar/(fac2-Tpsi(i))
892!^
893 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
894 & (log(blk_zt(ng)/zot(i))-tpsi(i))
895 tl_tstar(i)=-(tl_delt(i)-tl_deltc(i))*vonkar/(fac2-tpsi(i))-&
896 & (tl_fac2-tl_tpsi(i))*tstar(i)/(fac2-tpsi(i))
897 fac1=blk_zq(ng)/zoq(i)
898 tl_fac1=-tl_zoq(i)*fac1/zoq(i)
899 fac2=log(fac1)
900 tl_fac2=tl_fac1/fac1
901!^ Qstar(i)=-(delQ(i)-delQc(i))*vonKar/(fac2-Qpsi(i))
902!^
903 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
904 & (log(blk_zq(ng)/zoq(i))-qpsi(i))
905 tl_qstar(i)=-(tl_delq(i)-tl_delqc(i))*vonkar/(fac2-qpsi(i))-&
906 & (tl_fac2-tl_qpsi(i))*qstar(i)/(fac2-qpsi(i))
907
908!
909! Compute gustiness in wind speed.
910!
911 bf=-g/tairk(i)* &
912 & wstar(i)*(tstar(i)+0.61_r8*tairk(i)*qstar(i))
913 tl_bf=-g/tairk(i)* &
914 & (tl_wstar(i)*(tstar(i)+0.61_r8*tairk(i)*qstar(i))+ &
915 & wstar(i)*(tl_tstar(i)+0.61_r8*tairk(i)*tl_qstar(i)))
916 IF (bf.gt.0.0_r8) THEN
917 wgus(i)=blk_beta*(bf*blk_zabl)**r3
918 tl_wgus(i)=r3*blk_beta*tl_bf*blk_zabl* &
919 & (bf*blk_zabl)**(r3-1.0_r8)
920 ELSE
921 wgus(i)=0.2_r8
922 tl_wgus(i)=0.0_r8
923 END IF
924 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
925 IF (delw(i).ne.0.0_r8)THEN
926 tl_delw(i)=tl_wgus(i)*wgus(i)/delw(i)
927 ELSE
928 tl_delw(i)=0.0_r8
929 END IF
930
931# ifdef COOL_SKIN
932!
933!-----------------------------------------------------------------------
934! Cool Skin correction.
935!-----------------------------------------------------------------------
936!
937! Cool skin correction constants. Clam: part of Saunders constant
938! lambda; Cwet: slope of saturation vapor.
939!
940 clam=16.0_r8*g*blk_cpw*(rhosea(i)*blk_visw)**3/ &
941 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
942 tl_clam=48.0_r8*g*blk_cpw*tl_rhosea(i)*blk_visw* &
943 & (rhosea(i)*blk_visw)**2/ &
944 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
945!
946! Set initial guesses for cool-skin layer thickness (Hcool).
947!
948 hcool=0.001_r8
949!
950! Backgound sensible and latent heat.
951!
952 hsb=-rhoair(i)*blk_cpa*wstar(i)*tstar(i)
953 tl_hsb=-rhoair(i)*blk_cpa*(tl_wstar(i)*tstar(i)+ &
954 & wstar(i)*tl_tstar(i))
955 hlb=-rhoair(i)*hlv(i,j)*wstar(i)*qstar(i)
956 tl_hlb=-rhoair(i)*(tl_hlv(i,j)*wstar(i)*qstar(i)+ &
957 & hlv(i,j)*(tl_wstar(i)*qstar(i)+ &
958 & wstar(i)*tl_qstar(i)))
959!
960! Mean absoption in cool-skin layer.
961!
962 fc=0.065_r8+11.0_r8*hcool- &
963 & (1.0_r8-exp(-hcool*1250.0_r8))*6.6e-5_r8/hcool
964!
965! Total cooling at the interface.
966!
967 qcool=lrad(i,j)+hsb+hlb-srad(i,j)*fc
968 tl_qcool=tl_lrad(i,j)+tl_hsb+tl_hlb
969 qbouy=tcff(i)*qcool+scff(i)*hlb*blk_cpw/hlv(i,j)
970 tl_qbouy=tl_tcff(i)*qcool+tcff(i)*tl_qcool+ &
971 & (tl_scff(i)*hlb*blk_cpw+ &
972 & scff(i)*tl_hlb*blk_cpw)/hlv(i,j)- &
973 & tl_hlv(i,j)*scff(i)*hlb*blk_cpw/ &
974 & (hlv(i,j)*hlv(i,j))
975!
976! Compute temperature and moisture change.
977!
978 IF ((qcool.gt.0.0_r8).and.(qbouy.gt.0.0_r8)) THEN
979 fac1=(wstar(i)+eps)**4
980 tl_fac1=4.0_r8*tl_wstar(i)*(wstar(i)+eps)**3
981 fac2=clam*qbouy
982 tl_fac2=tl_clam*qbouy+clam*tl_qbouy
983 fac3=(fac2/fac1)**0.75_r8
984 tl_fac3=0.75_r8*(fac2/fac1)**(0.75_r8-1.0_r8)* &
985 & (tl_fac2/fac1-tl_fac1*fac2/(fac1*fac1))
986!^ lambd=6.0_r8/(1.0_r8+fac3)**r3
987!^
988 lambd=6.0_r8/ &
989 & (1.0_r8+(clam*qbouy/(wstar(i)+eps)**4)**0.75_r8)**r3
990 tl_lambd=-r3*6.0_r8*tl_fac3/(1.0_r8+fac3)**(r3+1.0_r8)
991 fac1=sqrt(rhoair(i)/rhosea(i))
992 IF (fac1.ne.0.0_r8) THEN
993 tl_fac1=-0.5_r8*tl_rhosea(i)*fac1/rhosea(i)
994 ELSE
995 tl_fac1=0.0_r8
996 END IF
997 fac2=fac1*wstar(i)+eps
998 tl_fac2=tl_fac1*wstar(i)+fac1*tl_wstar(i)
999!^ Hcool=lambd*blk_visw/fac2
1000!^
1001 hcool=lambd*blk_visw/(sqrt(rhoair(i)/rhosea(i))* &
1002 & wstar(i)+eps)
1003 tl_hcool=tl_lambd*blk_visw/fac2-tl_fac2*hcool/fac2
1004 deltc(i)=qcool*hcool/blk_tcw
1005 tl_deltc(i)=(tl_qcool*hcool+qcool*tl_hcool)/blk_tcw
1006 ELSE
1007 deltc(i)=0.0_r8
1008 tl_deltc(i)=0.0_r8
1009 END IF
1010 delqc(i)=cwet(i)*deltc(i)
1011 tl_delqc(i)=tl_cwet(i)*deltc(i)+cwet(i)*tl_deltc(i)
1012# endif
1013 END DO
1014 END DO
1015!
1016!-----------------------------------------------------------------------
1017! Compute Atmosphere/Ocean fluxes.
1018!-----------------------------------------------------------------------
1019!
1020 DO i=istr-1,iendr
1021!
1022! Compute transfer coefficients for momentum (Cd).
1023!
1024 wspeed=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
1025 IF (wspeed.ne.0.0_r8) THEN
1026 tl_wspeed=tl_wgus(i)*wgus(i)/wspeed
1027 ELSE
1028 tl_wspeed=0.0_r8
1029 END IF
1030 cd=wstar(i)*wstar(i)/(wspeed*wspeed+eps)
1031 tl_cd=2.0_r8*(tl_wstar(i)*wstar(i)/(wspeed*wspeed+eps)- &
1032 & tl_wspeed*wspeed*cd/(wspeed*wspeed+eps))
1033!
1034! Compute turbulent sensible heat flux (W/m2), Hs.
1035!
1036 hs=-blk_cpa*rhoair(i)*wstar(i)*tstar(i)
1037 tl_hs=-blk_cpa*rhoair(i)*(tl_wstar(i)*tstar(i)+ &
1038 & wstar(i)*tl_tstar(i))
1039!
1040! Compute sensible heat flux (W/m2) due to rainfall (kg/m2/s), Hsr.
1041!
1042 diffw=2.11e-5_r8*(tairk(i)/273.16_r8)**1.94_r8
1043 diffh=0.02411_r8*(1.0_r8+tairc(i)* &
1044 & (3.309e-3_r8-1.44e-6_r8*tairc(i)))/ &
1045 & (rhoair(i)*blk_cpa)
1046 cff=qair(i)*hlv(i,j)/(blk_rgas*tairk(i)*tairk(i))
1047 tl_cff=qair(i)*tl_hlv(i,j)/(blk_rgas*tairk(i)*tairk(i))
1048 fac=0.622_r8*(cff*hlv(i,j)*diffw)/(blk_cpa*diffh)
1049 tl_fac=0.622_r8*diffw*(tl_cff*hlv(i,j)+cff*tl_hlv(i,j))/ &
1050 & (blk_cpa*diffh)
1051!^ wet_bulb=1.0_r8/(1.0_r8+fac)
1052!^
1053 wet_bulb=1.0_r8/(1.0_r8+0.622_r8*(cff*hlv(i,j)*diffw)/ &
1054 & (blk_cpa*diffh))
1055 tl_wet_bulb=-tl_fac*wet_bulb*wet_bulb
1056 hsr=rain(i,j)*wet_bulb*blk_cpw* &
1057 & ((tseac(i)-tairc(i))+(qsea(i)-q(i))*hlv(i,j)/blk_cpa)
1058 tl_hsr=hsr*tl_wet_bulb/wet_bulb+ &
1059 & rain(i,j)*wet_bulb*blk_cpw* &
1060 & (tl_tseac(i)+tl_qsea(i)*hlv(i,j)/blk_cpa+ &
1061 & (qsea(i)-q(i))*tl_hlv(i,j)/blk_cpa)
1062 sheat(i,j)=(hs+hsr)
1063 tl_sheat(i,j)=(tl_hs+tl_hsr)
1064# ifdef MASKING
1065 sheat(i,j)=sheat(i,j)*rmask(i,j)
1066 tl_sheat(i,j)=tl_sheat(i,j)*rmask(i,j)
1067# endif
1068# ifdef WET_DRY
1069 sheat(i,j)=sheat(i,j)*rmask_wet(i,j)
1070 tl_sheat(i,j)=tl_sheat(i,j)*rmask_wet(i,j)
1071# endif
1072!
1073! Compute turbulent latent heat flux (W/m2), Hl.
1074!
1075 hl=-hlv(i,j)*rhoair(i)*wstar(i)*qstar(i)
1076 tl_hl=-tl_hlv(i,j)*rhoair(i)*wstar(i)*qstar(i)- &
1077 & hlv(i,j)*rhoair(i)*(tl_wstar(i)*qstar(i)+ &
1078 & wstar(i)*tl_qstar(i) )
1079!
1080! Compute Webb correction (Webb effect) to latent heat flux, Hlw.
1081!
1082 upvel=-1.61_r8*wstar(i)*qstar(i)- &
1083 & (1.0_r8+1.61_r8*q(i))*wstar(i)*tstar(i)/tairk(i)
1084 tl_upvel=-1.61_r8* &
1085 & (tl_wstar(i)*qstar(i)+wstar(i)*tl_qstar(i))- &
1086 & (1.0_r8+1.61_r8*q(i))*(tl_wstar(i)*tstar(i)+ &
1087 & wstar(i)*tl_tstar(i))/ &
1088 & tairk(i)
1089 hlw=rhoair(i)*hlv(i,j)*upvel*q(i)
1090 tl_hlw=rhoair(i)*q(i)*(tl_hlv(i,j)*upvel+ &
1091 & hlv(i,j)*tl_upvel)
1092 lheat(i,j)=(hl+hlw)
1093 tl_lheat(i,j)=(tl_hl+tl_hlw)
1094# ifdef MASKING
1095 lheat(i,j)=lheat(i,j)*rmask(i,j)
1096 tl_lheat(i,j)=tl_lheat(i,j)*rmask(i,j)
1097# endif
1098# ifdef WET_DRY
1099 lheat(i,j)=lheat(i,j)*rmask_wet(i,j)
1100 tl_lheat(i,j)=tl_lheat(i,j)*rmask_wet(i,j)
1101# endif
1102!
1103! Compute momentum flux (N/m2) due to rainfall (kg/m2/s).
1104!
1105 taur=0.85_r8*rain(i,j)*wmag(i)
1106!
1107! Compute wind stress components (N/m2), Tau.
1108!
1109 cff=rhoair(i)*cd*wspeed
1110 tl_cff=rhoair(i)*(tl_cd*wspeed+cd*tl_wspeed)
1111 taux(i,j)=(cff*uwind(i,j)+taur*sign(1.0_r8,uwind(i,j)))
1112 tl_taux(i,j)=tl_cff*uwind(i,j)
1113# ifdef MASKING
1114 taux(i,j)=taux(i,j)*rmask(i,j)
1115 tl_taux(i,j)=tl_taux(i,j)*rmask(i,j)
1116# endif
1117# ifdef WET_DRY
1118 taux(i,j)=taux(i,j)*rmask_wet(i,j)
1119 tl_taux(i,j)=tl_taux(i,j)*rmask_wet(i,j)
1120# endif
1121 tauy(i,j)=(cff*vwind(i,j)+taur*sign(1.0_r8,vwind(i,j)))
1122 tl_tauy(i,j)=tl_cff*vwind(i,j)
1123# ifdef MASKING
1124 tauy(i,j)=tauy(i,j)*rmask(i,j)
1125 tl_tauy(i,j)=tl_tauy(i,j)*rmask(i,j)
1126# endif
1127# ifdef WET_DRY
1128 tauy(i,j)=tauy(i,j)*rmask_wet(i,j)
1129 tl_tauy(i,j)=tl_tauy(i,j)*rmask_wet(i,j)
1130# endif
1131 END DO
1132 END DO
1133!
1134!=======================================================================
1135! Compute surface net heat flux and surface wind stress.
1136!=======================================================================
1137!
1138! Compute kinematic, surface, net heat flux (degC m/s). Notice that
1139! the signs of latent and sensible fluxes are reversed because fluxes
1140! calculated from the bulk formulations above are positive out of the
1141! ocean.
1142!
1143! For EMINUSP option, EVAP = LHeat (W/m2) / Hlv (J/kg) = kg/m2/s
1144! PREC = rain = kg/m2/s
1145!
1146! To convert these rates to m/s divide by freshwater density, rhow.
1147!
1148! Note that when the air is undersaturated in water vapor (Q < Qsea)
1149! the model will evaporate and LHeat > 0:
1150!
1151! LHeat positive out of the ocean
1152! evap positive out of the ocean
1153!
1154! Note that if evaporating, the salt flux is positive
1155! and if raining, the salt flux is negative
1156!
1157! Note that stflux(:,:,isalt) is the E-P flux. The fresh water flux
1158! is positive out of the ocean and the salt flux is positive into the
1159! ocean. It is multiplied by surface salinity when computing state
1160! variable stflx(:,:,isalt) in "set_vbc.F".
1161!
1162 hscale=1.0_r8/(rho0*cp)
1163 cff=1.0_r8/rhow
1164 DO j=jstrr,jendr
1165 DO i=istrr,iendr
1166!^ lrflx(i,j)=LRad(i,j)*Hscale
1167!^
1168 tl_lrflx(i,j)=tl_lrad(i,j)*hscale
1169!^ lhflx(i,j)=-LHeat(i,j)*Hscale
1170!^
1171 tl_lhflx(i,j)=-tl_lheat(i,j)*hscale
1172!^ shflx(i,j)=-SHeat(i,j)*Hscale
1173!^
1174 tl_shflx(i,j)=-tl_sheat(i,j)*hscale
1175!^ stflux(i,j,itemp)=(srflx(i,j)+lrflx(i,j)+ &
1176!^ & lhflx(i,j)+shflx(i,j))
1177!^
1178 tl_stflux(i,j,itemp)=(tl_lrflx(i,j)+ &
1179 & tl_lhflx(i,j)+tl_shflx(i,j))
1180# ifdef MASKING
1181!^ stflux(i,j,itemp)=stflx(i,j,itemp)*rmask(i,j)
1182!^
1183 tl_stflx(i,j,itemp)=tl_stflx(i,j,itemp)*rmask(i,j)
1184# endif
1185# ifdef WET_DRY
1186!^ stflux(i,j,itemp)=stflx(i,j,itemp)*rmask_wet(i,j)
1187!^
1188 tl_stflux(i,j,itemp)=tl_stflx(i,j,itemp)*rmask_wet(i,j)
1189# endif
1190# ifdef EMINUSP
1191 evap(i,j)=lheat(i,j)/hlv(i,j)
1192 tl_evap(i,j)=tl_lheat(i,j)/hlv(i,j)- &
1193 & tl_hlv(i,j)*evap(i,j)/hlv(i,j)
1194!^ stflux(i,j,isalt)=cff*(evap(i,j)-rain(i,j))
1195!^
1196 tl_stflx(i,j,isalt)=cff*tl_evap(i,j)
1197# ifdef MASKING
1198 evap(i,j)=evap(i,j)*rmask(i,j)
1199 tl_evap(i,j)=tl_evap(i,j)*rmask(i,j)
1200!^ stflux(i,j,isalt)=stflx(i,j,isalt)*rmask(i,j)
1201!^
1202 tl_stflux(i,j,isalt)=tl_stflx(i,j,isalt)*rmask(i,j)
1203# endif
1204# ifdef WET_DRY
1205 evap(i,j)=evap(i,j)*rmask_wet(i,j)
1206 tl_evap(i,j)=tl_evap(i,j)*rmask_wet(i,j)
1207!^ stflux(i,j,isalt)=stflx(i,j,isalt)*rmask_wet(i,j)
1208!^
1209 tl_stflux(i,j,isalt)=tl_stflx(i,j,isalt)*rmask_wet(i,j)
1210# endif
1211# endif
1212 END DO
1213 END DO
1214!
1215! Compute kinematic, surface wind stress (m2/s2).
1216!
1217 cff=0.5_r8/rho0
1218 DO j=jstrr,jendr
1219 DO i=istr,iendr
1220!^ sustr(i,j)=cff*(Taux(i-1,j)+Taux(i,j))
1221!^
1222 tl_sustr(i,j)=cff*(tl_taux(i-1,j)+tl_taux(i,j))
1223# ifdef MASKING
1224!^ sustr(i,j)=sustr(i,j)*umask(i,j)
1225!^
1226 tl_sustr(i,j)=tl_sustr(i,j)*umask(i,j)
1227# endif
1228# ifdef WET_DRY
1229!^ sustr(i,j)=sustr(i,j)*umask_wet(i,j)
1230!^
1231 tl_sustr(i,j)=tl_sustr(i,j)*umask_wet(i,j)
1232# endif
1233 END DO
1234 END DO
1235 DO j=jstr,jendr
1236 DO i=istrr,iendr
1237!^ svstr(i,j)=cff*(Tauy(i,j-1)+Tauy(i,j))
1238!^
1239 tl_svstr(i,j)=cff*(tl_tauy(i,j-1)+tl_tauy(i,j))
1240# ifdef MASKING
1241!^ svstr(i,j)=svstr(i,j)*vmask(i,j)
1242!^
1243 tl_svstr(i,j)=tl_svstr(i,j)*vmask(i,j)
1244# endif
1245# ifdef WET_DRY
1246!^ svstr(i,j)=svstr(i,j)*vmask_wet(i,j)
1247!^
1248 tl_svstr(i,j)=tl_svstr(i,j)*vmask_wet(i,j)
1249# endif
1250 END DO
1251 END DO
1252!
1253!-----------------------------------------------------------------------
1254! Exchange boundary data.
1255!-----------------------------------------------------------------------
1256!
1257 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1258!^ CALL exchange_r2d_tile (ng, tile, &
1259!^ & LBi, UBi, LBj, UBj, &
1260!^ & lrflx)
1261!^
1262 CALL exchange_r2d_tile (ng, tile, &
1263 & lbi, ubi, lbj, ubj, &
1264 & tl_lrflx)
1265!^ CALL exchange_r2d_tile (ng, tile, &
1266!^ & LBi, UBi, LBj, UBj, &
1267!^ & lhflx)
1268!^
1269 CALL exchange_r2d_tile (ng, tile, &
1270 & lbi, ubi, lbj, ubj, &
1271 & tl_lhflx)
1272!^ CALL exchange_r2d_tile (ng, tile, &
1273!^ & LBi, UBi, LBj, UBj, &
1274!^ & shflx)
1275!^
1276 CALL exchange_r2d_tile (ng, tile, &
1277 & lbi, ubi, lbj, ubj, &
1278 & tl_shflx)
1279!^ CALL exchange_r2d_tile (ng, tile, &
1280!^ & LBi, UBi, LBj, UBj, &
1281!^ & stflux(:,:,itemp))
1282!^
1283 CALL exchange_r2d_tile (ng, tile, &
1284 & lbi, ubi, lbj, ubj, &
1285 & tl_stflux(:,:,itemp))
1286# ifdef EMINUSP
1287!^ CALL exchange_r2d_tile (ng, tile, &
1288!^ & LBi, UBi, LBj, UBj, &
1289!^ & stflux(:,:,isalt))
1290!^
1291 CALL exchange_r2d_tile (ng, tile, &
1292 & lbi, ubi, lbj, ubj, &
1293 & tl_stflux(:,:,isalt))
1294!^ CALL exchange_r2d_tile (ng, tile, &
1295!^ & LBi, UBi, LBj, UBj, &
1296!^ & evap)
1297!^
1298 CALL exchange_r2d_tile (ng, tile, &
1299 & lbi, ubi, lbj, ubj, &
1300 & tl_evap)
1301# endif
1302!^ CALL exchange_u2d_tile (ng, tile, &
1303!^ & LBi, UBi, LBj, UBj, &
1304!^ & sustr)
1305!^
1306 CALL exchange_u2d_tile (ng, tile, &
1307 & lbi, ubi, lbj, ubj, &
1308 & tl_sustr)
1309!^ CALL exchange_v2d_tile (ng, tile, &
1310!^ & LBi, UBi, LBj, UBj, &
1311!^ & svstr)
1312!^
1313 CALL exchange_v2d_tile (ng, tile, &
1314 & lbi, ubi, lbj, ubj, &
1315 & tl_svstr)
1316 END IF
1317
1318# ifdef DISTRIBUTE
1319!^ CALL mp_exchange2d (ng, tile, iTLM, 4, &
1320!^ & LBi, UBi, LBj, UBj, &
1321!^ & NghostPoints, &
1322!^ & EWperiodic(ng), NSperiodic(ng), &
1323!^ & lrflx, lhflx, shflx, &
1324!^ & stflux(:,:,itemp))
1325!^
1326 CALL mp_exchange2d (ng, tile, itlm, 4, &
1327 & lbi, ubi, lbj, ubj, &
1328 & nghostpoints, &
1329 & ewperiodic(ng), nsperiodic(ng), &
1330 & tl_lrflx, tl_lhflx, tl_shflx, &
1331 & tl_stflux(:,:,itemp))
1332# ifdef EMINUSP
1333!^ CALL mp_exchange2d (ng, tile, iTLM, 2, &
1334!^ & LBi, UBi, LBj, UBj, &
1335!^ & NghostPoints, &
1336!^ & EWperiodic(ng), NSperiodic(ng), &
1337!^ & evap, &
1338!^ & stflux(:,:,isalt))
1339!^
1340 CALL mp_exchange2d (ng, tile, itlm, 2, &
1341 & lbi, ubi, lbj, ubj, &
1342 & nghostpoints, &
1343 & ewperiodic(ng), nsperiodic(ng), &
1344 & tl_evap, &
1345 & tl_stflux(:,:,isalt))
1346# endif
1347!^ CALL mp_exchange2d (ng, tile, iTLM, 2, &
1348!^ & LBi, UBi, LBj, UBj, &
1349!^ & NghostPoints, &
1350!^ & EWperiodic(ng), NSperiodic(ng), &
1351!^ & sustr, svstr)
1352!^
1353 CALL mp_exchange2d (ng, tile, itlm, 2, &
1354 & lbi, ubi, lbj, ubj, &
1355 & nghostpoints, &
1356 & ewperiodic(ng), nsperiodic(ng), &
1357 & tl_sustr, tl_svstr)
1358# endif
1359!
1360 RETURN
real(r8) function, public bulk_psiu(zol, pi)
Definition bulk_flux.F:1466
real(r8) function, public bulk_psit(zol, pi)
Definition bulk_flux.F:1531
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
integer, dimension(:), allocatable n
Definition mod_param.F:479
integer nghostpoints
Definition mod_param.F:710
real(dp) blk_dter
real(r8), dimension(:), allocatable blk_zt
real(r8), dimension(:), allocatable blk_zw
real(dp) vonkar
real(dp) blk_cpw
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
real(dp) cp
real(dp) blk_zabl
real(dp) blk_visw
real(dp) blk_tcw
real(dp) blk_rgas
integer isalt
integer itemp
real(dp) stefbo
real(dp) rhow
real(r8), dimension(:), allocatable blk_zq
real(dp) blk_beta
real(dp) g
real(dp) rho0
real(dp) blk_cpa
real(dp) emmiss
real(dp), parameter pi
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)

References mod_scalars::blk_beta, mod_scalars::blk_cpa, mod_scalars::blk_cpw, mod_scalars::blk_dter, mod_scalars::blk_rgas, mod_scalars::blk_tcw, mod_scalars::blk_visw, mod_scalars::blk_zabl, mod_scalars::blk_zq, mod_scalars::blk_zt, mod_scalars::blk_zw, bulk_flux_mod::bulk_psit(), bulk_flux_mod::bulk_psiu(), mod_scalars::cp, mod_scalars::emmiss, mod_scalars::ewperiodic, exchange_2d_mod::exchange_r2d_tile(), exchange_2d_mod::exchange_u2d_tile(), exchange_2d_mod::exchange_v2d_tile(), mod_scalars::g, mod_scalars::isalt, mod_scalars::itemp, mod_param::itlm, mp_exchange_mod::mp_exchange2d(), mod_param::nghostpoints, mod_scalars::nsperiodic, mod_scalars::pi, mod_scalars::rho0, mod_scalars::rhow, mod_scalars::stefbo, tl_bulk_psit(), tl_bulk_psiu(), and mod_scalars::vonkar.

Referenced by tl_bulk_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ tl_bulk_psit()

real(r8) function, public tl_bulk_flux_mod::tl_bulk_psit ( real(r8), intent(in) tl_zol,
real(r8), intent(in) zol,
real(dp), intent(in) pi )

Definition at line 1446 of file tl_bulk_flux.F.

1447!
1448!=======================================================================
1449! !
1450! This function evaluates the stability function for moisture and !
1451! heat by matching Kansas and free convection forms. The convective !
1452! form follows Fairall et al. (1996) with profile constants from !
1453! Grachev et al. (2000) BLM. The stable form is from Beljaars and !
1454! and Holtslag (1991). !
1455! !
1456!=======================================================================
1457!
1458 USE mod_kinds
1459!
1460! Function result
1461!
1462 real(r8) :: tl_bulk_psit
1463!
1464! Imported variable declarations.
1465!
1466 real(r8), intent(in) :: tl_ZoL, ZoL
1467 real(dp), intent(in) :: pi
1468!
1469! Local variable declarations.
1470!
1471 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
1472!
1473 real(r8) :: Fw, cff, psic, psik, x, y, fac
1474 real(r8) :: tl_Fw, tl_cff, tl_psic, tl_psik, tl_x, tl_y, tl_fac
1475!
1476!-----------------------------------------------------------------------
1477! Compute stability function, PSI.
1478!-----------------------------------------------------------------------
1479!
1480! Unstable conditions.
1481!
1482 IF (zol.lt.0.0_r8) THEN
1483 x=(1.0_r8-15.0_r8*zol)**0.5_r8
1484 IF (x.ne.0.0) THEN
1485 tl_x=-0.5_r8*15.0_r8*tl_zol/x
1486 ELSE
1487 tl_x=0.0_r8
1488 END IF
1489 psik=2.0_r8*log(0.5_r8*(1.0_r8+x))
1490 tl_psik=tl_x/(0.5_r8*(1.0_r8+x))
1491!
1492! For very unstable conditions, use free-convection (Fairall).
1493!
1494 cff=sqrt(3.0_r8)
1495 y=(1.0_r8-34.15_r8*zol)**r3
1496 tl_y=-r3*34.15_r8*tl_zol*(1.0_r8-34.15_r8*zol)**(r3-1.0_r8)
1497 fac=(1.0_r8+2.0_r8*y)/cff
1498 tl_fac=2.0_r8*tl_y/cff
1499 psic=1.5_r8*log(r3*(1.0_r8+y+y*y))- &
1500 & cff*atan(fac)+pi/cff
1501 tl_psic=tl_y*(1.0_r8+2.0_r8*y)*1.5_r8/(1.0_r8+y+y*y)- &
1502 & cff*tl_fac/(1.0_r8+fac*fac)
1503!! tl_psic=tl_y*(1.0_r8+2.0_r8*y)*1.5_r8/(1.0_r8+y+y*y)- &
1504!! & cff*tl_fac/SQRT(1.0_r8+fac*fac)
1505!
1506! Match Kansas and free-convection forms with weighting Fw.
1507!
1508 cff=zol*zol
1509 tl_cff=2.0_r8*tl_zol*zol
1510 fw=cff/(1.0_r8+cff)
1511 tl_fw=tl_cff/(1.0_r8+cff)-tl_cff*fw*fw/cff
1512 tl_bulk_psit=(1.0_r8-fw)*tl_psik-tl_fw*psik+ &
1513 & tl_fw*psic+fw*tl_psic
1514!
1515! Stable conditions.
1516!
1517 ELSE
1518 cff=min(50.0_r8,0.35_r8*zol)
1519 tl_cff=(0.5_r8-sign(0.5_r8,0.35_r8*zol-50.0_r8))*0.35_r8*tl_zol
1520 fac=exp(cff)
1521 tl_fac=tl_cff*fac
1522 tl_bulk_psit=-(3.0_r8*tl_zol*(1.0_r8+2.0_r8*zol)**0.5_r8+ &
1523 & 0.6667_r8*tl_zol/fac- &
1524 & tl_fac*0.6667_r8*(zol-14.28_r8)/(fac*fac))
1525 END IF
1526!
1527 RETURN

References tl_bulk_psit().

Referenced by rp_bulk_flux_mod::rp_bulk_flux_tile(), tl_bulk_flux_tile(), and tl_bulk_psit().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ tl_bulk_psiu()

real(r8) function, public tl_bulk_flux_mod::tl_bulk_psiu ( real(r8), intent(in) tl_zol,
real(r8), intent(in) zol,
real(dp), intent(in) pi )

Definition at line 1363 of file tl_bulk_flux.F.

1364!
1365!=======================================================================
1366! !
1367! This function evaluates the stability function for wind speed !
1368! by matching Kansas and free convection forms. The convective !
1369! form follows Fairall et al. (1996) with profile constants from !
1370! Grachev et al. (2000) BLM. The stable form is from Beljaars !
1371! and Holtslag (1991). !
1372! !
1373!=======================================================================
1374!
1375 USE mod_kinds
1376!
1377! Function result
1378!
1379 real(r8) :: tl_bulk_psiu
1380!
1381! Imported variable declarations.
1382!
1383 real(r8), intent(in) :: tl_ZoL, ZoL
1384 real(dp), intent(in) :: pi
1385!
1386! Local variable declarations.
1387!
1388 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
1389!
1390 real(r8) :: Fw, cff, psic, psik, x, y, fac
1391 real(r8) :: tl_Fw, tl_cff, tl_psic, tl_psik, tl_x, tl_y, tl_fac
1392!
1393!-----------------------------------------------------------------------
1394! Compute stability function, PSI.
1395!-----------------------------------------------------------------------
1396!
1397! Unstable conditions.
1398!
1399 IF (zol.lt.0.0_r8) THEN
1400 x=(1.0_r8-15.0_r8*zol)**0.25_r8
1401 tl_x=-0.25_r8*15.0_r8*tl_zol/((1.0_r8-15.0_r8*zol)**0.75_r8)
1402 psik=2.0_r8*log(0.5_r8*(1.0_r8+x))+ &
1403 & log(0.5_r8*(1.0_r8+x*x))- &
1404 & 2.0_r8*atan(x)+0.5_r8*pi
1405 tl_psik=tl_x/(0.5_r8*(1.0_r8+x))+ &
1406 & tl_x*x/(0.5_r8*(1.0_r8+x*x))- &
1407 & 2.0_r8*tl_x/(1.0_r8+x*x)
1408!! & 2.0_r8*tl_x/SQRT(1.0_r8+x*x)
1409!
1410! For very unstable conditions, use free-convection (Fairall).
1411!
1412 cff=sqrt(3.0_r8)
1413 y=(1.0_r8-10.15_r8*zol)**r3
1414 tl_y=-r3*10.15_r8*tl_zol*(1.0_r8-10.15_r8*zol)**(r3-1.0_r8)
1415 fac=(1.0_r8+2.0_r8*y)/cff
1416 tl_fac=2.0_r8*tl_y/cff
1417 psic=1.5_r8*log(r3*(1.0_r8+y+y*y))- &
1418 & cff*atan(fac)+pi/cff
1419 tl_psic=tl_y*(1.0_r8+2.0_r8*y)*1.5_r8/(1.0_r8+y+y*y)- &
1420 & cff*tl_fac/(1.0_r8+fac*fac)
1421!! & cff*tl_fac/SQRT(1.0_r8+fac*fac)
1422!
1423! Match Kansas and free-convection forms with weighting Fw.
1424!
1425 cff=zol*zol
1426 tl_cff=2.0_r8*tl_zol*zol
1427 fw=cff/(1.0_r8+cff)
1428 tl_fw=tl_cff/(1.0_r8+cff)-tl_cff*fw*fw/cff
1429 tl_bulk_psiu=(1.0_r8-fw)*tl_psik-tl_fw*psik+ &
1430 & tl_fw*psic+fw*tl_psic
1431!
1432! Stable conditions.
1433!
1434 ELSE
1435 cff=min(50.0_r8,0.35_r8*zol)
1436 tl_cff=(0.5_r8-sign(0.5_r8,0.35_r8*zol-50.0))*0.35_r8*tl_zol
1437 fac=exp(cff)
1438 tl_fac=tl_cff*fac
1439 tl_bulk_psiu=-(tl_zol+0.6667_r8*tl_zol/fac- &
1440 & tl_fac*0.6667_r8*(zol-14.28_r8)/(fac*fac))
1441 END IF
1442!
1443 RETURN

References tl_bulk_psiu().

Referenced by rp_bulk_flux_mod::rp_bulk_flux_tile(), tl_bulk_flux_tile(), and tl_bulk_psiu().

Here is the call graph for this function:
Here is the caller graph for this function: