ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
bulk_flux_mod Module Reference

Functions/Subroutines

subroutine, public bulk_flux (ng, tile)
 
subroutine bulk_flux_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, rmask, umask, vmask, rmask_wet, umask_wet, vmask_wet, alpha, beta, rho, t, u, v, hair, pair, tair, uwind, vwind, hwave, pwave_top, rain, lhflx, lrflx, shflx, srflx, evap, stflux, sustr, svstr)
 
real(r8) function, public bulk_psiu (zol, pi)
 
real(r8) function, public bulk_psit (zol, pi)
 

Function/Subroutine Documentation

◆ bulk_flux()

subroutine, public bulk_flux_mod::bulk_flux ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 106 of file bulk_flux.F.

107!***********************************************************************
108!
109 USE mod_stepping
110!
111! Imported variable declarations.
112!
113 integer, intent(in) :: ng, tile
114!
115! Local variable declarations.
116!
117 character (len=*), parameter :: MyFile = &
118 & __FILE__
119!
120# include "tile.h"
121!
122# ifdef PROFILE
123 CALL wclock_on (ng, inlm, 17, __line__, myfile)
124# endif
125 CALL bulk_flux_tile (ng, tile, &
126 & lbi, ubi, lbj, ubj, &
127 & imins, imaxs, jmins, jmaxs, &
128 & nrhs(ng), &
129# if defined ICE_MODEL && defined ICE_BULK_FLUXES
130 & liold(ng), linew(ng), &
131# endif
132# ifdef MASKING
133 & grid(ng) % rmask, &
134 & grid(ng) % umask, &
135 & grid(ng) % vmask, &
136# endif
137# ifdef WET_DRY
138 & grid(ng) % rmask_wet, &
139 & grid(ng) % umask_wet, &
140 & grid(ng) % vmask_wet, &
141# endif
142 & mixing(ng) % alpha, &
143 & mixing(ng) % beta, &
144 & ocean(ng) % rho, &
145 & ocean(ng) % t, &
146# ifdef WIND_MINUS_CURRENT
147 & ocean(ng) % u, &
148 & ocean(ng) % v, &
149# endif
150 & forces(ng) % Hair, &
151 & forces(ng) % Pair, &
152 & forces(ng) % Tair, &
153 & forces(ng) % Uwind, &
154 & forces(ng) % Vwind, &
155# ifdef CLOUDS
156 & forces(ng) % cloud, &
157# endif
158# if defined COARE_TAYLOR_YELLAND || defined DRENNAN
159 & forces(ng) % Hwave, &
160# endif
161# if defined COARE_TAYLOR_YELLAND || defined COARE_OOST || \
162 defined drennan
163 & forces(ng) % Pwave_top, &
164# endif
165# if !defined DEEPWATER_WAVES && \
166 (defined coare_taylor_yelland || defined coare_oost || \
167 defined drennan)
168 & forces(ng) % Lwavep, &
169# endif
170# if defined ICE_MODEL && defined ICE_BULK_FLUXES && \
171 defined ice_albedo && defined shortwave
172!! & FORCES(ng) % albedo, &
173 & forces(ng) % albedo_ice, &
174# endif
175 & forces(ng) % rain, &
176 & forces(ng) % lhflx, &
177 & forces(ng) % lrflx, &
178 & forces(ng) % shflx, &
179 & forces(ng) % srflx, &
180# if defined ICE_MODEL && defined ICE_BULK_FLUXES
181 & ice(ng) % Fi, &
182 & ice(ng) % Si, &
183 & forces(ng) % sustr_ai, &
184 & forces(ng) % svstr_ai, &
185 & forces(ng) % sustr_ao, &
186 & forces(ng) % svstr_ao, &
187 & forces(ng) % Qnet_ao, &
188# ifdef ICE_THERMO
189 & forces(ng) % Qnet_ai, &
190 & forces(ng) % srflx_ice, &
191 & forces(ng) % snow, &
192# endif
193# endif
194# ifdef EMINUSP
195 & forces(ng) % evap, &
196# endif
197 & forces(ng) % stflux, &
198 & forces(ng) % sustr, &
199 & forces(ng) % svstr)
200# ifdef PROFILE
201 CALL wclock_off (ng, inlm, 17, __line__, myfile)
202# endif
203!
204 RETURN
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
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

References bulk_flux_tile(), mod_forces::forces, mod_grid::grid, mod_ice::ice, mod_param::inlm, mod_mixing::mixing, mod_stepping::nrhs, mod_ocean::ocean, wclock_off(), and wclock_on().

Referenced by ad_main3d(), and main3d().

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

◆ bulk_flux_tile()

subroutine bulk_flux_mod::bulk_flux_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:,lbj:), intent(in) rmask,
real(r8), dimension(lbi:,lbj:), intent(in) umask,
real(r8), dimension(lbi:,lbj:), intent(in) vmask,
real(r8), dimension(lbi:,lbj:), intent(in) rmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) umask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) vmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) alpha,
real(r8), dimension(lbi:,lbj:), intent(in) beta,
real(r8), dimension(lbi:,lbj:,:), intent(in) rho,
real(r8), dimension(lbi:,lbj:,:,:,:), intent(in) t,
real(r8), dimension(lbi:,lbj:,:,:), intent(in) u,
real(r8), dimension(lbi:,lbj:,:,:), intent(in) v,
real(r8), dimension(lbi:,lbj:), intent(in) hair,
real(r8), dimension(lbi:,lbj:), intent(in) pair,
real(r8), dimension(lbi:,lbj:), intent(in) tair,
real(r8), dimension(lbi:,lbj:), intent(in) uwind,
real(r8), dimension(lbi:,lbj:), intent(in) vwind,
real(r8), dimension(lbi:,lbj:), intent(in) hwave,
real(r8), dimension(lbi:,lbj:), intent(in) pwave_top,
real(r8), dimension(lbi:,lbj:), intent(in) rain,
real(r8), dimension(lbi:,lbj:), intent(inout) lhflx,
real(r8), dimension(lbi:,lbj:), intent(inout) lrflx,
real(r8), dimension(lbi:,lbj:), intent(inout) shflx,
real(r8), dimension(lbi:,lbj:), intent(inout) srflx,
real(r8), dimension(lbi:,lbj:), intent(out) evap,
real(r8), dimension(lbi:,lbj:,:), intent(inout) stflux,
real(r8), dimension(lbi:,lbj:), intent(out) sustr,
real(r8), dimension(lbi:,lbj:), intent(out) svstr )
private

Definition at line 208 of file bulk_flux.F.

