ROMS
Loading...
Searching...
No Matches
tl_bulk_flux.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE tl_bulk_flux_mod
3#if (defined TANGENT || defined TL_IOMS) && defined BULK_FLUXES_NOT_YET
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine computes the bulk parameterization of surface wind !
13! stress and surface net heat fluxes. !
14! !
15! References: !
16! !
17! Fairall, C.W., E.F. Bradley, D.P. Rogers, J.B. Edson and G.S. !
18! Young, 1996: Bulk parameterization of air-sea fluxes for !
19! tropical ocean-global atmosphere Coupled-Ocean Atmosphere !
20! Response Experiment, JGR, 101, 3747-3764. !
21! !
22! Fairall, C.W., E.F. Bradley, J.S. Godfrey, G.A. Wick, J.B. !
23! Edson, and G.S. Young, 1996: Cool-skin and warm-layer !
24! effects on sea surface temperature, JGR, 101, 1295-1308. !
25! !
26! Liu, W.T., K.B. Katsaros, and J.A. Businger, 1979: Bulk !
27! parameterization of the air-sea exchange of heat and !
28! water vapor including the molecular constraints at !
29! the interface, J. Atmos. Sci, 36, 1722-1735. !
30! !
31! Adapted from COARE code written originally by David Rutgers and !
32! Frank Bradley. !
33! !
34! EMINUSP option for equivalent salt fluxes added by Paul Goodman !
35! (10/2004). !
36! !
37! Modified by Kate Hedstrom for COARE version 3.0 (03/2005). !
38! Modified by Jim Edson to correct specific hunidities. !
39! !
40! Reference: !
41! !
42! Fairall et al., 2003: J. Climate, 16, 571-591. !
43! !
44! HGA: This routine is deactivated and needs to be revised in the !
45! future if using it in the TLM. There is no need to recompute !
46! the NLM fluxes and pass them outside of this routine. The !
47! input fields are from data and not the governing equations. !
48! The computed surface net heat flux and surface freshwater flux !
49! are loaded into the data variable 'stflux' and used elsewhere !
50! when calculating the state variable stfflx. The other heat !
51! fluxes are for NLM output NetCDF files only. They are not !
52! needed in the TLM. Notice that the wind stress components !
53! are loaded to the state variables (7-Sep-2020). !
54! !
55!=======================================================================
56!
57 implicit none
58!
59 PRIVATE
60 PUBLIC :: tl_bulk_flux, tl_bulk_psiu, tl_bulk_psit
61!
62 CONTAINS
63!
64!***********************************************************************
65 SUBROUTINE tl_bulk_flux (ng, tile)
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
146 END SUBROUTINE tl_bulk_flux
147!
148!***********************************************************************
149 SUBROUTINE tl_bulk_flux_tile (ng, tile, &
150 & LBi, UBi, LBj, UBj, &
151 & IminS, ImaxS, JminS, JmaxS, &
152 & nrhs, &
153# ifdef MASKING
154 & rmask, umask, vmask, &
155# endif
156# ifdef WET_DRY
157 & rmask_wet, umask_wet, vmask_wet, &
158# endif
159 & alpha, tl_alpha, &
160 & beta, tl_beta, &
161 & rho, tl_rho, t, tl_t, &
162 & Hair, Pair, Tair, &
163 & Uwind, Vwind, &
164# ifdef CLOUDS
165 & cloud, &
166# endif
167# ifdef BBL_MODEL_NOT_YET
168 & Awave, Pwave, &
169# endif
170 & rain, &
171 & lhflx, tl_lhflx, &
172 & lrflx, tl_lrflx, &
173 & shflx, tl_shflx, &
174 & srflx, &
175 & stflux, tl_stflux, &
176# ifdef EMINUSP
177 & evap, tl_evap, &
178# endif
179 & sustr, tl_sustr, &
180 & svstr, tl_svstr)
181!***********************************************************************
182!
183 USE mod_param
184 USE mod_scalars
185!
186 USE bulk_flux_mod, ONLY : bulk_psiu, bulk_psit
187 USE exchange_2d_mod
188# ifdef DISTRIBUTE
189 USE mp_exchange_mod, ONLY : mp_exchange2d
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
1361 END SUBROUTINE tl_bulk_flux_tile
1362!
1363 FUNCTION tl_bulk_psiu (tl_ZoL, ZoL, pi)
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
1444 END FUNCTION tl_bulk_psiu
1445!
1446 FUNCTION tl_bulk_psit (tl_ZoL, ZoL, pi)
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
1528 END FUNCTION tl_bulk_psit
1529#endif
1530 END MODULE tl_bulk_flux_mod
bulk_flux_mod
Definition bulk_flux.F:2
bulk_flux_mod::bulk_psiu
real(r8) function, public bulk_psiu(zol, pi)
Definition bulk_flux.F:1466
bulk_flux_mod::bulk_psit
real(r8) function, public bulk_psit(zol, pi)
Definition bulk_flux.F:1531
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_r2d_tile
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:253
exchange_2d_mod::exchange_u2d_tile
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:440
exchange_2d_mod::exchange_v2d_tile
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:627
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_kinds
Definition mod_kinds.F:2
mod_kinds::r8
integer, parameter r8
Definition mod_kinds.F:28
mod_kinds::dp
integer, parameter dp
Definition mod_kinds.F:25
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_scalars
Definition mod_scalars.F:2
mod_scalars::blk_dter
real(dp) blk_dter
Definition mod_scalars.F:1501
mod_scalars::blk_zt
real(r8), dimension(:), allocatable blk_zt
Definition mod_scalars.F:1506
mod_scalars::blk_zw
real(r8), dimension(:), allocatable blk_zw
Definition mod_scalars.F:1507
mod_scalars::vonkar
real(dp) vonkar
Definition mod_scalars.F:469
mod_scalars::blk_cpw
real(dp) blk_cpw
Definition mod_scalars.F:1497
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::cp
real(dp) cp
Definition mod_scalars.F:456
mod_scalars::blk_zabl
real(dp) blk_zabl
Definition mod_scalars.F:1499
mod_scalars::blk_visw
real(dp) blk_visw
Definition mod_scalars.F:1503
mod_scalars::blk_tcw
real(dp) blk_tcw
Definition mod_scalars.F:1502
mod_scalars::blk_rgas
real(dp) blk_rgas
Definition mod_scalars.F:1498
mod_scalars::isalt
integer isalt
Definition mod_scalars.F:127
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::stefbo
real(dp) stefbo
Definition mod_scalars.F:460
mod_scalars::rhow
real(dp) rhow
Definition mod_scalars.F:462
mod_scalars::blk_zq
real(r8), dimension(:), allocatable blk_zq
Definition mod_scalars.F:1505
mod_scalars::blk_beta
real(dp) blk_beta
Definition mod_scalars.F:1500
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::blk_cpa
real(dp) blk_cpa
Definition mod_scalars.F:1496
mod_scalars::emmiss
real(dp) emmiss
Definition mod_scalars.F:461
mod_scalars::pi
real(dp), parameter pi
Definition mod_scalars.F:813
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
tl_bulk_flux_mod
Definition tl_bulk_flux.F:2
tl_bulk_flux_mod::tl_bulk_flux_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)
Definition tl_bulk_flux.F:181
tl_bulk_flux_mod::tl_bulk_psit
real(r8) function, public tl_bulk_psit(tl_zol, zol, pi)
Definition tl_bulk_flux.F:1447
tl_bulk_flux_mod::tl_bulk_psiu
real(r8) function, public tl_bulk_psiu(tl_zol, zol, pi)
Definition tl_bulk_flux.F:1364
tl_bulk_flux_mod::tl_bulk_flux
subroutine, public tl_bulk_flux(ng, tile)
Definition tl_bulk_flux.F:66
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3