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