260!***********************************************************************
261!
262! Imported variable declarations.
263!
264 integer, intent(in) :: ng, tile
265 integer, intent(in) :: LBi, UBi, LBj, UBj
266 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
267 integer, intent(in) :: nrhs
268# if defined ICE_MODEL && defined ICE_BULK_FLUXES
269 integer, intent(in) :: liold, linew
270# endif
271!
272# ifdef ASSUMED_SHAPE
273# ifdef MASKING
274 real(r8), intent(in ) :: rmask(LBi:,LBj:)
275 real(r8), intent(in ) :: umask(LBi:,LBj:)
276 real(r8), intent(in ) :: vmask(LBi:,LBj:)
277# endif
278# ifdef WET_DRY
279 real(r8), intent(in ) :: rmask_wet(LBi:,LBj:)
280 real(r8), intent(in ) :: umask_wet(LBi:,LBj:)
281 real(r8), intent(in ) :: vmask_wet(LBi:,LBj:)
282# endif
283 real(r8), intent(in ) :: alpha(LBi:,LBj:)
284 real(r8), intent(in ) :: beta(LBi:,LBj:)
285 real(r8), intent(in ) :: rho(LBi:,LBj:,:)
286 real(r8), intent(in ) :: t(LBi:,LBj:,:,:,:)
287# ifdef WIND_MINUS_CURRENT
288 real(r8), intent(in ) :: u(LBi:,LBj:,:,:)
289 real(r8), intent(in ) :: v(LBi:,LBj:,:,:)
290# endif
291 real(r8), intent(in ) :: Hair(LBi:,LBj:)
292 real(r8), intent(in ) :: Pair(LBi:,LBj:)
293 real(r8), intent(in ) :: Tair(LBi:,LBj:)
294 real(r8), intent(in ) :: Uwind(LBi:,LBj:)
295 real(r8), intent(in ) :: Vwind(LBi:,LBj:)
296# ifdef CLOUDS
297 real(r8), intent(in ) :: cloud(LBi:,LBj:)
298# endif
299# if defined COARE_TAYLOR_YELLAND || defined DRENNAN
300 real(r8), intent(in ) :: Hwave(LBi:,LBj:)
301# endif
302# if defined COARE_TAYLOR_YELLAND || defined COARE_OOST || \
303 defined drennan
304 real(r8), intent(in ) :: Pwave_top(LBi:,LBj:)
305# endif
306# if !defined DEEPWATER_WAVES && \
307 (defined coare_taylor_yelland || defined coare_oost || \
308 defined drennan)
309 real(r8), intent(in ) :: Lwavep(LBi:,LBj:)
310# endif
311# if defined ICE_MODEL && defined ICE_BULK_FLUXES && \
312 defined ice_albedo && defined shortwave
313 real(r8), intent(in ) :: albedo_ice(LBi:,LBj:)
314# endif
315 real(r8), intent(in ) :: rain(LBi:,LBj:)
316 real(r8), intent(inout) :: lhflx(LBi:,LBj:)
317 real(r8), intent(inout) :: lrflx(LBi:,LBj:)
318 real(r8), intent(inout) :: shflx(LBi:,LBj:)
319 real(r8), intent(inout) :: srflx(LBi:,LBj:)
320 real(r8), intent(inout) :: stflux(LBi:,LBj:,:)
321# if defined ICE_MODEL && defined ICE_BULK_FLUXES
322 real(r8), intent(inout) :: Fi(LBi:,LBj:,:)
323 real(r8), intent(inout) :: Si(LBi:,LBj:,:,:)
324 real(r8), intent(out ) :: sustr_ai(LBi:,LBj:)
325 real(r8), intent(out ) :: svstr_ai(LBi:,LBj:)
326 real(r8), intent(out ) :: sustr_ao(LBi:,LBj:)
327 real(r8), intent(out ) :: svstr_ao(LBi:,LBj:)
328 real(r8), intent(out ) :: Qnet_ao(LBi:,LBj:)
329# ifdef ICE_THERMO
330 real(r8), intent(out ) :: Qnet_ai(LBi:,LBj:)
331 real(r8), intent(out ) :: srflx_ice(LBi:,LBj:)
332 real(r8), intent(out ) :: snow(LBi:,LBj:)
333# endif
334# endif
335# ifdef EMINUSP
336 real(r8), intent(out ) :: evap(LBi:,LBj:)
337# endif
338 real(r8), intent(out ) :: sustr(LBi:,LBj:)
339 real(r8), intent(out ) :: svstr(LBi:,LBj:)
340# else
341# ifdef MASKING
342 real(r8), intent(in ) :: rmask(LBi:UBi,LBj:UBj)
343 real(r8), intent(in ) :: umask(LBi:UBi,LBj:UBj)
344 real(r8), intent(in ) :: vmask(LBi:UBi,LBj:UBj)
345# endif
346# ifdef WET_DRY
347 real(r8), intent(in ) :: rmask_wet(LBi:UBi,LBj:UBj)
348 real(r8), intent(in ) :: umask_wet(LBi:UBi,LBj:UBj)
349 real(r8), intent(in ) :: vmask_wet(LBi:UBi,LBj:UBj)
350# endif
351 real(r8), intent(in ) :: alpha(LBi:UBi,LBj:UBj)
352 real(r8), intent(in ) :: beta(LBi:UBi,LBj:UBj)
353 real(r8), intent(in ) :: rho(LBi:UBi,LBj:UBj,N(ng))
354 real(r8), intent(in ) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
355# ifdef WIND_MINUS_CURRENT
356 real(r8), intent(in ) :: u(LBi:UBi,LBj:UBj,N(ng),3)
357 real(r8), intent(in ) :: v(LBi:UBi,LBj:UBj,N(ng),3)
358# endif
359 real(r8), intent(in ) :: Hair(LBi:UBi,LBj:UBj)
360 real(r8), intent(in ) :: Pair(LBi:UBi,LBj:UBj)
361 real(r8), intent(in ) :: Tair(LBi:UBi,LBj:UBj)
362 real(r8), intent(in ) :: Uwind(LBi:UBi,LBj:UBj)
363 real(r8), intent(in ) :: Vwind(LBi:UBi,LBj:UBj)
364# ifdef CLOUDS
365 real(r8), intent(in ) :: cloud(LBi:UBi,LBj:UBj)
366# endif
367# if defined COARE_TAYLOR_YELLAND || \
368 defined drennan
369 real(r8), intent(in ) :: Hwave(LBi:UBi,LBj:UBj)
370# endif
371# if defined COARE_TAYLOR_YELLAND || defined COARE_OOST || \
372 defined drennan
373 real(r8), intent(in ) :: Pwave_top(LBi:UBi,LBj:UBj)
374# endif
375# if !defined DEEPWATER_WAVES && \
376 (defined coare_taylor_yelland || defined coare_oost || \
377 defined drennan)
378 real(r8), intent(in ) :: Lwavep(LBi:UBi,LBj:UBj)
379# endif
380# if defined ICE_MODEL && defined ICE_BULK_FLUXES && \
381 defined ice_albedo && defined shortwave
382 real(r8), intent(in ) :: albedo_ice(LBi:UBi,LBj:UBj)
383# endif
384 real(r8), intent(in ) :: rain(LBi:UBi,LBj:UBj)
385 real(r8), intent(inout) :: lhflx(LBi:UBi,LBj:UBj)
386 real(r8), intent(inout) :: lrflx(LBi:UBi,LBj:UBj)
387 real(r8), intent(inout) :: shflx(LBi:UBi,LBj:UBj)
388 real(r8), intent(inout) :: srflx(LBi:UBi,LBj:UBj)
389 real(r8), intent(inout) :: stflux(LBi:UBi,LBj:UBj,NT(ng))
390# if defined ICE_MODEL && defined ICE_BULK_FLUXES
391 real(r8), intent(inout) :: Fi(LBi:UBi,LBj:UBj,nIceF)
392 real(r8), intent(inout) :: Si(LBi:UBi,LBj:UBj,2,nIceS)
393 real(r8), intent(out ) :: sustr_ai(LBi:UBi,LBj:UBj)
394 real(r8), intent(out ) :: svstr_ai(LBi:UBi,LBj:UBj)
395 real(r8), intent(out ) :: sustr_ao(LBi:UBi,LBj:UBj)
396 real(r8), intent(out ) :: svstr_ao(LBi:UBi,LBj:UBj)
397 real(r8), intent(out ) :: Qnet_ao(LBi:UBi,LBj:UBj)
398# ifdef ICE_THERMO
399 real(r8), intent(out ) :: Qnet_ai(LBi:UBi,LBj:UBj)
400 real(r8), intent(out ) :: srflx_ice(LBi:UBi,LBj:UBj)
401 real(r8), intent(out ) :: snow(LBi:UBi,LBj:UBj)
402# endif
403# endif
404# ifdef EMINUSP
405 real(r8), intent(out ) :: evap(LBi:UBi,LBj:UBj)
406# endif
407 real(r8), intent(out ) :: sustr(LBi:UBi,LBj:UBj)
408 real(r8), intent(out ) :: svstr(LBi:UBi,LBj:UBj)
409# endif
410!
411! Local variable declarations.
412!
413 integer :: Iter, i, j, k
414# if defined ICE_MODEL && defined ICE_BULK_FLUXES
415 integer :: li_stp
416# endif
417 integer, parameter :: IterMax = 3
418!
419 real(r8), parameter :: eps = 1.0e-20_r8
420 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
421!
422 real(r8) :: Bf, Cd, Hl, Hlw, Hscale, Hs, Hsr, IER
423 real(r8) :: PairM, RH, Taur
424 real(r8) :: Wspeed, ZQoL, ZToL
425# if defined ICE_MODEL && defined ICE_BULK_FLUXES
426 real(r8) :: Qsw_i, Qlw_i, Qlh_i, Qsh_i
427 real(r8) :: Qsati, TiceK
428 real(r8) :: Ce, Ch, dq_i, le_i, slp, vap_p_i
429# endif
430 real(r8) :: cff, cff1, cff2, cff3, cff4
431 real(r8) :: diffh, diffw, oL, upvel
432 real(r8) :: twopi_inv, wet_bulb
433# ifdef LONGWAVE
434 real(r8) :: e_sat, vap_p
435# endif
436# ifdef COOL_SKIN
437 real(r8) :: Clam, Fc, Hcool, Hsb, Hlb, Qbouy, Qcool, lambd
438# endif
439
440 real(r8), dimension(IminS:ImaxS) :: CC
441 real(r8), dimension(IminS:ImaxS) :: Cd10
442 real(r8), dimension(IminS:ImaxS) :: Ch10
443 real(r8), dimension(IminS:ImaxS) :: Ct10
444 real(r8), dimension(IminS:ImaxS) :: charn
445 real(r8), dimension(IminS:ImaxS) :: Ct
446 real(r8), dimension(IminS:ImaxS) :: Cwave
447 real(r8), dimension(IminS:ImaxS) :: Cwet
448 real(r8), dimension(IminS:ImaxS) :: delQ
449 real(r8), dimension(IminS:ImaxS) :: delQc
450 real(r8), dimension(IminS:ImaxS) :: delT
451 real(r8), dimension(IminS:ImaxS) :: delTc
452 real(r8), dimension(IminS:ImaxS) :: delW
453 real(r8), dimension(IminS:ImaxS) :: L
454 real(r8), dimension(IminS:ImaxS) :: L10
455 real(r8), dimension(IminS:ImaxS) :: Q
456 real(r8), dimension(IminS:ImaxS) :: Qair
457 real(r8), dimension(IminS:ImaxS) :: Qpsi
458 real(r8), dimension(IminS:ImaxS) :: Qsea
459 real(r8), dimension(IminS:ImaxS) :: Qstar
460 real(r8), dimension(IminS:ImaxS) :: rhoAir
461 real(r8), dimension(IminS:ImaxS) :: rhoSea
462 real(r8), dimension(IminS:ImaxS) :: Ri
463 real(r8), dimension(IminS:ImaxS) :: Ribcu
464 real(r8), dimension(IminS:ImaxS) :: Rr
465 real(r8), dimension(IminS:ImaxS) :: Scff
466 real(r8), dimension(IminS:ImaxS) :: TairC
467 real(r8), dimension(IminS:ImaxS) :: TairK
468 real(r8), dimension(IminS:ImaxS) :: Tcff
469 real(r8), dimension(IminS:ImaxS) :: Tpsi
470 real(r8), dimension(IminS:ImaxS) :: TseaC
471 real(r8), dimension(IminS:ImaxS) :: TseaK
472 real(r8), dimension(IminS:ImaxS) :: Tstar
473 real(r8), dimension(IminS:ImaxS) :: u10
474 real(r8), dimension(IminS:ImaxS) :: VisAir
475 real(r8), dimension(IminS:ImaxS) :: WaveLength
476 real(r8), dimension(IminS:ImaxS) :: Wgus
477 real(r8), dimension(IminS:ImaxS) :: Wmag
478 real(r8), dimension(IminS:ImaxS) :: Wpsi
479 real(r8), dimension(IminS:ImaxS) :: Wstar
480 real(r8), dimension(IminS:ImaxS) :: Zetu
481 real(r8), dimension(IminS:ImaxS) :: Zo10
482 real(r8), dimension(IminS:ImaxS) :: ZoT10
483 real(r8), dimension(IminS:ImaxS) :: ZoL
484 real(r8), dimension(IminS:ImaxS) :: ZoQ
485 real(r8), dimension(IminS:ImaxS) :: ZoT
486 real(r8), dimension(IminS:ImaxS) :: ZoW
487 real(r8), dimension(IminS:ImaxS) :: ZWoL
488
489 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hlv
490 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LHeat
491 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LRad
492 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: SHeat
493 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: SRad
494 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taux
495 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauy
496 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uair
497 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vair
498# ifdef ICE_THERMO
499 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Coef_Ice
500 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: RHS_Ice
501 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Qai
502# endif
503# if defined ICE_MODEL && defined ICE_BULK_FLUXES
504 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taux_Ice
505 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauy_Ice
506 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ice_thick
507 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: snow_thick
508# endif
509
510# include "set_bounds.h"
511!
512!=======================================================================
513! Atmosphere-Ocean bulk fluxes parameterization.
514!=======================================================================
515# ifdef WIND_MINUS_CURRENT
516!
517! Modify near-surface (2m or 10m) effective winds by subtracting the
518! ocean surface current (J Wilkin). See:
519!
520! Bye, J.A.T. and J.-O. Wolff, 1999: Atmosphere-ocean momentum exchange
521! in general circulation models. J. Phys. Oceanogr. 29, 671-691.
522!
523 DO j=jstr-1,jend+1
524 DO i=istr-1,min(iend+1,lm(ng))
525 uair(i,j)=uwind(i,j)- &
526 & 0.5_r8*(u(i,j,n(ng),nrhs)+u(i+1,j,n(ng),nrhs))
527 END DO
528 IF (domain(ng)%Eastern_Edge(tile)) THEN
529 uair(iend+1,j)=uwind(iend+1,j)-u(iend,j,n(ng),nrhs)
530 END IF
531 END DO
532 DO i=istr-1,iend+1
533 DO j=jstr-1,min(jend+1,mm(ng))
534 vair(i,j)=vwind(i,j)- &
535 & 0.5_r8*(v(i,j,n(ng),nrhs)+v(i,j+1,n(ng),nrhs))
536 END DO
537 IF (domain(ng)%Northern_Edge(tile)) THEN
538 vair(i,jend+1)=vwind(i,jend+1)-v(i,jend,n(ng),nrhs)
539 END IF
540 END DO
541# else
542!
543! Load wind components to local arrays.
544!
545 DO j=jstr-1,jend+1
546 DO i=istr-1,iend+1
547 uair(i,j)=uwind(i,j)
548 vair(i,j)=vwind(i,j)
549 END DO
550 END DO
551# endif
552!
553!-----------------------------------------------------------------------
554! Compute Atmosphere-ocean fluxes using a bulk flux parameterization.
555!-----------------------------------------------------------------------
556!
557 hscale=rho0*cp ! Celsius m/s to W/m2
558 twopi_inv=0.5_r8/pi
559# if defined ICE_MODEL && defined ICE_BULK_FLUXES
560 IF (perfectrst(ng) .and. iic(ng).eq.ntstart(ng)) THEN
561 li_stp=liold
562 ELSE
563 li_stp=linew
564 END IF
565# endif
566!
567 j_loop : DO j=jstr-1,jendr
568 DO i=istr-1,iendr
569!
570! Input bulk parameterization fields.
571!
572! (At input the flux data is converted to Celsius m/s. Here, we need to
573! convert back to W/m2 using 'Hscale').
574!
575 wmag(i)=sqrt(uair(i,j)*uair(i,j)+vair(i,j)*vair(i,j))
576 pairm=pair(i,j)
577 tairc(i)=tair(i,j) ! Celsius
578 tairk(i)=tairc(i)+273.16_r8 ! Kelvin
579 tseac(i)=t(i,j,n(ng),nrhs,itemp) ! Celsius
580 tseak(i)=tseac(i)+273.16_r8 ! Kelvin
581 rhosea(i)=rho(i,j,n(ng))+1000.0_r8
582 rh=hair(i,j)
583 srad(i,j)=srflx(i,j)*hscale
584# if defined ICE_MODEL && defined ICE_BULK_FLUXES && defined ICE_THERMO
585 srflx_ice(i,j)=srad(i,j)
586# endif
587 tcff(i)=alpha(i,j)
588 scff(i)=beta(i,j)
589!
590! Initialize.
591!
592 deltc(i)=0.0_r8
593 delqc(i)=0.0_r8
594 lheat(i,j)=lhflx(i,j)*hscale ! latent heat data, if any
595 sheat(i,j)=shflx(i,j)*hscale ! sensible heat data, if any
596 taur=0.0_r8
597 taux(i,j)=0.0_r8
598 tauy(i,j)=0.0_r8
599!
600!-----------------------------------------------------------------------
601! Compute net longwave radiation (W/m2), LRad.
602!-----------------------------------------------------------------------
603
604# if defined LONGWAVE
605!
606! Use Berliand (1952) formula to calculate net longwave radiation.
607! The equation for saturation vapor pressure is from Gill (Atmosphere-
608! Ocean Dynamics, pp 606). Here the coefficient in the cloud term
609! is assumed constant, but it is a function of latitude varying from
610! 1.0 at poles to 0.5 at the Equator).
611!
612 cff=(0.7859_r8+0.03477_r8*tairc(i))/ &
613 & (1.0_r8+0.00412_r8*tairc(i))
614 e_sat=10.0_r8**cff ! saturation vapor pressure (hPa or mbar)
615 vap_p=e_sat*rh ! water vapor pressure (hPa or mbar)
616 cff2=tairk(i)*tairk(i)*tairk(i)
617 cff1=cff2*tairk(i)
618 lrad(i,j)=-emmiss*stefbo* &
619 & (cff1*(0.39_r8-0.05_r8*sqrt(vap_p))* &
620 & (1.0_r8-0.6823_r8*cloud(i,j)*cloud(i,j))+ &
621 & cff2*4.0_r8*(tseak(i)-tairk(i)))
622
623# elif defined LONGWAVE_OUT
624!
625! Treat input longwave data as downwelling radiation only and add
626! outgoing IR from model sea surface temperature.
627!
628 lrad(i,j)=lrflx(i,j)*hscale- &
629 & emmiss*stefbo*tseak(i)*tseak(i)*tseak(i)*tseak(i)
630
631# else
632 lrad(i,j)=lrflx(i,j)*hscale
633# endif
634# ifdef MASKING
635 lrad(i,j)=lrad(i,j)*rmask(i,j)
636# endif
637# ifdef WET_DRY
638 lrad(i,j)=lrad(i,j)*rmask_wet(i,j)
639# endif
640!
641!-----------------------------------------------------------------------
642! Compute specific humidities (kg/kg).
643!
644! note that Qair(i) is the saturation specific humidity at Tair
645! Q(i) is the actual specific humidity
646! Qsea(i) is the saturation specific humidity at Tsea
647!
648! Saturation vapor pressure in mb is first computed and then
649! converted to specific humidity in kg/kg
650!
651! The saturation vapor pressure is computed from Teten formula
652! using the approach of Buck (1981):
653!
654! Esat(mb) = (1.0007_r8+3.46E-6_r8*PairM(mb))*6.1121_r8*
655! EXP(17.502_r8*TairC(C)/(240.97_r8+TairC(C)))
656!
657! The ambient vapor is found from the definition of the
658! Relative humidity:
659!
660! RH = W/Ws*100 ~ E/Esat*100 E = RH/100*Esat if RH is in %
661! E = RH*Esat if RH fractional
662!
663! The specific humidity is then found using the relationship:
664!
665! Q = 0.622 E/(P + (0.622-1)e)
666!
667! Q(kg/kg) = 0.62197_r8*(E(mb)/(PairM(mb)-0.378_r8*E(mb)))
668!
669!-----------------------------------------------------------------------
670!
671! Compute air saturation vapor pressure (mb), using Teten formula.
672!
673 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
674 & exp(17.502_r8*tairc(i)/(240.97_r8+tairc(i)))
675!
676! Compute specific humidity at Saturation, Qair (kg/kg).
677!
678 qair(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff+eps))
679!
680! Compute specific humidity, Q (kg/kg).
681!
682 IF (rh.lt.2.0_r8) THEN !RH fraction
683 cff=cff*rh !Vapor pres (mb)
684 q(i)=0.62197_r8* &
685 & (cff/(pairm-0.378_r8*cff+eps)) !Spec hum (kg/kg)
686 ELSE !RH input was actually specific humidity in g/kg
687 q(i)=rh/1000.0_r8 !Spec Hum (kg/kg)
688 END IF
689!
690! Compute water saturation vapor pressure (mb), using Teten formula.
691!
692 cff=(1.0007_r8+3.46e-6_r8*pairm)*6.1121_r8* &
693 & exp(17.502_r8*tseac(i)/(240.97_r8+tseac(i)))
694!
695! Vapor Pressure reduced for salinity (Kraus and Businger, 1994, pp42).
696!
697 cff=cff*0.98_r8
698!
699! Compute Qsea (kg/kg) from vapor pressure.
700!
701 qsea(i)=0.62197_r8*(cff/(pairm-0.378_r8*cff))
702!
703!-----------------------------------------------------------------------
704! Compute Monin-Obukhov similarity parameters for wind (Wstar),
705! heat (Tstar), and moisture (Qstar), Liu et al. (1979).
706!-----------------------------------------------------------------------
707!
708! Moist air density (kg/m3).
709!
710 rhoair(i)=pairm*100.0_r8/(blk_rgas*tairk(i)* &
711 & (1.0_r8+0.61_r8*q(i)))
712!
713! Kinematic viscosity of dry air (m2/s), Andreas (1989).
714!
715 visair(i)=1.326e-5_r8* &
716 & (1.0_r8+tairc(i)*(6.542e-3_r8+tairc(i)* &
717 & (8.301e-6_r8-4.84e-9_r8*tairc(i))))
718!
719! Compute latent heat of vaporization (J/kg) at sea surface, Hlv.
720!
721# ifdef REGCM_COUPLING
722 hlv(i,j)=2.5104e+6_r8 ! to be consistent with RegCM
723# else
724 hlv(i,j)=(2.501_r8-0.00237_r8*tseac(i))*1.0e+6_r8
725# endif
726!
727! Assume that wind is measured relative to sea surface and include
728! gustiness.
729!
730 wgus(i)=0.5_r8
731 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
732 delq(i)=qsea(i)-q(i)
733 delt(i)=tseac(i)-tairc(i)
734!
735! Neutral coefficients.
736!
737 zow(i)=0.0001_r8
738 u10(i)=delw(i)*log(10.0_r8/zow(i))/log(blk_zw(ng)/zow(i))
739 wstar(i)=0.035_r8*u10(i)
740 zo10(i)=0.011_r8*wstar(i)*wstar(i)/g+ &
741 & 0.11_r8*visair(i)/wstar(i)
742 cd10(i)=(vonkar/log(10.0_r8/zo10(i)))**2
743 ch10(i)=0.00115_r8
744 ct10(i)=ch10(i)/sqrt(cd10(i))
745 zot10(i)=10.0_r8/exp(vonkar/ct10(i))
746 cd=(vonkar/log(blk_zw(ng)/zo10(i)))**2
747!
748! Compute Richardson number.
749!
750 ct(i)=vonkar/log(blk_zt(ng)/zot10(i)) ! T transfer coefficient
751 cc(i)=vonkar*ct(i)/cd
752 deltc(i)=0.0_r8
753# ifdef COOL_SKIN
754 deltc(i)=blk_dter
755# endif
756 ribcu(i)=-blk_zw(ng)/(blk_zabl*0.004_r8*blk_beta**3)
757 ri(i)=-g*blk_zw(ng)*((delt(i)-deltc(i))+ &
758 & 0.61_r8*tairk(i)*delq(i))/ &
759 & (tairk(i)*delw(i)*delw(i)+eps)
760 IF (ri(i).lt.0.0_r8) THEN
761 zetu(i)=cc(i)*ri(i)/(1.0_r8+ri(i)/ribcu(i)) ! Unstable
762 ELSE
763 zetu(i)=cc(i)*ri(i)/(1.0_r8+3.0_r8*ri(i)/cc(i)) ! Stable
764 END IF
765 l10(i)=blk_zw(ng)/zetu(i)
766!
767! First guesses for Monon-Obukhov similarity scales.
768!
769 wstar(i)=delw(i)*vonkar/(log(blk_zw(ng)/zo10(i))- &
770 & bulk_psiu(blk_zw(ng)/l10(i),pi))
771 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
772 & (log(blk_zt(ng)/zot10(i))- &
773 & bulk_psit(blk_zt(ng)/l10(i),pi))
774 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
775 & (log(blk_zq(ng)/zot10(i))- &
776 & bulk_psit(blk_zq(ng)/l10(i),pi))
777
778# if defined COARE_30
779!
780! Momentum roughness length based on COARE 3.0 (Fairall et al 2003, JPO)
781!
782 IF (delw(i).gt.18.0_r8) THEN
783 charn(i)=0.018_r8
784 ELSE IF ((10.0_r8.lt.delw(i)).and.(delw(i).le.18.0_r8)) THEN
785 charn(i)=0.011_r8+ &
786 & 0.125_r8*(0.018_r8-0.011_r8)*(delw(i)-10.0_r8)
787 ELSE
788 charn(i)=0.011_r8
789 END IF
790
791# else
792!
793! Momentum roughness length based on COARE 3.5 (Edson et al 2013, JPO)
794!
795 charn(i)=min(0.028_r8,-0.005_r8+0.0017_r8*delw(i))
796
797# endif
798
799# if defined COARE_OOST || defined COARE_TAYLOR_YELLAND || \
800 defined drennan
801# if defined DEEPWATER_WAVES
802 cwave(i)=g*max(pwave_top(i,j),eps)*twopi_inv
803 wavelength(i)=cwave(i)*max(pwave_top(i,j),eps)
804# else
805 cwave(i)=lwavep(i,j)/max(pwave_top(i,j),eps)
806 wavelength(i)=lwavep(i,j)
807# endif
808# endif
809 END DO
810!
811! Iterate until convergence. It usually converges within 3 iterations.
812# if defined COARE_OOST || defined COARE_TAYLOR_YELLAND
813! Use wave info if we have it, two different options.
814# endif
815!
816 DO iter=1,itermax
817 DO i=istr-1,iendr
818# ifdef COARE_OOST
819 zow(i)=(25.0_r8/pi)*wavelength(i)* &
820 & (wstar(i)/cwave(i))**4.5_r8+ &
821 & 0.11_r8*visair(i)/(wstar(i)+eps)
822# elif defined COARE_TAYLOR_YELLAND
823 zow(i)=1200.0_r8*hwave(i,j)* &
824 & (hwave(i,j)/wavelength(i))**4.5_r8+ &
825 & 0.11_r8*visair(i)/(wstar(i)+eps)
826# elif defined DRENNAN
827 zow(i)=(3.35_r8)*hwave(i,j)* &
828 & min((wstar(i)/cwave(i)),0.1_r8)**3.4_r8+ &
829 & 0.11_r8*visair(i)/(wstar(i)+eps)
830# else
831 zow(i)=charn(i)*wstar(i)*wstar(i)/g+ &
832 & 0.11_r8*visair(i)/(wstar(i)+eps)
833# endif
834# if defined DRAGLIM_DAVIS
835 zow(i)=max(zow(i),1.27e-7_r8) ! Davis et al. (2008)
836 zow(i)=min(zow(i),2.85e-3_r8)
837# endif
838 rr(i)=zow(i)*wstar(i)/visair(i)
839!
840#if defined COARE_30
841!
842! Moisture and thermal roughness lengths, COARE 3.0
843!
844 zoq(i)=min(1.15e-4_r8,5.5e-5_r8/rr(i)**0.6_r8)
845#else
846!
847! In COARE 3.5, moisture and thermal roughness lengths are adjusted
848! to reflect COARE 3.0 fluxes (J. Edson, pers. comm, 11/2023)
849!
850 zoq(i)=min(1.6e-4_r8,5.8e-5_r8/(rr(i)**0.72_r8))
851# endif
852 zot(i)=zoq(i)
853!
854! Compute Monin-Obukhov stability parameter, Z/L.
855!
856 zol(i)=vonkar*g*blk_zw(ng)* &
857 & (tstar(i)*(1.0_r8+0.61_r8*q(i))+ &
858 & 0.61_r8*tairk(i)*qstar(i))/ &
859 & (tairk(i)*wstar(i)*wstar(i)* &
860 & (1.0_r8+0.61_r8*q(i))+eps)
861 l(i)=blk_zw(ng)/(zol(i)+eps)
862!
863! Evaluate stability functions at Z/L.
864!
865 wpsi(i)=bulk_psiu(zol(i),pi)
866 tpsi(i)=bulk_psit(blk_zt(ng)/l(i),pi)
867 qpsi(i)=bulk_psit(blk_zq(ng)/l(i),pi)
868# ifdef COOL_SKIN
869 cwet(i)=0.622_r8*hlv(i,j)*qsea(i)/ &
870 & (blk_rgas*tseak(i)*tseak(i))
871 delqc(i)=cwet(i)*deltc(i)
872# endif
873!
874! Compute wind scaling parameters, Wstar.
875!
876 wstar(i)=max(eps,delw(i)*vonkar/ &
877 & (log(blk_zw(ng)/zow(i))-wpsi(i)))
878 tstar(i)=-(delt(i)-deltc(i))*vonkar/ &
879 & (log(blk_zt(ng)/zot(i))-tpsi(i))
880 qstar(i)=-(delq(i)-delqc(i))*vonkar/ &
881 & (log(blk_zq(ng)/zoq(i))-qpsi(i))
882!
883! Compute gustiness in wind speed.
884!
885 bf=-g/tairk(i)* &
886 & wstar(i)*(tstar(i)+0.61_r8*tairk(i)*qstar(i))
887 IF (bf.gt.0.0_r8) THEN
888 wgus(i)=blk_beta*(bf*blk_zabl)**r3
889 ELSE
890 wgus(i)=0.2_r8
891 END IF
892 delw(i)=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
893
894# ifdef COOL_SKIN
895!
896!-----------------------------------------------------------------------
897! Cool Skin correction.
898!-----------------------------------------------------------------------
899!
900! Cool skin correction constants. Clam: part of Saunders constant
901! lambda; Cwet: slope of saturation vapor.
902!
903 clam=16.0_r8*g*blk_cpw*(rhosea(i)*blk_visw)**3.0_r8/ &
904 & (blk_tcw*blk_tcw*rhoair(i)*rhoair(i))
905!
906! Set initial guesses for cool-skin layer thickness (Hcool).
907!
908 hcool=0.001_r8
909!
910! Backgound sensible and latent heat.
911!
912 hsb=-rhoair(i)*blk_cpa*wstar(i)*tstar(i)
913 hlb=-rhoair(i)*hlv(i,j)*wstar(i)*qstar(i)
914!
915! Mean absoption in cool-skin layer.
916!
917 fc=0.065_r8+11.0_r8*hcool- &
918 & (1.0_r8-exp(-hcool*1250.0_r8))*6.6e-5_r8/hcool
919!
920! Total cooling at the interface.
921!
922 qcool=lrad(i,j)+hsb+hlb-srad(i,j)*fc
923 qbouy=tcff(i)*qcool+scff(i)*hlb*blk_cpw/hlv(i,j)
924!
925! Compute temperature and moisture change.
926!
927 IF ((qcool.gt.0.0_r8).and.(qbouy.gt.0.0_r8)) THEN
928 lambd=6.0_r8/ &
929 & (1.0_r8+ &
930 & (clam*qbouy/(wstar(i)+eps)**4.0_r8)**0.75_r8)**r3
931 hcool=lambd*blk_visw/(sqrt(rhoair(i)/rhosea(i))* &
932 & wstar(i)+eps)
933 deltc(i)=qcool*hcool/blk_tcw
934 ELSE
935 deltc(i)=0.0_r8
936 END IF
937 delqc(i)=cwet(i)*deltc(i)
938# endif
939 END DO
940 END DO
941!
942!-----------------------------------------------------------------------
943! Compute Atmosphere/Ocean fluxes.
944!-----------------------------------------------------------------------
945!
946 DO i=istr-1,iendr
947!
948! Compute nondimensional transfer coefficients for stress (Cd),
949! sensible heat (Ch), and latent heat (Ce).
950!
951 wspeed=sqrt(wmag(i)*wmag(i)+wgus(i)*wgus(i))
952! Cd=Wstar(i)*Wstar(i)/(Wspeed*Wspeed+eps) ! used prior to 02/2024
953# if defined ICE_MODEL && defined ICE_BULK_FLUXES && defined ICE_THERMO
954 ch=wstar(i)*tstar(i)/(-wspeed*delt(i)+0.0098_r8*blk_zt(ng))
955 ce=wstar(i)*qstar(i)/(-wspeed*delq(i)+eps)
956# endif
957!
958! Compute turbulent sensible heat flux (W/m2), Hs.
959!
960 hs=-blk_cpa*rhoair(i)*wstar(i)*tstar(i)
961!
962! Compute sensible heat flux (W/m2) due to rainfall (kg/m2/s), Hsr.
963!
964 diffw=2.11e-5_r8*(tairk(i)/273.16_r8)**1.94_r8
965 diffh=0.02411_r8*(1.0_r8+tairc(i)* &
966 & (3.309e-3_r8-1.44e-6_r8*tairc(i)))/ &
967 & (rhoair(i)*blk_cpa+eps)
968 cff=qair(i)*hlv(i,j)/(blk_rgas*tairk(i)*tairk(i))
969 wet_bulb=1.0_r8/(1.0_r8+0.622_r8*(cff*hlv(i,j)*diffw)/ &
970 & (blk_cpa*diffh))
971 hsr=abs(rain(i,j))*wet_bulb*blk_cpw* &
972 & ((tseac(i)-tairc(i))+(qsea(i)-q(i))*hlv(i,j)/blk_cpa)
973# ifndef REGCM_COUPLING
974 sheat(i,j)=(hs+hsr)
975# endif
976# ifdef MASKING
977 sheat(i,j)=sheat(i,j)*rmask(i,j)
978# endif
979# ifdef WET_DRY
980 sheat(i,j)=sheat(i,j)*rmask_wet(i,j)
981# endif
982!
983! Compute turbulent latent heat flux (W/m2), Hl.
984!
985 hl=-hlv(i,j)*rhoair(i)*wstar(i)*qstar(i)
986
987# if defined ICE_MODEL && defined ICE_BULK_FLUXES && defined ICE_THERMO
988!
989! Subtract loss of heat associated with latent heat flux from ocean to
990! melting snow (SMD).
991!
992 IF (rain(i,j).lt.0.0_r8) THEN
993 hl=hl+3.347e+5_r8*rain(i,j)
994 END IF
995# endif
996!
997! Compute Webb correction (Webb effect) to latent heat flux, Hlw.
998!
999 upvel=-1.61_r8*wstar(i)*qstar(i)- &
1000 & (1.0_r8+1.61_r8*q(i))*wstar(i)*tstar(i)/tairk(i)
1001 hlw=rhoair(i)*hlv(i,j)*upvel*q(i)
1002# ifndef REGCM_COUPLING
1003 lheat(i,j)=(hl+hlw)
1004# endif
1005# ifdef MASKING
1006 lheat(i,j)=lheat(i,j)*rmask(i,j)
1007# endif
1008# ifdef WET_DRY
1009 lheat(i,j)=lheat(i,j)*rmask_wet(i,j)
1010# endif
1011!
1012! Horizontal momentum flux (N/m2) due to rain (kg/m2/s) impact
1013! traveling at 85% of air velocity
1014!
1015 taur=0.85_r8*abs(rain(i,j))*wmag(i)
1016!
1017! Sum of stresses (N/m2) expressed as friction velocities.
1018! Normalization by Wmag (m/s) so that stress components (Taux,
1019! Tauy in N/m2) are cff times the respective wind components
1020! (Uair,Vair in m/s).
1021!
1022 cff=rhoair(i)*(wstar(i)*wstar(i)+taur/rhoair(i))/(wmag(i)+eps)
1023!
1024! In prior versions incorrect normalization included gustiness in
1025! Wspeed causing the vector stress magnitude to be too be low by ~3%
1026! (Corrected J. Wilkin and F. Parela-Roman 02/2024)
1027!
1028 taux(i,j)=cff*uair(i,j)
1029# ifdef MASKING
1030 taux(i,j)=taux(i,j)*rmask(i,j)
1031# endif
1032# ifdef WET_DRY
1033 taux(i,j)=taux(i,j)*rmask_wet(i,j)
1034# endif
1035 tauy(i,j)=cff*vair(i,j)
1036# ifdef MASKING
1037 tauy(i,j)=tauy(i,j)*rmask(i,j)
1038# endif
1039# ifdef WET_DRY
1040 tauy(i,j)=tauy(i,j)*rmask_wet(i,j)
1041# endif
1042# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1043 taux_ice(i,j)=rhoair(i)*cd_ai(ng)*uwind(i,j)*wspeed
1044 tauy_ice(i,j)=rhoair(i)*cd_ai(ng)*vwind(i,j)*wspeed
1045
1046# ifdef ICE_THERMO
1047!
1048! Correct ocean net short-wave radiation for ice cover and convert to
1049! kinematic units.
1050!
1051 coef_ice(i,j)=0.0_r8
1052 rhs_ice(i,j)=0.0_r8
1053!
1054! Sensible heat flux over ice.
1055!
1056 IF (si(i,j,li_stp,isaice).gt.min_ai(ng)) THEN
1057 ticek=fi(i,j,icisst)+273.16_r8
1058 ELSE
1059 ticek=t(i,j,n(ng),nrhs,itemp)+273.16_r8
1060 ENDIF
1061 cff=rhoair(i)*blk_cpa*ch*delw(i) ! (W/m2)/K
1062 qsh_i=cff*(tairk(i)+0.0098_r8*blk_zt(ng)-ticek) ! W/m2
1063 coef_ice(i,j)=coef_ice(i,j)+cff ! (W/m2)/K
1064 rhs_ice(i,j)=rhs_ice(i,j)+ &
1065 & cff*(tairk(i)+0.0098_r8*blk_zt(ng)) ! W/m2
1066!
1067! Add in the sensible heat transfer associated with warm rain or cold
1068! snow contacting the ice usinf the wet bulb factor in the sensible
1069! heat flux for water calculation above. This, of course, is a bit off
1070! if the actual surface is snow or meltwater covered. But in those
1071! cases there are inaccuracies in other terms as well (SMD).
1072!
1073 cff=wet_bulb*2093.0_r8
1074 qsh_i=qsh_i+ &
1075 & abs(rain(i,j))*cff* &
1076 & (ticek-tairc(i)+(qsea(i)-q(i))*hlv(i,j)/blk_cpa)
1077 coef_ice(i,j)=coef_ice(i,j) - &
1078 & abs(rain(i,j))*cff
1079!
1080! Use TairK here, unlike in above calculation of Hsr, since
1081! Fi(:,:,icIsst) will be calculated in Kelvin (SMD, 10.20.20).
1082!
1083 rhs_ice(i,j)=rhs_ice(i,j)+ &
1084 & abs(rain(i,j))*cff* &
1085 & (-tairk(i)+(qsea(i)-q(i))*hlv(i,j)/blk_cpa)
1086!
1087! Latent heat flux over ice.
1088!
1089 le_i=2.834e+6_r8
1090!
1091! Compute saturation specific humidity over the ice (Qsati) from
1092! Parkinson and Washington (1979).
1093!
1094 slp=pair(i,j)*100.0_r8 ! Convert back to Pa from mb
1095 cff=611.0_r8* &
1096 & 10.0_r8**(9.5_r8*(ticek-273.16_r8)/(ticek-7.66_r8))
1097!! Qsati=eps*cff/(slp-(1.0_r8-.62197_r8)*cff)
1098 qsati=0.62197_r8*cff/(slp-(1.0_r8-0.62197_r8)*cff)
1099 dq_i=q(i)-qsati
1100 qlh_i=rhoair(i)*ce*delw(i)* &
1101 & le_i*dq_i
1102!
1103! Add in the Latent heat gain from rain freezing on ice/snow surface.
1104! If meltwater is already present then assume rain just pools,
1105! otherwise rain may freeze or ice may melt (SMD).
1106!
1107 IF ((rain(i,j).gt.0.0_r8).and. &
1108 & (si(i,j,li_stp,ishmel).eq.0.0_r8)) THEN
1109 qlh_i=qlh_i+3.347e+5_r8*rain(i,j)
1110 END IF
1111 rhs_ice(i,j)=rhs_ice(i,j)+qlh_i
1112!
1113! Short-wave radiation (W/m**2) to ice.
1114!
1115 qsw_i=(1.0_r8-albedo_ice(i,j))*srflx_ice(i,j)
1116!
1117! A question remains as to what fraction of the solar radiation should
1118! be considered to contribute to the surface temperature of the
1119! ice/snow/meltwater? Leaving unchanged for now, assuming entire
1120! 'Qsw_i' contributes to ice surface temperature although it should
1121! probably be something less (SMD).
1122!
1123 rhs_ice(i,j)=rhs_ice(i,j)+qsw_i
1124!
1125! Net longwave radiation from ice.
1126!
1127 cff=q(i)/(1.0_r8+q(i))
1128 vap_p_i=slp*cff/(0.62197_r8+cff)
1129
1130# ifdef LONGWAVE
1131 cff1=0.39_r8-0.005_r8*sqrt(vap_p_i)
1132 cff2=1.0_r8-0.65_r8*cloud(i,j)
1133 cff3=stefbo*emmiss*(tairk(i)*tairk(i)*tairk(i))
1134 cff4=cff3*tairk(i)
1135 qlw_i=cff4*cff1*cff2+ &
1136 & 4.0_r8*cff3*(ticek-tairk(i))
1137 coef_ice(i,j)=coef_ice(i,j)+ &
1138 & 4.0_r8*cff3
1139 rhs_ice(i,j)=rhs_ice(i,j)- &
1140 & cff4*(cff1*cff2-4.0_r8)
1141
1142# elif defined LONGWAVE_OUT
1143!
1144! Option where the incoming longwave radiation (downward) is given as
1145! input (lrflx). The following should be consistent, with the implicit
1146! (pseudo) approach to solving for 'Fi(:,:,icIsst). Here, 'Qlw_i' is
1147! the flux from ice to atmosphere, thus the negative sign (SMD).
1148!
1149 cff3=stefbo*emmiss*(ticek**3)
1150 cff4=cff3*ticek
1151 qlw_i=-(lrflx(i,j)*hscale-cff4)
1152 coef_ice(i,j)=coef_ice(i,j)+ &
1153 & 4.0_r8*cff3
1154 rhs_ice(i,j)=rhs_ice(i,j)+ &
1155 & lrflx(i,j)*hscale+3.0_r8*cff4
1156
1157# else
1158 cff3=stefbo*emmiss*(tairk(i)**3)
1159 qlw_i=lrflx(i,j)*hscale
1160 coef_ice(i,j)=coef_ice(i,j)+ &
1161 & 4.0_r8*cff3
1162 rhs_ice(i,j)=rhs_ice(i,j)- &
1163 & lrflx(i,j)*hscale- &
1164 & 4.0_r8*cff3*(ticek-2.0_r8*tairk(i))
1165# endif
1166!
1167! Net heat flux from ice to atmosphere. The heat equation terms are
1168! computed with quantities in Kelvin to facilitate the computation of
1169! seaice surface temperature, Fi(:,:,icIsst), in Kelivin in module
1170! 'ice_mk.h' and then converted to Celsius afterwards (SMD).
1171!
1172 qai(i,j)=-qsh_i-qlh_i-qsw_i+qlw_i
1173# endif
1174# endif
1175 END DO
1176 END DO j_loop
1177!
1178!=======================================================================
1179! Compute surface net heat flux and surface wind stress.
1180!=======================================================================
1181!
1182! Compute kinematic, surface, net heat flux (degC m/s). Notice that
1183! the signs of latent and sensible fluxes are reversed because fluxes
1184! calculated from the bulk formulations above are positive out of the
1185! ocean.
1186!
1187! For EMINUSP option, EVAP = LHeat (W/m2) / Hlv (J/kg) = kg/m2/s
1188! PREC = rain = kg/m2/s
1189!
1190! To convert these rates to m/s divide by freshwater density, rhow.
1191!
1192! Note that when the air is undersaturated in water vapor (Q < Qsea)
1193! the model will evaporate and LHeat > 0:
1194!
1195! LHeat positive out of the ocean
1196! evap positive out of the ocean
1197!
1198! Note that if evaporating, the salt flux is positive
1199! and if raining, the salt flux is negative
1200!
1201! Note that stflux(:,:,isalt) is the E-P flux. The fresh water flux
1202! is positive out of the ocean and the salt flux is positive into the
1203! ocean. It is multiplied by surface salinity when computing state
1204! variable stflx(:,:,isalt) in "set_vbc.F".
1205!
1206 hscale=1.0_r8/(rho0*cp) ! W/2 to Celsius m/s
1207 cff=1.0_r8/rhow
1208 DO j=jstrr,jendr
1209 DO i=istrr,iendr
1210# if defined ALBEDO && defined SHORTWAVE
1211!! srflx(i,j) = srflx(i,j)*(1.0_r8-albedo(i,j))
1212# endif
1213# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1214!
1215! Limit shortwave radiation to be used in computing penetrative
1216! radiation by presence of ice alone.
1217!
1218 cff1=1.0_r8/(si(i,j,li_stp,isaice)+eps)
1219 ice_thick(i,j) =si(i,j,li_stp,ishice)*cff1
1220 snow_thick(i,j)=si(i,j,li_stp,ishsno)*cff1
1221
1222 IF (ice_thick(i,j).le.0.1_r8) THEN
1223 cff2=ice_thick(i,j)*5.0_r8
1224 ELSE
1225 cff2=0.1_r8*(5.0_r8+(ice_thick(i,j)-0.1_r8))
1226 ENDIF
1227
1228 cff2=cff2+snow_thick(i,j)*20.0_r8
1229 srflx(i,j)=(1.0_r8-si(i,j,li_stp,isaice))*srflx(i,j)+ &
1230 & si(i,j,li_stp,isaice)*srflx(i,j)*exp(-cff2)
1231!
1232! Net precipitation (m/s). Negative values of 'rain' indicate frozen
1233! precipitation (SMD).
1234!
1235 IF (rain(i,j).ge.0.0_r8) THEN
1236 snow(i,j)=0.0_r8
1237 ELSE
1238 snow(i,j)=abs(rain(i,j))/snowdryrho(ng)
1239 END IF
1240# endif
1241!
1242! Load surface fluxes to state arrays.
1243!
1244# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1245 fi(i,j,icqcon)=coef_ice(i,j) !(W/m2)/K
1246 fi(i,j,icqrhs)=rhs_ice(i,j) ! W/m2
1247# endif
1248 lrflx(i,j)=lrad(i,j)*hscale ! C m/s
1249 lhflx(i,j)=-lheat(i,j)*hscale ! C m/s
1250 shflx(i,j)=-sheat(i,j)*hscale ! C m/s
1251 stflux(i,j,itemp)=(srflx(i,j)+lrflx(i,j)+ &
1252 & lhflx(i,j)+shflx(i,j)) ! C m/s
1253
1254# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1255 qnet_ao(i,j)=-stflux(i,j,itemp)/hscale ! W/m2
1256 qnet_ai(i,j)=qai(i,j) ! W/m2
1257# endif
1258# ifdef MASKING
1259 stflux(i,j,itemp)=stflux(i,j,itemp)*rmask(i,j)
1260# endif
1261# ifdef WET_DRY
1262 stflux(i,j,itemp)=stflux(i,j,itemp)*rmask_wet(i,j)
1263# endif
1264# ifdef EMINUSP
1265 evap(i,j)=lheat(i,j)/(hlv(i,j)+eps) ! kg/m2/s
1266# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1267 evap(i,j)=(1.0_r8-si(i,j,li_stp,isaice))*evap(i,j)
1268# endif
1269# ifdef MASKING
1270 evap(i,j)=evap(i,j)*rmask(i,j)
1271# endif
1272# ifdef WET_DRY
1273 evap(i,j)=evap(i,j)*rmask_wet(i,j)
1274# endif
1275# ifdef SALINITY
1276 stflux(i,j,isalt)=cff*(evap(i,j)-rain(i,j)) ! m/s
1277# ifdef MASKING
1278 stflux(i,j,isalt)=stflux(i,j,isalt)*rmask(i,j)
1279# endif
1280# ifdef WET_DRY
1281 stflux(i,j,isalt)=stflux(i,j,isalt)*rmask_wet(i,j)
1282# endif
1283# endif
1284# endif
1285 END DO
1286 END DO
1287!
1288! Compute kinematic, surface wind stress (m2/s2).
1289!
1290 cff=0.5_r8/rho0
1291 DO j=jstrr,jendr
1292 DO i=istr,iendr
1293 sustr(i,j)=cff*(taux(i-1,j)+taux(i,j))
1294# ifdef MASKING
1295 sustr(i,j)=sustr(i,j)*umask(i,j)
1296# endif
1297# ifdef WET_DRY
1298 sustr(i,j)=sustr(i,j)*umask_wet(i,j)
1299# endif
1300# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1301 sustr_ao(i,j)=sustr(i,j)
1302 sustr_ai(i,j)=taux_ice(i,j)
1303# endif
1304 END DO
1305 END DO
1306 DO j=jstr,jendr
1307 DO i=istrr,iendr
1308 svstr(i,j)=cff*(tauy(i,j-1)+tauy(i,j))
1309# ifdef MASKING
1310 svstr(i,j)=svstr(i,j)*vmask(i,j)
1311# endif
1312# ifdef WET_DRY
1313 svstr(i,j)=svstr(i,j)*vmask_wet(i,j)
1314# endif
1315# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1316 svstr_ao(i,j)=svstr(i,j)
1317 svstr_ai(i,j)=tauy_ice(i,j)
1318# endif
1319 END DO
1320 END DO
1321!
1322!-----------------------------------------------------------------------
1323! Exchange boundary data.
1324!-----------------------------------------------------------------------
1325!
1326 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1327 CALL exchange_r2d_tile (ng, tile, &
1328 & lbi, ubi, lbj, ubj, &
1329 & lrflx)
1330
1331 CALL exchange_r2d_tile (ng, tile, &
1332 & lbi, ubi, lbj, ubj, &
1333 & lhflx)
1334
1335 CALL exchange_r2d_tile (ng, tile, &
1336 & lbi, ubi, lbj, ubj, &
1337 & shflx)
1338
1339 CALL exchange_r2d_tile (ng, tile, &
1340 & lbi, ubi, lbj, ubj, &
1341 & stflux(:,:,itemp))
1342
1343# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1344 CALL exchange_r2d_tile (ng, tile, &
1345 & lbi, ubi, lbj, ubj, &
1346 & srflx)
1347
1348 CALL exchange_r2d_tile (ng, tile, &
1349 & lbi, ubi, lbj, ubj, &
1350 & qnet_ao)
1351
1352# ifdef ICE_THERMO
1353 CALL exchange_r2d_tile (ng, tile, &
1354 & lbi, ubi, lbj, ubj, &
1355 & qnet_ai)
1356
1357 CALL exchange_r2d_tile (ng, tile, &
1358 & lbi, ubi, lbj, ubj, &
1359 & fi(:,:,icqcon))
1360
1361 CALL exchange_r2d_tile (ng, tile, &
1362 & lbi, ubi, lbj, ubj, &
1363 & fi(:,:,icqrhs))
1364
1365 CALL exchange_r2d_tile (ng, tile, &
1366 & lbi, ubi, lbj, ubj, &
1367 & snow)
1368# endif
1369
1370 CALL exchange_r2d_tile (ng, tile, &
1371 & lbi, ubi, lbj, ubj, &
1372 & sustr_ai)
1373
1374 CALL exchange_r2d_tile (ng, tile, &
1375 & lbi, ubi, lbj, ubj, &
1376 & svstr_ai)
1377# endif
1378
1379# ifdef EMINUSP
1380 CALL exchange_r2d_tile (ng, tile, &
1381 & lbi, ubi, lbj, ubj, &
1382 & evap)
1383# ifdef SALINITY
1384 CALL exchange_r2d_tile (ng, tile, &
1385 & lbi, ubi, lbj, ubj, &
1386 & stflux(:,:,isalt))
1387# endif
1388# endif
1389 CALL exchange_u2d_tile (ng, tile, &
1390 & lbi, ubi, lbj, ubj, &
1391 & sustr)
1392
1393 CALL exchange_v2d_tile (ng, tile, &
1394 & lbi, ubi, lbj, ubj, &
1395 & svstr)
1396# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1397 CALL exchange_u2d_tile (ng, tile, &
1398 & lbi, ubi, lbj, ubj, &
1399 & sustr_ao)
1400 CALL exchange_v2d_tile (ng, tile, &
1401 & lbi, ubi, lbj, ubj, &
1402 & svstr_ao)
1403# endif
1404 END IF
1405
1406# ifdef DISTRIBUTE
1407!
1408 CALL mp_exchange2d (ng, tile, inlm, 4, &
1409 & lbi, ubi, lbj, ubj, &
1410 & nghostpoints, &
1411 & ewperiodic(ng), nsperiodic(ng), &
1412 & lrflx, lhflx, shflx, &
1413 & stflux(:,:,itemp))
1414
1415# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1416 CALL mp_exchange2d (ng, tile, inlm, 4, &
1417 & lbi, ubi, lbj, ubj, &
1418 & nghostpoints, &
1419 & ewperiodic(ng), nsperiodic(ng), &
1420 & srflx, qnet_ao, sustr_ai, svstr_ai)
1421
1422# ifdef ICE_THERMO
1423 CALL mp_exchange2d (ng, tile, inlm, 4, &
1424 & lbi, ubi, lbj, ubj, &
1425 & nghostpoints, &
1426 & ewperiodic(ng), nsperiodic(ng), &
1427 fi(:,:,icqcon), &
1428 & fi(:,:,icqrhs), &
1429 & qnet_ai, snow)
1430# endif
1431# endif
1432
1433# ifdef EMINUSP
1434 CALL mp_exchange2d (ng, tile, inlm, 1, &
1435 & lbi, ubi, lbj, ubj, &
1436 & nghostpoints, &
1437 & ewperiodic(ng), nsperiodic(ng), &
1438 & evap)
1439# ifdef SALINITY
1440 CALL mp_exchange2d (ng, tile, inlm, 1, &
1441 & lbi, ubi, lbj, ubj, &
1442 & nghostpoints, &
1443 & ewperiodic(ng), nsperiodic(ng), &
1444 & stflux(:,:,isalt))
1445# endif
1446# endif
1447 CALL mp_exchange2d (ng, tile, inlm, 2, &
1448 & lbi, ubi, lbj, ubj, &
1449 & nghostpoints, &
1450 & ewperiodic(ng), nsperiodic(ng), &
1451 & sustr, svstr)
1452
1453# if defined ICE_MODEL && defined ICE_BULK_FLUXES
1454 CALL mp_exchange2d (ng, tile, inlm, 2, &
1455 & lbi, ubi, lbj, ubj, &
1456 & nghostpoints, &
1457 & ewperiodic(ng), nsperiodic(ng), &
1458 & sustr_ao, svstr_ao)
1459# endif
1460# endif
1461!
1462 RETURN

References mod_scalars::blk_beta, mod_scalars::blk_cpa, mod_scalars::blk_cpw, mod_scalars::blk_dter, mod_scalars::blk_rgas, mod_scalars::blk_tcw, mod_scalars::blk_visw, mod_scalars::blk_zabl, mod_scalars::blk_zq, mod_scalars::blk_zt, mod_scalars::blk_zw, bulk_psit(), bulk_psiu(), mod_ice::cd_ai, mod_scalars::cp, mod_param::domain, mod_scalars::emmiss, mod_scalars::ewperiodic, exchange_2d_mod::exchange_r2d_tile(), exchange_2d_mod::exchange_u2d_tile(), exchange_2d_mod::exchange_v2d_tile(), mod_scalars::g, mod_ice::icisst, mod_ice::icqcon, mod_ice::icqrhs, mod_scalars::iic, mod_param::inlm, mod_ice::isaice, mod_scalars::isalt, mod_ice::ishice, mod_ice::ishmel, mod_ice::ishsno, mod_scalars::itemp, mod_param::lm, mod_ice::min_ai, mod_param::mm, mp_exchange_mod::mp_exchange2d(), mod_param::nghostpoints, mod_scalars::nsperiodic, mod_scalars::ntstart, mod_scalars::perfectrst, mod_scalars::pi, mod_scalars::rho0, mod_scalars::rhow, mod_ice::snowdryrho, mod_scalars::stefbo, and mod_scalars::vonkar.

Referenced by bulk_flux().

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

◆ bulk_psit()

real(r8) function, public bulk_flux_mod::bulk_psit ( real(r8), intent(in) zol,
real(dp), intent(in) pi )

Definition at line 1530 of file bulk_flux.F.

1531!
1532!=======================================================================
1533! !
1534! This function evaluates the stability function for moisture and !
1535! heat by matching Kansas and free convection forms. The convective !
1536! form follows Fairall et al. (1996) with profile constants from !
1537! Grachev et al. (2000) BLM. The stable form is from Beljaars and !
1538! and Holtslag (1991). !
1539!
1540!=======================================================================
1541! !
1542!
1543 USE mod_kinds
1544!
1545! Function result
1546!
1547 real(r8) :: bulk_psit
1548!
1549! Imported variable declarations.
1550!
1551 real(dp), intent(in) :: pi
1552 real(r8), intent(in) :: ZoL
1553!
1554! Local variable declarations.
1555!
1556 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
1557
1558 real(r8) :: Fw, cff, psic, psik, x, y
1559!
1560!-----------------------------------------------------------------------
1561! Compute stability function, PSI.
1562!-----------------------------------------------------------------------
1563!
1564! Unstable conditions.
1565!
1566 IF (zol.lt.0.0_r8) THEN
1567 x=(1.0_r8-15.0_r8*zol)**0.5_r8
1568 psik=2.0_r8*log(0.5_r8*(1.0_r8+x))
1569!
1570! For very unstable conditions, use free-convection (Fairall).
1571!
1572 cff=sqrt(3.0_r8)
1573 y=(1.0_r8-34.15_r8*zol)**r3
1574 psic=1.5_r8*log(r3*(1.0_r8+y+y*y))- &
1575 & cff*atan((1.0_r8+2.0_r8*y)/cff)+pi/cff
1576!
1577! Match Kansas and free-convection forms with weighting Fw.
1578!
1579 cff=zol*zol
1580 fw=cff/(1.0_r8+cff)
1581 bulk_psit=(1.0_r8-fw)*psik+fw*psic
1582!
1583! Stable conditions.
1584!
1585 ELSE
1586 cff=min(50.0_r8,0.35_r8*zol)
1587 bulk_psit=-((1.0_r8+2.0_r8*zol)**1.5_r8+ &
1588 & 0.6667_r8*(zol-14.28_r8)/exp(cff)+8.525_r8)
1589 END IF
1590!
1591 RETURN
mod_kinds
Definition mod_kinds.F:2

References bulk_psit(), and mod_scalars::pi.

Referenced by ad_bulk_flux_mod::ad_bulk_flux_tile(), bulk_flux_tile(), bulk_psit(), rp_bulk_flux_mod::rp_bulk_flux_tile(), and tl_bulk_flux_mod::tl_bulk_flux_tile().

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

◆ bulk_psiu()

real(r8) function, public bulk_flux_mod::bulk_psiu ( real(r8), intent(in) zol,
real(dp), intent(in) pi )

Definition at line 1465 of file bulk_flux.F.

1466!
1467!=======================================================================
1468! !
1469! This function evaluates the stability function for wind speed !
1470! by matching Kansas and free convection forms. The convective !
1471! form follows Fairall et al. (1996) with profile constants from !
1472! Grachev et al. (2000) BLM. The stable form is from Beljaars !
1473! and Holtslag (1991). !
1474! !
1475!=======================================================================
1476!
1477 USE mod_kinds
1478!
1479! Function result
1480!
1481 real(r8) :: bulk_psiu
1482!
1483! Imported variable declarations.
1484!
1485 real(dp), intent(in) :: pi
1486 real(r8), intent(in) :: ZoL
1487!
1488! Local variable declarations.
1489!
1490 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
1491
1492 real(r8) :: Fw, cff, psic, psik, x, y
1493!
1494!-----------------------------------------------------------------------
1495! Compute stability function, PSI.
1496!-----------------------------------------------------------------------
1497!
1498! Unstable conditions.
1499!
1500 IF (zol.lt.0.0_r8) THEN
1501 x=(1.0_r8-15.0_r8*zol)**0.25_r8
1502 psik=2.0_r8*log(0.5_r8*(1.0_r8+x))+ &
1503 & log(0.5_r8*(1.0_r8+x*x))- &
1504 & 2.0_r8*atan(x)+0.5_r8*pi
1505!
1506! For very unstable conditions, use free-convection (Fairall).
1507!
1508 cff=sqrt(3.0_r8)
1509 y=(1.0_r8-10.15_r8*zol)**r3
1510 psic=1.5_r8*log(r3*(1.0_r8+y+y*y))- &
1511 & cff*atan((1.0_r8+2.0_r8*y)/cff)+pi/cff
1512!
1513! Match Kansas and free-convection forms with weighting Fw.
1514!
1515 cff=zol*zol
1516 fw=cff/(1.0_r8+cff)
1517 bulk_psiu=(1.0_r8-fw)*psik+fw*psic
1518!
1519! Stable conditions.
1520!
1521 ELSE
1522 cff=min(50.0_r8,0.35_r8*zol)
1523 bulk_psiu=-((1.0_r8+zol)+0.6667_r8*(zol-14.28_r8)/ &
1524 & exp(cff)+8.525_r8)
1525 END IF
1526!
1527 RETURN

References bulk_psiu(), and mod_scalars::pi.

Referenced by ad_bulk_flux_mod::ad_bulk_flux_tile(), bulk_flux_tile(), bulk_psiu(), rp_bulk_flux_mod::rp_bulk_flux_tile(), and tl_bulk_flux_mod::tl_bulk_flux_tile().

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