!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module vmix_kpp 6,18
!BOP
! !MODULE: vmix_kpp
!
! !DESCRIPTION:
! This module contains routines for initializing and computing
! vertical mixing coefficients for the KPP parameterization
! (see Large, McWilliams and Doney, Reviews of Geophysics, 32, 363
! November 1994).
!
! !REVISION HISTORY:
! SVN:$Id: vmix_kpp.F90 22881 2010-05-11 04:23:39Z njn01 $
! !USES:
use kinds_mod
use blocks
use distribution
use domain_size
use domain
use constants
use grid
use broadcast
use io
use state_mod
use exit_mod
use sw_absorption
use tavg, only: define_tavg_field, tavg_requested, accumulate_tavg_field
use io_types, only: stdout
use communicate, only: my_task, master_task
use tidal_mixing, only: TIDAL_COEF, tidal_mix_max, ltidal_mixing
use registry
use prognostic
use time_management
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: init_vmix_kpp, &
vmix_coeffs_kpp, &
add_kpp_sources, &
smooth_hblt, &
linertial
! !PUBLIC DATA MEMBERS:
real (r8), dimension(:,:,:), allocatable, public :: &
HMXL, &! mixed layer depth
KPP_HBLT, &! boundary layer depth
BOLUS_SP ! scaled eddy-induced (bolus) speed used in inertial
! mixing parameterization
real (r8), public :: &
bckgrnd_vdc2 ! variation in diffusivity
!EOP
!BOC
!-----------------------------------------------------------------------
!
! mixing constants
!
!-----------------------------------------------------------------------
real (r8) :: &
rich_mix ! coefficient for rich number term
real (r8), dimension(:,:,:,:), allocatable :: &
bckgrnd_vvc, &! background value for viscosity
bckgrnd_vdc ! background value for diffusivity
logical (log_kind) :: &
lrich, &! flag for computing Ri-dependent mixing
ldbl_diff, &! flag for computing double-diffusive mixing
lshort_wave, &! flag for computing short-wave forcing
lcheckekmo, &! check Ekman, Monin-Obhukov depth limit
llangmuir, &! flag for using Langmuir parameterization
linertial, &! flag for using inertial mixing parameterization
lccsm_control_compatible !flag for backwards compatibility with ccsm4 control
integer (int_kind) :: &
num_v_smooth_Ri ! num of times to vertically smooth Ri
real (r8), parameter :: &
epssfc = 0.1_r8 ! non-dimensional extent of sfc layer
real (r8) :: &
Prandtl ! Prandtl number
real (r8), dimension(:,:,:), allocatable :: &
FSTOKES ! ratio of stokes velocity to ustar used in Langmuir
! parameterization
!-----------------------------------------------------------------------
!
! non-local mixing (counter-gradient mixing), treated as source term
!
!-----------------------------------------------------------------------
real (r8), dimension(:,:,:,:,:), allocatable :: &
KPP_SRC ! non-local mixing (treated as source term)
!-----------------------------------------------------------------------
!
! parameters for subroutine bldepth: computes bndy layer depth
!
! concv = ratio of interior buoyancy frequency to
! buoyancy frequency at entrainment depth
! parameter statement sets the minimum value.
!
!-----------------------------------------------------------------------
real (r8), dimension(km) :: &
Ricr ! crit Rich no. for bndy layer depth
! as a function of vertical resolution
real (r8), parameter :: &
cekman = 0.7_r8, &! coefficient for Ekman depth
cmonob = 1.0_r8, &! coefficient for Monin-Obukhov depth
concv = 1.7_r8, &! minimum allowed value
hbf = 1.0_r8 ! frac of layer affected by sol forcing
!-----------------------------------------------------------------------
!
! parameters for subroutine ri_iwmix which computes
! vertical mixing coefficients below boundary layer due to shear
! instabilities, internal waves and convection
!
!-----------------------------------------------------------------------
real (r8), parameter :: &
Riinfty = 0.8_r8, &! Rich. no. limit for shear instab.
BVSQcon = c0 ! Brunt-Vaisala square cutoff(s**-2)
!-----------------------------------------------------------------------
!
! parameters for subroutine ddmix (double-diffusive mixing)
!
!-----------------------------------------------------------------------
real (r8), parameter :: &
Rrho0 = 2.55_r8, &! limit for double-diff density ratio
dsfmax = 1.0_r8 ! max diffusivity for salt fingering
!-----------------------------------------------------------------------
!
! parameters for subroutine blmix: mixing within boundary layer
!
! cstar = proportionality coefficient for nonlocal transport
! cg = non-dimensional coefficient for counter-gradient term
!
!-----------------------------------------------------------------------
real (r8), parameter :: &
cstar = 10.0_r8 ! coeff for nonlocal transport
real (r8) :: &
cg, &! coefficient for counter-gradient term
Vtc ! resolution and buoyancy independent part of the
! turbulent velocity shear coefficient (for bulk Ri no)
!-----------------------------------------------------------------------
!
! parameters for velocity scale function (from Large et al.)
!
!-----------------------------------------------------------------------
real (r8), parameter :: &
zeta_m = -0.2_r8, &
zeta_s = -1.0_r8, &
c_m = 8.38_r8, &
c_s = 98.96_r8, &
a_m = 1.26_r8, &
a_s = -28.86_r8
!-----------------------------------------------------------------------
!
! common vertical grid arrays used by KPP mixing
!
!-----------------------------------------------------------------------
real (r8), dimension(:), allocatable :: &
zgrid, &! depth at cell interfaces
hwide ! layer thickness at interfaces
!-----------------------------------------------------------------------
!
! ids for tavg diagnostics computed in this module
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_QSW_HBL, &! tavg id for solar short-wave heat flux in bndry layer
tavg_VDC_BCK, &! tavg id for bckgrnd vertical tracer diffusivity
tavg_VVC_BCK, &! tavg id for bckgrnd vertical momentum viscosity
tavg_KVMIX, &! tavg id for tidal+bckgrnd vertical tracer diffusivity
tavg_KVMIX_M, &! tavg id for tidal+bckgrnd vertical momentum viscosity
tavg_TPOWER
integer (int_kind), dimension(nt) :: &
tavg_KPP_SRC ! tavg id for KPP_SRC for each tracer
real (r8), dimension(:,:,:,:), allocatable :: &
TIDAL_DIFF ! diffusivity due to tidal mixing
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_vmix_kpp
! !INTERFACE:
subroutine init_vmix_kpp(VDC, VVC) 1,29
! !DESCRIPTION:
! Initializes constants and reads options for the KPP parameterization.
!
! !REVISION HISTORY:
! same as module
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(:,:,:,:), intent(inout) :: &
VVC ! viscosity for momentum diffusion
real (r8), dimension(:,:,0:,:,:),intent(inout) :: &
VDC ! diffusivity for tracer diffusion
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables and namelist inputs
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
n, &! local dummy index for tracer
k, &! local dummy index for vertical lvl
i, j, iblock, &! local dummy indexes
nml_error ! namelist i/o error flag
real (r8) :: &
bckgrnd_vdc1, &! background diffusivity (Ledwell)
bckgrnd_vdc_eq, &! equatorial diffusivity (Gregg)
bckgrnd_vdc_psim, &! Max. PSI induced diffusivity (MacKinnon)
bckgrnd_vdc_psin, &! PSI diffusivity in northern hemisphere
bckgrnd_vdc_psis, &! PSI diffusivity in southern hemisphere
bckgrnd_vdc_ban, &! Banda Sea diffusivity (Gordon)
bckgrnd_vdc_dpth, &! depth at which diff equals vdc1
bckgrnd_vdc_linv, &! inverse length for transition region
tlatd, &! tlat * radian (degrees)
tlond ! tlon * radian (degrees)
logical (log_kind) :: &
lhoriz_varying_bckgrnd
namelist /vmix_kpp_nml/bckgrnd_vdc1, bckgrnd_vdc2, &
bckgrnd_vdc_eq, bckgrnd_vdc_psim, &
bckgrnd_vdc_ban, &
bckgrnd_vdc_dpth, bckgrnd_vdc_linv, &
Prandtl, rich_mix, &
num_v_smooth_Ri, lrich, ldbl_diff, &
lshort_wave, lcheckekmo, &
lhoriz_varying_bckgrnd, llangmuir, &
linertial
character (16), parameter :: &
fmt_real = '(a30,2x,1pe12.5)'
character (11), parameter :: &
fmt_log = '(a30,2x,l7)'
character (11), parameter :: &
fmt_int = '(a30,2x,i5)'
character (char_len) :: &
string, string2
!-----------------------------------------------------------------------
!
! set defaults for mixing coefficients, then read them from namelist
!
!-----------------------------------------------------------------------
bckgrnd_vdc1 = 0.1_r8
bckgrnd_vdc2 = c0
bckgrnd_vdc_eq = 0.01_r8
bckgrnd_vdc_psim = 0.13_r8
bckgrnd_vdc_ban = c1
bckgrnd_vdc_dpth = 2500.0e02_r8
bckgrnd_vdc_linv = 4.5e-05_r8
Prandtl = 10.0_r8
rich_mix = 50.0_r8
lrich = .true.
ldbl_diff = .false.
lshort_wave = .false.
lcheckekmo = .false.
lhoriz_varying_bckgrnd = .false.
llangmuir = .false.
linertial = .false.
num_v_smooth_Ri = 1
if (my_task == master_task) then
open (nml_in, file=nml_filename, status='old',iostat=nml_error)
if (nml_error /= 0) then
nml_error = -1
else
nml_error = 1
endif
do while (nml_error > 0)
read(nml_in, nml=vmix_kpp_nml, iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
endif
call broadcast_scalar(nml_error, master_task)
if (nml_error /= 0) then
call exit_POP(sigAbort,'ERROR (init_vmix_kpp) reading vmix_kpp_nml')
endif
if (my_task == master_task) then
write(stdout,delim_fmt)
write(stdout,blank_fmt)
write(stdout,*) ' vmix_kpp_nml namelist settings:'
write(stdout,blank_fmt)
write(stdout,vmix_kpp_nml)
write(stdout,blank_fmt)
write(stdout,fmt_real) ' bckgrnd_vdc1 =', bckgrnd_vdc1
write(stdout,fmt_real) ' bckgrnd_vdc2 =', bckgrnd_vdc2
write(stdout,fmt_real) ' bckgrnd_vdc_dpth =', bckgrnd_vdc_dpth
write(stdout,fmt_real) ' bckgrnd_vdc_linv =', bckgrnd_vdc_linv
write(stdout,fmt_real) ' bckgrnd_vdc_eq =', bckgrnd_vdc_eq
write(stdout,fmt_real) ' bckgrnd_vdc_psim =', bckgrnd_vdc_psim
write(stdout,fmt_real) ' bckgrnd_vdc_ban =', bckgrnd_vdc_ban
write(stdout,fmt_real) ' Prandtl =', Prandtl
write(stdout,fmt_real) ' rich_mix =', rich_mix
write(stdout,fmt_log ) ' Ri mixing =', lrich
write(stdout,fmt_log ) ' double-diff =', ldbl_diff
write(stdout,fmt_log ) ' short_wave =', lshort_wave
write(stdout,fmt_log ) ' lcheckekmo =', lcheckekmo
write(stdout,fmt_int ) ' num_smooth_Ri =', num_v_smooth_Ri
write(stdout,fmt_log ) ' lhoriz_varying_bckgrnd =', lhoriz_varying_bckgrnd
write(stdout,fmt_log ) ' langmuir parameterization =', llangmuir
write(stdout,fmt_log ) ' inertial mixing param. =', linertial
endif
call broadcast_scalar(bckgrnd_vdc1, master_task)
call broadcast_scalar(bckgrnd_vdc2, master_task)
call broadcast_scalar(bckgrnd_vdc_dpth, master_task)
call broadcast_scalar(bckgrnd_vdc_linv, master_task)
call broadcast_scalar(bckgrnd_vdc_eq , master_task)
call broadcast_scalar(bckgrnd_vdc_psim, master_task)
call broadcast_scalar(bckgrnd_vdc_ban , master_task)
call broadcast_scalar(Prandtl, master_task)
call broadcast_scalar(rich_mix, master_task)
call broadcast_scalar(num_v_smooth_Ri, master_task)
call broadcast_scalar(lrich, master_task)
call broadcast_scalar(ldbl_diff, master_task)
call broadcast_scalar(lshort_wave, master_task)
call broadcast_scalar(lcheckekmo, master_task)
call broadcast_scalar(lhoriz_varying_bckgrnd,master_task)
call broadcast_scalar(llangmuir, master_task)
call broadcast_scalar(linertial, master_task)
!-----------------------------------------------------------------------
!
! determine if this case must be backwards compatible with ccsm4 control
!
!-----------------------------------------------------------------------
lccsm_control_compatible = registry_match('lccsm_control_compatible')
!-----------------------------------------------------------------------
!
! define some non-dimensional constants
!
!-----------------------------------------------------------------------
Vtc = sqrt(0.2_r8/c_s/epssfc)/vonkar**2
cg = cstar*vonkar*(c_s*vonkar*epssfc)**p33
!-----------------------------------------------------------------------
!
! define vertical grid coordinates and cell widths
! compute horizontally or vertically varying background
! (internal wave) diffusivity and viscosity
!
! the vertical profile has the functional form
!
! BCKGRND_VDC(z) = vdc1 + vdc2*atan((|z| - dpth)/L) or
! = vdc1 + vdc2*atan((|z| - dpth)*linv)
!
! where
!
! vdc1 = vertical diffusivity at |z|=D (cm^2/s)
! vdc2 = amplitude of variation (cm^2/s)
! linv = inverse length scale for transition region (1/cm)
! dpth = depth where diffusivity is vdc1 (cm)
!
! the viscosity has the same form but multiplied by a constant
! Prandtl number
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! if tidal mixing is on, use vertically constant (internal wave)
! background diffusivity and viscosity (see consistency testing in POP_check)
!
!-----------------------------------------------------------------------
if (bckgrnd_vdc2 /= c0 .and. my_task == master_task) then
write (stdout,blank_fmt)
write (stdout,'(a43)') &
'Vertical Profile for background diffusivity'
endif
!-----------------------------------------------------------------------
!
! error checking
!
!-----------------------------------------------------------------------
if (lhoriz_varying_bckgrnd .and. bckgrnd_vdc2 /= c0) then
call exit_POP (sigAbort, &
'ERROR (init_vmix_kpp): lhoriz_varying_bckgrnd .and. bckgrnd_vdc2 /= c0')
endif
!-----------------------------------------------------------------------
!
! initialize grid info (only need one block since the vertical grid is
! the same across blocks)
!
!-----------------------------------------------------------------------
allocate (zgrid(0:km+1), hwide(0:km+1))
zgrid(0) = eps
hwide(0) = eps
do k=1,km
zgrid(k) = -zt(k)
hwide(k) = dz(k)
enddo
zgrid(km+1) = -zw(km)
hwide(km+1) = eps
allocate (bckgrnd_vvc(nx_block,ny_block,km,nblocks_clinic))
allocate (bckgrnd_vdc(nx_block,ny_block,km,nblocks_clinic))
if (lhoriz_varying_bckgrnd) then
k = 1
do iblock=1,nblocks_clinic
do j=1,ny_block
do i=1,nx_block
tlatd = TLAT(i,j,iblock)*radian
tlond = TLON(i,j,iblock)*radian
bckgrnd_vdc_psis= bckgrnd_vdc_psim*exp(-(0.4_r8*(tlatd+28.9_r8))**c2)
bckgrnd_vdc_psin= bckgrnd_vdc_psim*exp(-(0.4_r8*(tlatd-28.9_r8))**c2)
bckgrnd_vdc(i,j,k,iblock)=bckgrnd_vdc_eq+bckgrnd_vdc_psin+bckgrnd_vdc_psis
if ( tlatd .lt. -10.0_r8 ) then
bckgrnd_vdc(i,j,k,iblock) = bckgrnd_vdc(i,j,k,iblock) + bckgrnd_vdc1
elseif ( tlatd .le. 10.0_r8 ) then
bckgrnd_vdc(i,j,k,iblock) = bckgrnd_vdc(i,j,k,iblock) + &
bckgrnd_vdc1 * (tlatd/10.0_r8)**c2
else
bckgrnd_vdc(i,j,k,iblock)=bckgrnd_vdc(i,j,k,iblock) + bckgrnd_vdc1
endif
!----------------
! North Banda Sea
!----------------
if ( (tlatd .lt. -1.0_r8) .and. (tlatd .gt. -4.0_r8) .and. &
(tlond .gt. 103.0_r8) .and. (tlond .lt. 134.0_r8)) then
bckgrnd_vdc(i,j,k,iblock) = bckgrnd_vdc_ban
endif
!-----------------
! Middle Banda Sea
!-----------------
if ( (tlatd .le. -4.0_r8) .and. (tlatd .gt. -7.0_r8) .and. &
(tlond .gt. 106.0_r8) .and. (tlond .lt. 140.0_r8)) then
bckgrnd_vdc(i,j,k,iblock) = bckgrnd_vdc_ban
endif
!----------------
! South Banda Sea
!----------------
if ( (tlatd .le. -7.0_r8) .and. (tlatd .gt. -8.3_r8) .and. &
(tlond .gt. 111.0_r8) .and. (tlond .lt. 142.0_r8)) then
bckgrnd_vdc(i,j,k,iblock) = bckgrnd_vdc_ban
endif
bckgrnd_vvc(i,j,k,iblock) = Prandtl*bckgrnd_vdc(i,j,k,iblock)
end do ! i
end do ! j
end do ! iblock
do k=2,km
bckgrnd_vdc(:,:,k,:) = bckgrnd_vdc(:,:,1,:)
bckgrnd_vvc(:,:,k,:) = bckgrnd_vvc(:,:,1,:)
enddo
else
!-----------------------------------------------------------------------
!
! only need one block since the vertical grid is the same across blocks
!
!-----------------------------------------------------------------------
do k=1,km
bckgrnd_vdc(:,:,k,:) = bckgrnd_vdc1 + bckgrnd_vdc2* &
atan(bckgrnd_vdc_linv* &
(zw(k)-bckgrnd_vdc_dpth))
bckgrnd_vvc(:,:,k,:) = Prandtl*bckgrnd_vdc(:,:,k,:)
if (bckgrnd_vdc2 /= c0 .and. my_task == master_task) then
write (stdout,'(2x,e12.6)') bckgrnd_vdc(1,1,k,1)
endif
end do
endif ! lhoriz_varying_bckgrnd
!-----------------------------------------------------------------------
!
! compute crit Rich number for bndy layer depth as a function
! of model vertical resolution. note that hwide is in cm.
!
!-----------------------------------------------------------------------
Ricr(1:km) = 0.3_r8
!-----------------------------------------------------------------------
!
! allocate and initialize KPP-specific arrays
!
!-----------------------------------------------------------------------
allocate (HMXL (nx_block,ny_block,nblocks_clinic), &
KPP_HBLT (nx_block,ny_block,nblocks_clinic), &
KPP_SRC (nx_block,ny_block,km,nt,nblocks_clinic))
HMXL = c0
KPP_HBLT = c0
KPP_SRC = c0
VDC = c0
VVC = c0
if ( ltidal_mixing ) then
allocate ( TIDAL_DIFF(nx_block,ny_block,km,nblocks_clinic) )
TIDAL_DIFF = c0
endif
!-----------------------------------------------------------------------
!
! allocate and initialize the eddy-induced speed array
!
!-----------------------------------------------------------------------
if ( linertial ) then
allocate ( BOLUS_SP(nx_block,ny_block,nblocks_clinic) )
BOLUS_SP = c0
endif
!-----------------------------------------------------------------------
!
! allocate and initialize ratio of stokes velocity to ustar
!
!-----------------------------------------------------------------------
allocate (FSTOKES(nx_block,ny_block,nblocks_clinic))
do iblock=1,nblocks_clinic
FSTOKES(:,:,iblock) = 11._r8 - MAX( c5*cos(c3*TLAT(:,:,iblock)) , c0 )
enddo
!-----------------------------------------------------------------------
!
! define tavg fields computed from vmix_kpp module routines
!
!-----------------------------------------------------------------------
string = 'Solar Short-Wave Heat Flux in bndry layer'
call define_tavg_field(tavg_QSW_HBL,'QSW_HBL',2, &
long_name=trim(string), &
units='watt/m^2', grid_loc='2110', &
coordinates='TLONG TLAT time')
if (ltidal_mixing) then
string = 'Vertical diabatic diffusivity due to Tidal Mixing + background'
else
string = 'Vertical diabatic diffusivity due to background'
endif
call define_tavg_field(tavg_KVMIX,'KVMIX',3, &
long_name=trim(string), &
units='centimeter^2/s', &
grid_loc='3113', &
coordinates ='TLONG TLAT z_w_bot time' )
string = 'Vertical diabatic diffusivity due to background'
call define_tavg_field(tavg_VDC_BCK,'VDC_BCK',3, &
long_name=trim(string), &
units='centimeter^2/s', &
grid_loc='3113', &
coordinates ='TLONG TLAT z_w_bot time' )
if (ltidal_mixing) then
string = 'Vertical viscosity due to Tidal Mixing + background'
else
string = 'Vertical viscosity due to background'
endif
call define_tavg_field(tavg_KVMIX_M,'KVMIX_M',3, &
long_name=trim(string), &
units='centimeter^2/s', &
grid_loc='3113', &
coordinates ='TLONG TLAT z_w_bot time' )
string = 'Vertical viscosity due to background'
call define_tavg_field(tavg_VVC_BCK,'VVC_BCK',3, &
long_name=trim(string), &
units='centimeter^2/s', &
grid_loc='3113', &
coordinates ='TLONG TLAT z_w_bot time' )
string = 'Energy Used by Vertical Mixing'
call define_tavg_field(tavg_TPOWER,'TPOWER',3, &
long_name=trim(string), &
units='erg/s', &
grid_loc='3113', &
coordinates ='TLONG TLAT z_w_bot time' )
do n = 1,nt
string = 'KPP_SRC_' /&
&/ trim(tracer_d(n)%short_name)
string2 = trim(tracer_d(n)%short_name) /&
&/ ' tendency from KPP non local mixing term'
call define_tavg_field(tavg_KPP_SRC(n),trim(string),3, &
long_name=trim(string2), &
units=trim(tracer_d(n)%tend_units), &
scale_factor=tracer_d(n)%scale_factor,&
grid_loc='3111', &
coordinates ='TLONG TLAT z_t time' )
enddo
!-----------------------------------------------------------------------
!EOC
call flushm (stdout)
end subroutine init_vmix_kpp
!***********************************************************************
!BOP
! !IROUTINE: vmix_coeffs_kpp
! !INTERFACE:
subroutine vmix_coeffs_kpp(VDC, VVC, TRCR, UUU, VVV, RHOMIX, STF, SHF_QSW, & 2,8
this_block, convect_diff, convect_visc, &
SMF, SMFT)
! !DESCRIPTION:
! This is the main driver routine which calculates the vertical
! mixing coefficients for the KPP mixing scheme as outlined in
! Large, McWilliams and Doney, Reviews of Geophysics, 32, 363
! (November 1994). The non-local mixing is also computed here, but
! is treated as a source term in baroclinic.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TRCR ! tracers at current time
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UUU,VVV ! velocities at current time
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
RHOMIX ! density at mix time
real (r8), dimension(nx_block,ny_block,nt), intent(in) :: &
STF ! surface forcing for all tracers
real (r8), dimension(nx_block,ny_block), intent(in) :: &
SHF_QSW ! short-wave forcing
real (r8), dimension(nx_block,ny_block,2), intent(in), optional :: &
SMF, &! surface momentum forcing at U points
SMFT ! surface momentum forcing at T points
! *** either one or the other (not
! *** both) should be passed
real (r8), intent(in) :: &
convect_diff, &! diffusivity to mimic convection
convect_visc ! viscosity to mimic convection
type (block), intent(in) :: &
this_block ! block information for current block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km), intent(inout) :: &
VVC ! viscosity for momentum diffusion
real (r8), dimension(nx_block,ny_block,0:km+1,2),intent(inout) :: &
VDC ! diffusivity for tracer diffusion
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, &! vertical level index
i,j, &! horizontal loop indices
n, &! tracer index
mt2, &! index for separating temp from other trcrs
bid ! local block address for this block
integer (int_kind), dimension(nx_block,ny_block) :: &
KBL ! index of first lvl below hbl
real (r8), dimension(nx_block,ny_block) :: &
USTAR, &! surface friction velocity
BFSFC, &! surface buoyancy forcing
WORK1,WORK2,&! temporary storage
FCON, &! convection temporary
STABLE ! = 1 for stable forcing; = 0 for unstable forcing
real (r8), dimension(nx_block,ny_block,km) :: &
DBLOC, &! buoyancy difference between adjacent levels
DBSFC, &! buoyancy difference between level and surface
GHAT ! non-local mixing coefficient
real (r8), dimension(nx_block,ny_block,0:km+1) :: &
VISC ! local temp for viscosity
!-----------------------------------------------------------------------
!
! initialize
!
!-----------------------------------------------------------------------
bid = this_block%local_id
if (.not. present(SMF) .and. .not. present(SMFT)) &
call exit_POP(sigAbort, &
'ERROR KPP: must supply either SMF or SMFT')
!-----------------------------------------------------------------------
!
! compute buoyancy differences at each vertical level.
!
!-----------------------------------------------------------------------
call buoydiff(DBLOC, DBSFC, TRCR, this_block)
!-----------------------------------------------------------------------
!
! compute mixing due to shear instability, internal waves and
! convection
!
!-----------------------------------------------------------------------
call ri_iwmix(DBLOC, VISC, VDC, UUU, VVV, RHOMIX,&
convect_diff, convect_visc, this_block)
!-----------------------------------------------------------------------
!
! compute double diffusion if desired
!
!-----------------------------------------------------------------------
if (ldbl_diff) call ddmix(VDC, TRCR, this_block)
!-----------------------------------------------------------------------
!
! compute boundary layer depth
!
!-----------------------------------------------------------------------
if (present(SMFT)) then
call bldepth (DBLOC, DBSFC, TRCR, UUU, VVV, STF, SHF_QSW, &
KPP_HBLT(:,:,bid), USTAR, BFSFC, STABLE, KBL, &
this_block, SMFT=SMFT)
else
call bldepth (DBLOC, DBSFC, TRCR, UUU, VVV, STF, SHF_QSW, &
KPP_HBLT(:,:,bid), USTAR, BFSFC, STABLE, KBL, &
this_block, SMF=SMF)
endif
!-----------------------------------------------------------------------
!
! compute boundary layer diffusivities
!
!-----------------------------------------------------------------------
call blmix(VISC, VDC, KPP_HBLT(:,:,bid), USTAR, BFSFC, STABLE, &
KBL, GHAT, this_block)
!-----------------------------------------------------------------------
!
! consider interior convection:
!
! compute function of Brunt-Vaisala squared for convection.
!
! use either a smooth
!
! WORK1 = N**2, FCON is function of N**2
! FCON = 0 for N**2 > 0
! FCON = [1-(1-WORK1/BVSQcon)**2]**3 for BVSQcon < N**2 < 0
! FCON = 1 for N**2 < BVSQcon
!
! or a step function. The smooth function has been used with
! BVSQcon = -0.2e-4_dbl_kind.
!
! after convection, average viscous coeffs to U-grid and reset sea
! floor values
!
!-----------------------------------------------------------------------
do k=1,km-1
if (partial_bottom_cells) then
WORK1 = DBLOC(:,:,k)/(p5*(DZT(:,:,k ,bid) + &
DZT(:,:,k+1,bid)))
else
WORK1 = DBLOC(:,:,k)/(zgrid(k) - zgrid(k+1))
end if
if (BVSQcon /= c0) then
WORK2 = min(c1-(max(WORK1,BVSQcon))/BVSQcon, c1)
FCON = (c1 - WORK2*WORK2)**3
else
where (WORK1 > c0)
FCON = c0
elsewhere
FCON = c1
end where
endif
!*** add convection and reset sea floor values to zero
do j=1,ny_block
do i=1,nx_block
if ( k >= KBL(i,j) ) then
VISC(i,j,k) = VISC(i,j,k) + convect_visc * FCON(i,j)
VDC(i,j,k,1) = VDC(i,j,k,1) + convect_diff * FCON(i,j)
VDC(i,j,k,2) = VDC(i,j,k,2) + convect_diff * FCON(i,j)
endif
if (k >= KMT(i,j,bid)) then
VISC(i,j,k ) = c0
VDC (i,j,k,1) = c0
VDC (i,j,k,2) = c0
endif
end do
end do
!*** now average visc to U grid
call tgrid_to_ugrid(WORK2,VISC(:,:,k),bid)
VVC(:,:,k) = merge(WORK2, c0, (k < KMU(:,:,bid)))
enddo
VDC(:,:,km,:) = c0
VVC(:,:,km) = c0
!-----------------------------------------------------------------------
!
! add ghatp term from previous computation to right-hand-side
! source term on current row
!
!-----------------------------------------------------------------------
do n=1,nt
mt2=min(n,2)
KPP_SRC(:,:,1,n,bid) = STF(:,:,n)/dz(1) &
*(-VDC(:,:,1,mt2)*GHAT(:,:,1))
if (partial_bottom_cells) then
do k=2,km
KPP_SRC(:,:,k,n,bid) = STF(:,:,n)/DZT(:,:,k,bid) &
*( VDC(:,:,k-1,mt2)*GHAT(:,:,k-1) &
-VDC(:,:,k ,mt2)*GHAT(:,:,k ))
enddo
else
do k=2,km
KPP_SRC(:,:,k,n,bid) = STF(:,:,n)/dz(k) &
*( VDC(:,:,k-1,mt2)*GHAT(:,:,k-1) &
-VDC(:,:,k ,mt2)*GHAT(:,:,k ))
enddo
endif
enddo
!-----------------------------------------------------------------------
!
! compute diagnostic mixed layer depth (cm) using a max buoyancy
! gradient criterion. Use USTAR and BFSFC as temps.
!
!-----------------------------------------------------------------------
USTAR = c0
where (KMT(:,:,bid) == 1)
HMXL(:,:,bid) = zt(1)
elsewhere
HMXL(:,:,bid) = c0
endwhere
if (partial_bottom_cells) then
do k=2,km
where (k <= KMT(:,:,bid))
STABLE = zt(k-1) + p5*(DZT(:,:,k-1,bid) + DZT(:,:,k,bid))
USTAR = max(DBSFC(:,:,k)/STABLE,USTAR)
HMXL(:,:,bid) = STABLE
endwhere
enddo
VISC(:,:,1) = c0
do k=2,km
where (USTAR > c0 )
VISC(:,:,k) = (DBSFC(:,:,k)-DBSFC(:,:,k-1))/ &
(p5*(DZT(:,:,k,bid) + DZT(:,:,k-1,bid)))
end where
where ( VISC(:,:,k) >= USTAR .and. &
(VISC(:,:,k)-VISC(:,:,k-1)) /= c0 .and. &
USTAR > c0 ) ! avoid divide by zero
BFSFC = (VISC(:,:,k) - USTAR)/ &
(VISC(:,:,k)-VISC(:,:,k-1))
! tqian
! HMXL(:,:,bid) = (zt(k-1) + p5*DZT(:,:,k-1,bid))*(c1-BFSFC) &
! + (zt(k-1) - p5*DZT(:,:,k-1,bid))*BFSFC
HMXL(:,:,bid) = (zt(k-1) + p25*(DZT(:,:,k-1,bid)+DZT(:,:,k,bid)))*(c1-BFSFC) &
+ (zt(k-1) - p25*(DZT(:,:,k-2,bid)+DZT(:,:,k-1,bid)))*BFSFC
USTAR(:,:) = c0
endwhere
enddo
else
do k=2,km
where (k <= KMT(:,:,bid))
USTAR = max(DBSFC(:,:,k)/zt(k),USTAR)
HMXL(:,:,bid) = zt(k)
endwhere
enddo
VISC(:,:,1) = c0
do k=2,km
where (USTAR > c0 )
VISC(:,:,k) = (DBSFC(:,:,k)-DBSFC(:,:,k-1))/ &
(zt(k) - zt(k-1))
end where
where ( VISC(:,:,k) >= USTAR .and. &
(VISC(:,:,k)-VISC(:,:,k-1)) /= c0 .and. &
USTAR > c0 ) ! avoid divide by zero
BFSFC = (VISC(:,:,k) - USTAR)/ &
(VISC(:,:,k)-VISC(:,:,k-1))
HMXL(:,:,bid) = -p5*(zgrid(k ) + zgrid(k-1))*(c1-BFSFC) &
-p5*(zgrid(k-1) + zgrid(k-2))*BFSFC
USTAR(:,:) = c0
endwhere
enddo
endif
!-----------------------------------------------------------------------
!EOC
end subroutine vmix_coeffs_kpp
!***********************************************************************
!BOP
! !IROUTINE: ri_iwmix
! !INTERFACE:
subroutine ri_iwmix(DBLOC, VISC, VDC, UUU, VVV, RHOMIX, & 1,11
convect_diff, convect_visc, this_block)
! !DESCRIPTION:
! Computes viscosity and diffusivity coefficients for the interior
! ocean due to shear instability (richardson number dependent),
! internal wave activity, and to static instability (Ri < 0).
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UUU ! U velocities at current time
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
VVV, &! V velocities at current time
DBLOC ! buoyancy difference between adjacent levels
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
RHOMIX ! density at mix time
real (r8), intent(in) :: &
convect_diff, &! diffusivity to mimic convection
convect_visc ! viscosity to mimic convection
type (block), intent(in) :: &
this_block ! block information for current block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,0:km+1,2), intent(inout) :: &
VDC ! diffusivity for tracer diffusion
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,0:km+1), intent(out) :: &
VISC ! viscosity
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, &! index for vertical levels
i,j, &! horizontal loop indices
bid, &! local block index
n ! vertical smoothing index
real (r8), dimension(nx_block,ny_block) :: &
KVMIX, &! vertical diffusivity
KVMIX_M ! vertical viscosity
real (r8), dimension(nx_block,ny_block) :: &
VSHEAR, &! (local velocity shear)^2
RI_LOC, &! local Richardson number
FRI, &! function of Ri for shear
FCON, &
WORK1
!-----------------------------------------------------------------------
!
! compute mixing at each level
!
!-----------------------------------------------------------------------
bid = this_block%local_id
KVMIX = c0
KVMIX_M = c0
VISC(:,:,0) = c0
do k = 1,km
!-----------------------------------------------------------------------
!
! compute velocity shear squared and average to T points:
! VSHEAR = (UUU(k)-UUU(k+1))**2+(VVV(k)-VVV(k+1))**2
! Use FRI here as a temporary.
!
!-----------------------------------------------------------------------
if (k < km) then
FRI = (UUU(:,:,k)-UUU(:,:,k+1))**2 + &
(VVV(:,:,k)-VVV(:,:,k+1))**2
if (partial_bottom_cells) then
FRI = FRI/(p5*(DZU(:,:,k ,bid) + DZU(:,:,k+1,bid)))**2
endif
call ugrid_to_tgrid(VSHEAR,FRI,bid)
else
VSHEAR = c0
endif
!-----------------------------------------------------------------------
!
! compute local richardson number
! use visc array as temporary Ri storage to be smoothed
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
if (k < km) then
RI_LOC = DBLOC(:,:,k)/ &
(VSHEAR + eps/(p5*(DZT(:,:,k ,bid) + &
DZT(:,:,k+1,bid)))**2)/ &
(p5*(DZT(:,:,k,bid) + DZT(:,:,k+1,bid)))
else
RI_LOC = DBLOC(:,:,k)/(VSHEAR + &
eps/(p5*DZT(:,:,k,bid))**2)/(p5*DZT(:,:,k,bid))
end if
else
RI_LOC = DBLOC(:,:,k)*(zgrid(k)-zgrid(k+1))/(VSHEAR + eps)
end if
VISC(:,:,k) = merge(RI_LOC, VISC(:,:,k-1), k <= KMT(:,:,bid))
enddo
!-----------------------------------------------------------------------
!
! vertically smooth Ri num_v_smooth_Ri times with 1-2-1 weighting
! result again stored temporarily in VISC and use RI_LOC and FRI
! as temps
!
!-----------------------------------------------------------------------
do n = 1,num_v_smooth_Ri
FRI = p25 * VISC(:,:,1)
VISC(:,:,km+1) = VISC(:,:,km)
do k=1,km
!DIR$ NODEP
do j=1,ny_block
!DIR$ NODEP
do i=1,nx_block
RI_LOC(i,j) = VISC(i,j,k)
if (KMT(i,j,bid) >= 3) then
VISC(i,j,k) = FRI(i,j) + p5*RI_LOC(i,j) + p25*VISC(i,j,k+1)
endif
FRI(i,j) = p25*RI_LOC(i,j)
enddo
enddo
enddo
enddo
!-----------------------------------------------------------------------
!
! now that we have a smoothed Ri field, finish computing coeffs
! at each level
!
!-----------------------------------------------------------------------
if ( ltidal_mixing ) TIDAL_DIFF(:,:,:,bid) = c0
do k = 1,km
!-----------------------------------------------------------------------
!
! if Ri-number mixing requested,
! evaluate function of Ri for shear instability:
! for 0 < Ri < Riinfty, function = (1 - (Ri/Riinfty)**2)**3
! for Ri > Riinfty, function = 0
! for Ri < 0 , function = 1
! compute contribution due to shear instability
! VISC holds smoothed Ri at k, but replaced by real VISC
!
! otherwise only use iw
! convection is added later
!
!-----------------------------------------------------------------------
if ( ltidal_mixing ) then
!-----------------------------------------------------------------------
!
! consider the internal wave mixing first. rich_mix is used as the
! upper limit for internal wave mixing coefficient. bckgrnd_vvc
! was already multiplied by Prandtl.
!
! NOTE: no partial_bottom_cell implementation at this time
!
!-----------------------------------------------------------------------
WORK1 = DBLOC(:,:,k)/(zgrid(k) - zgrid(k+1))
where (WORK1 > c0)
TIDAL_DIFF(:,:,k,bid) = TIDAL_COEF(:,:,k,bid)/WORK1
endwhere
if (.not. lccsm_control_compatible) then ! this step breaks backwards compatibility
where ( k > 2 .and. ( k == KMT(:,:,bid)-1 .or. k == KMT(:,:,bid)-2 ) )
TIDAL_DIFF(:,:,k,bid) = max( TIDAL_DIFF(:,:,k, bid), &
TIDAL_DIFF(:,:,k-1,bid) )
endwhere
endif
WORK1 = Prandtl*min(bckgrnd_vvc(:,:,k,bid)/Prandtl &
+ TIDAL_DIFF(:,:,k,bid), tidal_mix_max)
if ( k < km ) then
KVMIX_M(:,:) = WORK1(:,:)
endif
if ( k < km ) then
VDC(:,:,k,2) = min(bckgrnd_vdc(:,:,k,bid) + TIDAL_DIFF(:,:,k,bid), &
tidal_mix_max)
KVMIX(:,:) = VDC(:,:,k,2)
endif
if (lrich) then
FRI = min((max(VISC(:,:,k),c0))/Riinfty, c1)
VISC(:,:,k) = WORK1 + rich_mix*(c1 - FRI*FRI)**3
if ( k < km ) then
VDC(:,:,k,2) = VDC(:,:,k,2) + rich_mix*(c1 - FRI*FRI)**3
VDC(:,:,k,1) = VDC(:,:,k,2)
endif
else
VISC(:,:,k) = WORK1
if ( k < km ) then
VDC(:,:,k,1) = VDC(:,:,k,2)
endif
endif
else ! .not. ltidal_mixing
if ( k < km ) then
KVMIX(:,:) = bckgrnd_vdc(:,:,k,bid)
KVMIX_M(:,:) = bckgrnd_vvc(:,:,k,bid)
endif
if (lrich) then
FRI = min((max(VISC(:,:,k),c0))/Riinfty, c1)
VISC(:,:,k ) = bckgrnd_vvc(:,:,k,bid) + &
rich_mix*(c1 - FRI*FRI)**3
if ( k < km ) then
VDC (:,:,k,2) = bckgrnd_vdc(:,:,k,bid) + &
rich_mix*(c1 - FRI*FRI)**3
VDC(:,:,k,1) = VDC(:,:,k,2)
endif
else
VISC(:,:,k ) = bckgrnd_vvc(:,:,k,bid)
if ( k < km ) then
VDC (:,:,k,2) = bckgrnd_vdc(:,:,k,bid)
VDC(:,:,k,1) = VDC(:,:,k,2)
endif
endif
endif ! ltidal_mixing
!-----------------------------------------------------------------------
!
! set seafloor values to zero
!
!-----------------------------------------------------------------------
!DIR$ NODEP
!DIR$ COLLAPSE
do j=1,ny_block
!DIR$ NODEP
do i=1,nx_block
if ( k >= KMT(i,j,bid) ) then
VISC(i,j,k ) = c0
VDC (i,j,k,1) = c0
VDC (i,j,k,2) = c0
endif
end do
end do
if (tavg_requested(tavg_KVMIX)) then
! k index shifted because KVMIX is at cell bottom
! while output axis is at cell top
call accumulate_tavg_field(KVMIX,tavg_KVMIX,bid,k)
endif
if (tavg_requested(tavg_KVMIX_M)) then
! k index shifted because KVMIX_M is at cell bottom
! while output axis is at cell top
call accumulate_tavg_field(KVMIX_M,tavg_KVMIX_M,bid,k)
endif
if (tavg_requested(tavg_VDC_BCK)) then
! k index shifted because bckgrnd_vdc is at cell bottom
! while output axis is at cell top
call accumulate_tavg_field(bckgrnd_vdc(:,:,k,bid),tavg_VDC_BCK,bid,k)
endif
if (tavg_requested(tavg_VVC_BCK)) then
! k index shifted because bckgrnd_vdc is at cell bottom
! while output axis is at cell top
call accumulate_tavg_field(bckgrnd_vvc(:,:,k,bid),tavg_VVC_BCK,bid,k)
endif
if (tavg_requested(tavg_TPOWER)) then
WORK1(:,:) = KVMIX(:,:)*RHOMIX(:,:,k)*DBLOC(:,:,k)/ &
(zgrid(k) - zgrid(k+1))
call accumulate_tavg_field(WORK1,tavg_TPOWER,bid,k)
endif
!-----------------------------------------------------------------------
!
! move to next level.
!
!-----------------------------------------------------------------------
end do
!-----------------------------------------------------------------------
!
! fill extra coefficients for blmix
!
!-----------------------------------------------------------------------
VISC(:,:,0 ) = c0
VDC (:,:,0,:) = c0
VISC(:,:,km+1 ) = c0
VDC (:,:,km+1,:) = c0
!-----------------------------------------------------------------------
!EOC
end subroutine ri_iwmix
!***********************************************************************
!BOP
! !IROUTINE: bldepth
! !INTERFACE:
subroutine bldepth (DBLOC, DBSFC, TRCR, UUU, VVV, STF, SHF_QSW, & 2,16
HBLT, USTAR, BFSFC, STABLE, KBL, &
this_block, SMF, SMFT)
! !DESCRIPTION:
! This routine computes the ocean boundary layer depth defined as
! the shallowest depth where the bulk Richardson number is equal to
! a critical value, Ricr.
!
! NOTE: bulk richardson numbers are evaluated by computing
! differences between values at zgrid(kl) $< 0$ and surface
! reference values. currently, the reference values are equal
! to the values in the surface layer. when using higher
! vertical grid resolution, these reference values should be
! computed as the vertical averages from the surface down to
! epssfc*zgrid(kl).
!
! This routine also computes where surface forcing is stable
! or unstable (STABLE)
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TRCR ! tracers at current time
real (r8), dimension(nx_block,ny_block,km), intent(in) :: &
UUU,VVV, &! velocities at current time
DBLOC, &! buoyancy difference between adjacent levels
DBSFC ! buoyancy difference between level and surface
real (r8), dimension(nx_block,ny_block,nt), intent(in) :: &
STF ! surface forcing for all tracers
real (r8), dimension(nx_block,ny_block), intent(in) :: &
SHF_QSW ! short-wave forcing
real (r8), dimension(nx_block,ny_block,2), intent(in), optional :: &
SMF, &! surface momentum forcing at U points
SMFT ! surface momentum forcing at T points
! *** either one or the other (not
! *** both) should be passed
type (block), intent(in) :: &
this_block ! block information for current block
! !OUTPUT PARAMETERS:
integer (int_kind), dimension(nx_block,ny_block), intent(out) :: &
KBL ! index of first lvl below hbl
real (r8), dimension(nx_block,ny_block), intent(out) :: &
HBLT, &! boundary layer depth
BFSFC, &! Bo+radiation absorbed to d
STABLE, &! =1 stable forcing; =0 unstab
USTAR ! surface friction velocity
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
i,j, &! loop indices
bid, &! local block index
kupper, kup, kdn, ktmp, kl ! vertical level indices
real (r8), dimension(nx_block,ny_block) :: &
VSHEAR, &! (velocity shear re sfc)^2
SIGMA, &! d/hbl
WM, WS, &! turb vel scale functions
BO, &! surface buoyancy forcing
BOSOL, &! radiative buoyancy forcing
TALPHA, &! temperature expansion coeff
SBETA, &! salinity expansion coeff
RHO1, &! density at the surface
WORK, &! temp array
ZKL, &! depth at current z level
B_FRQNCY, &! buoyancy frequency
RSH_HBLT, &! resolved shear contribution to HBLT (fraction)
HLANGM, &! Langmuir depth
HEKMAN, &! Eckman depth limit
HLIMIT ! limit to mixed-layer depth
! (= min(HEKMAN,HMONOB))
real (r8), dimension(nx_block,ny_block,3) :: &
RI_BULK, &! Bulk Ri number at 3 lvls
HMONOB ! Monin-Obukhov depth limit
real (r8) :: &
absorb_frac, &! shortwave absorption frac
sqrt_arg, &! dummy sqrt argument
z_upper, z_up ! upper depths for RI_BULK interpolation
real (r8) :: &
a_co, b_co, c_co ! coefficients of the quadratic equation
! $(a_{co}z^2+b_{co}|z|+c_{co}=Ri_b) used to
! find the boundary layer depth. when
! finding the roots, c_co = c_co - Ricr
real (r8) :: &
slope_up ! slope of the above quadratic equation
! at zup. this is used as a boundary
! condition to determine the coefficients.
!-----------------------------------------------------------------------
!
! compute friction velocity USTAR. compute on U-grid and average
! to T-grid.
!
!-----------------------------------------------------------------------
bid = this_block%local_id
if (present(SMFT)) then
USTAR = sqrt(sqrt(SMFT(:,:,1)**2 + SMFT(:,:,2)**2))
else
WORK = sqrt(sqrt(SMF(:,:,1)**2 + SMF(:,:,2)**2))
call ugrid_to_tgrid(USTAR,WORK,bid)
endif
!-----------------------------------------------------------------------
!
! compute density and expansion coefficients at surface
!
!-----------------------------------------------------------------------
WORK = merge(-c2,TRCR(:,:,1,1),TRCR(:,:,1,1) < -c2)
call state(1,1,WORK,TRCR(:,:,1,2),this_block, &
RHOFULL=RHO1, DRHODT=TALPHA, DRHODS=SBETA)
!-----------------------------------------------------------------------
!
! compute turbulent and radiative sfc buoyancy forcing
!
!-----------------------------------------------------------------------
do j=1,ny_block
do i=1,nx_block
if (RHO1(i,j) /= c0) then
BO (i,j) = grav*(-TALPHA(i,j)*STF(i,j,1) - &
SBETA (i,j)*STF(i,j,2))/RHO1(i,j)
BOSOL(i,j) = -grav*TALPHA(i,j)*SHF_QSW(i,j)/RHO1(i,j)
else
BO (i,j) = c0
BOSOL(i,j) = c0
endif
end do
end do
!-----------------------------------------------------------------------
!
! Find bulk Richardson number at every grid level until > Ricr
! max values when Ricr never satisfied are KBL = KMT and
! HBLT = -zgrid(KMT)
!
! NOTE: the reference depth is -epssfc/2.*zgrid(i,k), but the
! reference u,v,t,s values are simply the surface layer
! values and not the averaged values from 0 to 2*ref.depth,
! which is necessary for very fine grids(top layer < 2m
! thickness)
!
!
! Initialize hbl and kbl to bottomed out values
! Initialize HEKMAN and HLIMIT (= HMONOB until reset) to model bottom
! Initialize Monin Obukhov depth to value at z_up
! Set HMONOB=-zgrid(km) if unstable
!
!-----------------------------------------------------------------------
kupper = 1
kup = 2
kdn = 3
z_upper = c0
z_up = zgrid(1)
RI_BULK(:,:,kupper) = c0
RI_BULK(:,:,kup) = c0
KBL = merge(KMT(:,:,bid), 1, (KMT(:,:,bid) > 1))
HLANGM = c0
do kl=1,km
if (partial_bottom_cells) then
if (kl > 1) then
ZKL = -zgrid(kl-1) + p5*(DZT(:,:,kl ,bid) + &
DZT(:,:,kl-1,bid))
else
ZKL = -zgrid(1)
endif
else
ZKL = -zgrid(kl)
endif
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
if (kl == KBL(i,j)) HBLT(i,j) = ZKL(i,j)
end do
end do
enddo
if ( lcheckekmo ) then
HEKMAN = -zgrid(km) + eps
HLIMIT = -zgrid(km) + eps
if ( lshort_wave ) then
select case (sw_absorption_type)
case ('top-layer')
BFSFC = BO + BOSOL
case ('jerlov')
call sw_absorb_frac(-z_up,absorb_frac)
BFSFC = BO + BOSOL * (c1 - absorb_frac)
case ('chlorophyll')
call sw_trans_chl(1,this_block)
BFSFC = BO + BOSOL*(c1-TRANS(:,:,bid))
end select
else
BFSFC = BO
endif
STABLE = merge(c1, c0, BFSFC >= c0)
BFSFC = BFSFC + STABLE*eps
WORK = STABLE * cmonob*USTAR*USTAR*USTAR/vonkar/BFSFC &
+ (STABLE -c1)*zgrid(km)
HMONOB(:,:,kup) = merge( -z_up+eps, WORK, WORK <= -z_up )
endif
RSH_HBLT = c0
!-----------------------------------------------------------------------
!
! compute velocity shear squared on U-grid and use the maximum
! of the four surrounding U-grid values for the T-grid.
!
!-----------------------------------------------------------------------
do kl = 2,km
WORK = (UUU(:,:,1)-UUU(:,:,kl))**2 + &
(VVV(:,:,1)-VVV(:,:,kl))**2
if (partial_bottom_cells) then
WORK = WORK/(-zgrid(kl-1) + &
p5*(DZU(:,:,kl ,bid) + &
DZU(:,:,kl-1,bid) - &
DZU(:,:,1 ,bid)))**2
ZKL = -zgrid(kl-1) + p5*(DZT(:,:,kl ,bid) + &
DZT(:,:,kl-1,bid))
else
ZKL = -zgrid(kl)
endif
VSHEAR(:,1) = c0
do j=2,ny_block
VSHEAR(1,j) = c0
do i=2,nx_block
VSHEAR(i,j) = max(WORK(i,j ), WORK(i-1,j ), &
WORK(i,j-1), WORK(i-1,j-1))
enddo
enddo
!-----------------------------------------------------------------------
!
! compute bfsfc= Bo + radiative contribution down to hbf * hbl
! add epsilon to BFSFC to ensure never = 0
!
!-----------------------------------------------------------------------
if (lshort_wave) then
select case (sw_absorption_type)
case ('top-layer')
BFSFC = BO + BOSOL
case ('jerlov')
do j=1,ny_block
do i=1,nx_block
call sw_absorb_frac(ZKL(i,j), absorb_frac)
BFSFC(i,j) = BO(i,j) + BOSOL(i,j)*(c1 - absorb_frac)
enddo
enddo
case ('chlorophyll')
call sw_trans_chl(2*kl-1,this_block)
BFSFC = BO + BOSOL*(c1-TRANS(:,:,bid))
end select
else
BFSFC = BO
endif
STABLE = MERGE(c1, c0, BFSFC >= c0)
BFSFC = BFSFC + STABLE*eps
!-----------------------------------------------------------------------
!
! compute the Ekman and Monin Obukhov depths using above stability
!
!-----------------------------------------------------------------------
if (lcheckekmo) then
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
if ( STABLE(i,j) > p5 .and. HEKMAN(i,j) >= -zgrid(km) ) then
HEKMAN(i,j) = max(ZKL(i,j), &
cekman*USTAR(i,j)/(abs(FCORT(i,j,bid))+eps))
endif
end do
end do
HMONOB(:,:,kdn) = STABLE*cmonob*USTAR*USTAR*USTAR/vonkar/BFSFC + &
(STABLE-c1)*zgrid(km)
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
if (HMONOB(i,j,kdn) <= ZKL(i,j) .and. &
HMONOB(i,j,kup) > -z_up) then
WORK(i,j) = (HMONOB(i,j,kdn) - HMONOB(i,j,kup))/ &
(z_up + ZKL(i,j))
HLIMIT(i,j) = (HMONOB(i,j,kdn) - WORK(i,j)*ZKL(i,j))/ &
(c1 - WORK(i,j))
endif
end do
end do
endif
!-----------------------------------------------------------------------
!
! compute velocity scales at sigma, for hbl = -zgrid(kl)
!
!-----------------------------------------------------------------------
SIGMA = epssfc
call wscale(SIGMA, ZKL, USTAR, BFSFC, 2, WM, WS)
!-----------------------------------------------------------------------
!
! compute the turbulent shear contribution to RI_BULK and store
! in WM.
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
if (kl < km) then
B_FRQNCY = sqrt( &
p5*(DBLOC(:,:,kl) + abs(DBLOC(:,:,kl)) + eps2)/ &
(p5*(DZT(:,:,kl,bid) + DZT(:,:,kl+1,bid))) )
else
B_FRQNCY = sqrt( &
p5*(DBLOC(:,:,kl) + abs(DBLOC(:,:,kl)) + eps2)/ &
DZT(:,:,kl,bid) )
end if
else
B_FRQNCY = sqrt( &
p5*(DBLOC(:,:,kl) + abs(DBLOC(:,:,kl)) + eps2)/ &
(zgrid(kl)-zgrid(kl+1)) )
endif
WM = ZKL*WS*B_FRQNCY* &
( (Vtc/Ricr(kl))*max(2.1_r8 - 200.0_r8*B_FRQNCY,concv) )
!-----------------------------------------------------------------------
!
! compute bulk Richardson number at new level
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
WORK = merge( DBSFC(:,:,kl)/(-zgrid(kl-1)+ &
p5*(DZT(:,:,kl-1,bid) + &
DZT(:,:,kl ,bid) - &
DZT(:,:,1 ,bid))), &
c0, KMT(:,:,bid) >= kl)
WM = WM/(-zgrid(kl-1) + &
p5*(DZT(:,:,kl-1,bid) + &
DZT(:,:,kl ,bid) - &
DZT(:,:,1 ,bid)))**2
RI_BULK(:,:,kdn) = WORK/(VSHEAR+WM+eps/(-zgrid(kl-1)+ &
p5*(DZU(:,:,kl,bid) + &
DZU(:,:,kl-1,bid) - &
DZU(:,:,1,bid)))**2)
else
WORK = MERGE( (zgrid(1)-zgrid(kl))*DBSFC(:,:,kl), &
c0, KMT(:,:,bid) >= kl)
if ( linertial ) then
RI_BULK(:,:,kdn) = WORK/(VSHEAR+WM+USTAR*BOLUS_SP(:,:,bid)+eps)
else
RI_BULK(:,:,kdn) = WORK/(VSHEAR+WM+eps)
endif
endif
!-----------------------------------------------------------------------
!
! find hbl where Rib = Ricr. if possible, use a quadratic
! interpolation. if not, linearly interpolate. the quadratic
! equation coefficients are determined using the slope and
! Ri_bulk at z_up and Ri_bulk at zgrid(kl). the slope at
! z_up is computed linearly between z_upper and z_up.
!
! compute Langmuir depth always
!-----------------------------------------------------------------------
do j=1,ny_block
do i=1,nx_block
if ( KBL(i,j) == KMT(i,j,bid) .and. &
RI_BULK(i,j,kdn) > Ricr(kl) ) then
slope_up = (RI_BULK(i,j,kupper) - RI_BULK(i,j,kup))/ &
(z_up - z_upper)
a_co = (RI_BULK(i,j,kdn) - RI_BULK(i,j,kup) - &
slope_up*(ZKL(i,j) + z_up) )/(z_up + ZKL(i,j))**2
b_co = slope_up + c2 * a_co * z_up
c_co = RI_BULK(i,j,kup) + &
z_up*(a_co*z_up + slope_up) - Ricr(kl)
sqrt_arg = b_co**2 - c4*a_co*c_co
if ( ( abs(b_co) > eps .and. abs(a_co)/abs(b_co) <= eps ) &
.or. sqrt_arg <= c0 ) then
HBLT(i,j) = -z_up + (z_up + ZKL(i,j)) * &
(Ricr(kl) - RI_BULK(i,j,kup))/ &
(RI_BULK(i,j,kdn) - RI_BULK(i,j,kup))
else
HBLT(i,j) = (-b_co + sqrt(sqrt_arg)) / (c2*a_co)
endif
KBL(i,j) = kl
RSH_HBLT(i,j) = (VSHEAR(i,j)*Ricr(kl)/ &
(DBSFC(i,j,kl)+eps))/HBLT(i,j)
HLANGM(i,j) = USTAR(i,j) * SQRT( FSTOKES(i,j,bid)*ZKL(i,j)/(DBSFC(i,j,kl)+eps) ) &
/ 0.9_r8
endif
enddo
enddo
!-----------------------------------------------------------------------
!
! swap klevel indices and move to next level
!
!-----------------------------------------------------------------------
ktmp = kupper
kupper = kup
kup = kdn
kdn = ktmp
z_upper = z_up
z_up = zgrid(kl)
end do
!-----------------------------------------------------------------------
!
! apply Langmuir parameterization if requested
!
!-----------------------------------------------------------------------
if ( llangmuir ) then
do kl = km,2,-1
where ( HLANGM > HBLT .and. &
HLANGM > -zgrid(kl-1) .and. &
HLANGM <= ZKL )
HBLT = HLANGM
KBL = kl
end where
enddo
endif
!-----------------------------------------------------------------------
!
! first combine Ekman and Monin-Obukhov depth limits. then apply
! these restrictions to HBLT. note that HLIMIT is set to -zgrid(km)
! in unstable forcing.
!
!-----------------------------------------------------------------------
if ( lcheckekmo ) then
where ( HEKMAN < HLIMIT ) HLIMIT = HEKMAN
do kl = 2,km
where ( HLIMIT < HBLT .and. &
HLIMIT > -zgrid(kl-1) .and. &
HLIMIT <= ZKL )
HBLT = HLIMIT
KBL = kl
end where
enddo
endif
!-----------------------------------------------------------------------
!
! apply a Gaussian filter
!
!-----------------------------------------------------------------------
call smooth_hblt (.true., .false., bid, HBLT=HBLT, KBL=KBL)
!-----------------------------------------------------------------------
!
! correct stability and buoyancy forcing for SW up to boundary layer
!
!-----------------------------------------------------------------------
if (lshort_wave) then
select case (sw_absorption_type)
case ('top-layer')
BFSFC = BO + BOSOL
! QSW_HBL = SHF_QSW
if (tavg_requested(tavg_QSW_HBL)) then
WORK = SHF_QSW/hflux_factor
call accumulate_tavg_field(WORK,tavg_QSW_HBL,bid,1)
endif
case ('jerlov')
do j = 1,ny_block
do i = 1,nx_block
call sw_absorb_frac(HBLT(i,j),absorb_frac)
BFSFC(i,j) = BO(i,j) + BOSOL(i,j)*(c1 - absorb_frac)
enddo
enddo
if (tavg_requested(tavg_QSW_HBL)) then
!QSW_HBL(i,j) = SHF_QSW(i,j)*(c1-absorb_frac) ! boundary layer sw
WORK = SHF_QSW*(c1-absorb_frac)/hflux_factor
call accumulate_tavg_field(WORK,tavg_QSW_HBL,bid,1)
endif
case ('chlorophyll')
ZTRANS(:,:,bid) = HBLT(:,:)
call sw_trans_chl(0,this_block)
BFSFC = BO + BOSOL*(c1-TRANS(:,:,bid))
if (tavg_requested(tavg_QSW_HBL)) then
!QSW_HBL = SHF_QSW *(c1-TRANS(:,:,bid)) ! boundary layer sw heating
WORK = SHF_QSW*(c1-TRANS(:,:,bid))/hflux_factor
call accumulate_tavg_field(WORK,tavg_QSW_HBL,bid,1)
endif
end select
endif
STABLE = MERGE(c1, c0, BFSFC >= c0)
BFSFC = BFSFC + STABLE * eps ! ensures bfsfc never=0
!-----------------------------------------------------------------------
!EOC
end subroutine bldepth
!***********************************************************************
!BOP
! !IROUTINE: blmix
! !INTERFACE:
subroutine blmix(VISC, VDC, HBLT, USTAR, BFSFC, STABLE, & 1,3
KBL, GHAT, this_block)
! !DESCRIPTION:
! This routine computes mixing coefficients within boundary layer
! which depend on surface forcing and the magnitude and gradient
! of interior mixing below the boundary layer (matching). These
! quantities have been computed in other routines.
!
! Caution: if mixing bottoms out at hbl = -zgrid(km) then
! fictitious layer at km+1 is needed with small but finite width
! hwide(km+1).
!
! !REVISION HISTORY:
! same as module
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,0:km+1), intent(inout) :: &
VISC ! interior mixing coeff on input
! combined interior/bndy layer coeff output
real (r8), dimension(nx_block,ny_block,0:km+1,2),intent(inout) :: &
VDC ! diffusivity for tracer diffusion
! !INPUT PARAMETERS:
integer (int_kind), dimension(nx_block,ny_block), intent(in) :: &
KBL ! index of first lvl below hbl
real (r8), dimension(nx_block,ny_block), intent(in) :: &
HBLT, & ! boundary layer depth
BFSFC, & ! surface buoyancy forcing
STABLE, & ! =1 stable forcing; =0 unstab
USTAR ! surface friction velocity
type (block), intent(in) :: &
this_block ! block information for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km),intent(out) :: &
GHAT ! non-local mixing coefficient
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k,kp1, &! dummy k level index
i,j, &! horizontal indices
bid ! local block index
integer (int_kind), dimension(nx_block,ny_block) :: &
KN ! klvl closest to HBLT
real (r8), dimension(nx_block,ny_block,km,3) :: &
BLMC ! bndy layer mixing coefs
real (r8), dimension(nx_block,ny_block,3) :: &
GAT1, &! shape function at sigma=1
DAT1, &! derivative of shape function
DKM1 ! bndy layer difs at kbl-1 lvl
real (r8), dimension(nx_block,ny_block) :: &
WM,WS, &! turbulent velocity scales
CASEA, &! =1 in case A, =0 in case B
SIGMA, &! normalized depth (d/hbl)
VISCH, &! viscosity at hbl
DIFTH, &! temp diffusivity at hbl
DIFSH, &! tracer diffusivity at hbl
DELHAT, R, DVDZUP, DVDZDN, &
VISCP, DIFTP, DIFSP, F1, &
WORK1,WORK2
!-----------------------------------------------------------------------
!
! compute velocity scales at hbl
!
!-----------------------------------------------------------------------
bid = this_block%local_id
SIGMA = epssfc
call wscale(SIGMA, HBLT, USTAR, BFSFC, 3, WM, WS)
!-----------------------------------------------------------------------
!
! determine caseA = 0 if closer to KBL than KBL-1
! KN is then the closest klevel to HBLT
!
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
k = KBL(i,j)
if (k == 1) then
CASEA(i,j) = p5 + SIGN(p5, -zgrid(0)-HBLT(i,j))
else
CASEA(i,j) = p5 + SIGN(p5, &
-zgrid(k-1)+p5*DZT(i,j,k-1,bid)-HBLT(i,j))
endif
enddo
enddo
else
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
k = KBL(i,j)
CASEA(i,j) = p5 + SIGN(p5, -zgrid(k)-p5*hwide(k)-HBLT(i,j))
enddo
enddo
endif
KN = NINT(CASEA)*(KBL-1) + (1-NINT(CASEA))*KBL
!-----------------------------------------------------------------------
!
! find the interior viscosities and derivatives at hbl by
! interpolating derivative values at vertical interfaces. compute
! matching conditions for shape function.
!
!-----------------------------------------------------------------------
F1 = STABLE*c5*BFSFC/(USTAR**4+eps)
do k=1,km
if (partial_bottom_cells) then
if (k == 1) then
WORK1 = c0
else
WORK1 = DZT(:,:,k-1,bid)
end if
if (k == km) then
WORK2 = eps
else
WORK2 = DZT(:,:,k+1,bid)
end if
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
if (k == KN(i,j)) then
DELHAT(i,j) = - zgrid(k-1) + DZT(i,j,k,bid) + &
p5*WORK1(i,j) - HBLT(i,j)
R (i,j) = c1 - DELHAT(i,j) /DZT(i,j,k,bid)
DVDZUP(i,j) = (VISC(i,j,k-1) - VISC(i,j,k ))/DZT(i,j,k,bid)
DVDZDN(i,j) = (VISC(i,j,k ) - VISC(i,j,k+1))/WORK2(i,j)
VISCP (i,j) = p5*( (c1-R(i,j))* &
(DVDZUP(i,j) + abs(DVDZUP(i,j))) + &
R(i,j) * &
(DVDZDN(i,j) + abs(DVDZDN(i,j))) )
DVDZUP(i,j) = (VDC(i,j,k-1,2) - VDC(i,j,k ,2))/DZT(i,j,k,bid)
DVDZDN(i,j) = (VDC(i,j,k ,2) - VDC(i,j,k+1,2))/WORK2(i,j)
DIFSP (i,j) = p5*( (c1-R(i,j))* &
(DVDZUP(i,j) + abs(DVDZUP(i,j))) + &
R(i,j) * &
(DVDZDN(i,j) + abs(DVDZDN(i,j))) )
DVDZUP(i,j) = (VDC(i,j,k-1,1) - VDC(i,j,k ,1))/DZT(i,j,k,bid)
DVDZDN(i,j) = (VDC(i,j,k ,1) - VDC(i,j,k+1,1))/WORK2(i,j)
DIFTP (i,j) = p5*( (c1-R(i,j))* &
(DVDZUP(i,j) + abs(DVDZUP(i,j))) + &
R(i,j) * &
(DVDZDN(i,j) + abs(DVDZDN(i,j))) )
VISCH(i,j) = VISC(i,j,k) + VISCP(i,j)*DELHAT(i,j)
DIFSH(i,j) = VDC(i,j,k,2) + DIFSP(i,j)*DELHAT(i,j)
DIFTH(i,j) = VDC(i,j,k,1) + DIFTP(i,j)*DELHAT(i,j)
GAT1(i,j,1) = VISCH(i,j) / HBLT(i,j) /(WM(i,j)+eps)
DAT1(i,j,1) = -VISCP(i,j)/(WM(i,j)+eps) + F1(i,j)*VISCH(i,j)
GAT1(i,j,2) = DIFSH(i,j) / HBLT(i,j) /(WS(i,j)+eps)
DAT1(i,j,2) = -DIFSP(i,j)/(WS(i,j)+eps) + F1(i,j)*DIFSH(i,j)
GAT1(i,j,3) = DIFTH(i,j) / HBLT(i,j) /(WS(i,j)+eps)
DAT1(i,j,3) = -DIFTP(i,j)/(WS(i,j)+eps) + F1(i,j)*DIFTH(i,j)
endif
end do
end do
else
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
if (k == KN(i,j)) then
DELHAT(i,j) = p5*hwide(k) - zgrid(k) - HBLT(i,j)
R (i,j) = c1 - DELHAT(i,j) / hwide(k)
DVDZUP(i,j) = (VISC(i,j,k-1) - VISC(i,j,k ))/hwide(k)
DVDZDN(i,j) = (VISC(i,j,k ) - VISC(i,j,k+1))/hwide(k+1)
VISCP (i,j) = p5*( (c1-R(i,j))* &
(DVDZUP(i,j) + abs(DVDZUP(i,j))) + &
R(i,j) * &
(DVDZDN(i,j) + abs(DVDZDN(i,j))) )
DVDZUP(i,j) = (VDC(i,j,k-1,2) - VDC(i,j,k ,2))/hwide(k)
DVDZDN(i,j) = (VDC(i,j,k ,2) - VDC(i,j,k+1,2))/hwide(k+1)
DIFSP (i,j) = p5*( (c1-R(i,j))* &
(DVDZUP(i,j) + abs(DVDZUP(i,j))) + &
R(i,j) * &
(DVDZDN(i,j) + abs(DVDZDN(i,j))) )
DVDZUP(i,j) = (VDC(i,j,k-1,1) - VDC(i,j,k ,1))/hwide(k)
DVDZDN(i,j) = (VDC(i,j,k ,1) - VDC(i,j,k+1,1))/hwide(k+1)
DIFTP (i,j) = p5*( (c1-R(i,j))* &
(DVDZUP(i,j) + abs(DVDZUP(i,j))) + &
R(i,j) * &
(DVDZDN(i,j) + abs(DVDZDN(i,j))) )
VISCH(i,j) = VISC(i,j,k) + VISCP(i,j)*DELHAT(i,j)
DIFSH(i,j) = VDC(i,j,k,2) + DIFSP(i,j)*DELHAT(i,j)
DIFTH(i,j) = VDC(i,j,k,1) + DIFTP(i,j)*DELHAT(i,j)
GAT1(i,j,1) = VISCH(i,j) / HBLT(i,j) /(WM(i,j)+eps)
DAT1(i,j,1) = -VISCP(i,j)/(WM(i,j)+eps) + &
F1(i,j)*VISCH(i,j)
GAT1(i,j,2) = DIFSH(i,j) / HBLT(i,j) /(WS(i,j)+eps)
DAT1(i,j,2) = -DIFSP(i,j)/(WS(i,j)+eps) + &
F1(i,j)*DIFSH(i,j)
GAT1(i,j,3) = DIFTH(i,j) / HBLT(i,j) /(WS(i,j)+eps)
DAT1(i,j,3) = -DIFTP(i,j)/(WS(i,j)+eps) + &
F1(i,j)*DIFTH(i,j)
endif
end do
end do
endif ! pbc
enddo
DAT1 = min(DAT1,c0)
!-----------------------------------------------------------------------
!
! compute the dimensionless shape functions and diffusivities
! at the grid interfaces. also compute function for non-local
! transport term (GHAT).
!
!-----------------------------------------------------------------------
do k = 1,km
if (partial_bottom_cells) then
if (k > 1) then
SIGMA = (-zgrid(k-1) + p5*DZT(:,:,k-1,bid) + &
DZT(:,:,k,bid)) / HBLT
else
SIGMA = (-zgrid(k) + p5*hwide(k)) / HBLT
end if
else
SIGMA = (-zgrid(k) + p5*hwide(k)) / HBLT
endif
F1 = min(SIGMA,epssfc)
call wscale(F1, HBLT, USTAR, BFSFC, 3, WM, WS)
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
BLMC(i,j,k,1) = HBLT(i,j)*WM(i,j)*SIGMA(i,j)* &
(c1 + SIGMA(i,j)*((SIGMA(i,j)-c2) + &
(c3-c2*SIGMA(i,j))*GAT1(i,j,1) + &
(SIGMA(i,j)-c1)*DAT1(i,j,1)))
BLMC(i,j,k,2) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* &
(c1 + SIGMA(i,j)*((SIGMA(i,j)-c2) + &
(c3-c2*SIGMA(i,j))*GAT1(i,j,2) + &
(SIGMA(i,j)-c1)*DAT1(i,j,2)))
BLMC(i,j,k,3) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* &
(c1 + SIGMA(i,j)*((SIGMA(i,j)-c2) + &
(c3-c2*SIGMA(i,j))*GAT1(i,j,3) + &
(SIGMA(i,j)-c1)*DAT1(i,j,3)))
GHAT(i,j,k) = (c1-STABLE(i,j))* cg/(WS(i,j)*HBLT(i,j) +eps)
end do
end do
end do
!-----------------------------------------------------------------------
!
! find diffusivities at kbl-1 grid level
!
!-----------------------------------------------------------------------
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
k = KBL(i,j) - 1
SIGMA(i,j) = -zgrid(k)/HBLT(i,j)
enddo
enddo
F1 = min(SIGMA,epssfc)
call wscale(F1, HBLT, USTAR, BFSFC, 3, WM, WS)
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
DKM1(i,j,1) = HBLT(i,j)*WM(i,j)*SIGMA(i,j)* &
(c1+SIGMA(i,j)*((SIGMA(i,j)-c2) + &
(c3-c2*SIGMA(i,j))*GAT1(i,j,1) + &
(SIGMA(i,j)-c1)*DAT1(i,j,1)))
DKM1(i,j,2) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* &
(c1+SIGMA(i,j)*((SIGMA(i,j)-c2) + &
(c3-c2*SIGMA(i,j))*GAT1(i,j,2) + &
(SIGMA(i,j)-c1)*DAT1(i,j,2)))
DKM1(i,j,3) = HBLT(i,j)*WS(i,j)*SIGMA(i,j)* &
(c1+SIGMA(i,j)*((SIGMA(i,j)-c2) + &
(c3-c2*SIGMA(i,j))*GAT1(i,j,3) + &
(SIGMA(i,j)-c1)*DAT1(i,j,3)))
end do
end do
!-----------------------------------------------------------------------
!
! compute the enhanced mixing
!
!-----------------------------------------------------------------------
!DIR$ NOVECTOR
do k=1,km-1
if (partial_bottom_cells) then
if (k == 1) then
WORK1 = -p5*DZT(:,:,k,bid)
else
WORK1 = zgrid(k-1) - p5*(DZT(:,:,k-1,bid) + &
DZT(:,:,k ,bid))
end if
where (k == (KBL - 1)) &
DELHAT = (HBLT + WORK1)/(p5*(DZT(:,:,k ,bid) + &
DZT(:,:,k+1,bid)))
else
where (k == (KBL - 1)) &
DELHAT = (HBLT + zgrid(k))/(zgrid(k)-zgrid(k+1))
endif
!DIR$ COLLAPSE
do j=1,ny_block
do i=1,nx_block
if (k == (KBL(i,j) - 1)) then
BLMC(i,j,k,1) = (c1-DELHAT(i,j))*VISC(i,j,k) + &
DELHAT(i,j) *( &
(c1-DELHAT(i,j))**2 *DKM1(i,j,1) + &
DELHAT(i,j)**2 *(CASEA(i,j)*VISC(i,j,k) + &
(c1-CASEA(i,j))*BLMC(i,j,k,1)))
BLMC(i,j,k,2) = (c1-DELHAT(i,j))*VDC(i,j,k,2) + &
DELHAT(i,j)*( &
(c1-DELHAT(i,j))**2 *DKM1(i,j,2) + &
DELHAT(i,j)**2 *(CASEA(i,j)*VDC(i,j,k,2) + &
(c1-CASEA(i,j))*BLMC(i,j,k,2)))
BLMC(i,j,k,3) = (c1-DELHAT(i,j))*VDC(i,j,k,1) + &
DELHAT(i,j) *( &
(c1-DELHAT(i,j))**2 *DKM1(i,j,3) + &
DELHAT(i,j)**2 *(CASEA(i,j)*VDC(i,j,k,1) + &
(c1-CASEA(i,j))*BLMC(i,j,k,3)))
GHAT(i,j,k) = (c1-CASEA(i,j)) * GHAT(i,j,k)
endif
end do
end do
end do
!-----------------------------------------------------------------------
!
! combine interior and boundary layer coefficients and nonlocal term
!
!-----------------------------------------------------------------------
!DIR$ NOVECTOR
do k=1,km
!DIR$ NODEP
do j=1,ny_block
!DIR$ NODEP
do i=1,nx_block
if (k < KBL(i,j)) then
VISC(i,j,k) = BLMC(i,j,k,1)
VDC(i,j,k,2) = BLMC(i,j,k,2)
VDC(i,j,k,1) = BLMC(i,j,k,3)
else
GHAT(i,j,k) = c0
endif
end do
end do
enddo
!-----------------------------------------------------------------------
!EOC
end subroutine blmix
!***********************************************************************
!BOP
! !IROUTINE: wscale
! !INTERFACE:
subroutine wscale(SIGMA, HBL, USTAR, BFSFC, m_or_s, WM, WS) 4
! !DESCRIPTION:
! Computes turbulent velocity scales.
!
! For $\zeta \geq 0,
! w_m = w_s = \kappa U^\star/(1+5\zeta)$
!
! For $\zeta_m \leq \zeta < 0,
! w_m = \kappa U^\star (1-16\zeta)^{1\over 4}$
!
! For $\zeta_s \leq \zeta < 0,
! w_s = \kappa U^\star (1-16\zeta)^{1\over 2}$
!
! For $\zeta < \zeta_m,
! w_m = \kappa U^\star (a_m - c_m\zeta)^{1\over 3}$
!
! For $\zeta < \zeta_s,
! w_s = \kappa U^\star (a_s - c_s\zeta)^{1\over 3}$
!
! where $\kappa$ is the von Karman constant.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
m_or_s ! flag =1 for wm only, 2 for ws, 3 for both
real (r8), dimension(nx_block,ny_block), intent(in) :: &
SIGMA, &! normalized depth (d/hbl)
HBL, &! boundary layer depth
BFSFC, &! surface buoyancy forcing
USTAR ! surface friction velocity
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(out) :: &
WM, &! turb velocity scales: momentum
WS ! turb velocity scales: tracer
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: i,j ! dummy loop indices
real (r8), dimension(nx_block,ny_block) :: &
ZETA, &! d/L or sigma*hbl/L(monin-obk)
ZETAH ! sigma*hbl*vonkar*BFSFC or ZETA = ZETAH/USTAR**3
!-----------------------------------------------------------------------
!
! compute zetah and zeta - surface layer is special case
!
!-----------------------------------------------------------------------
ZETAH = SIGMA*HBL*vonkar*BFSFC
ZETA = ZETAH/(USTAR**3 + eps)
!-----------------------------------------------------------------------
!
! compute velocity scales for momentum
!
!-----------------------------------------------------------------------
if (m_or_s == 1 .or. m_or_s == 3) then
do j=1,ny_block
do i=1,nx_block
if (ZETA(i,j) >= c0) then ! stable region
WM(i,j) = vonkar*USTAR(i,j)/(c1 + c5*ZETA(i,j))
else if (ZETA(i,j) >= zeta_m) then
WM(i,j) = vonkar*USTAR(i,j)*(c1 - c16*ZETA(i,j))**p25
else
WM(i,j) = vonkar*(a_m*(USTAR(i,j)**3)-c_m*ZETAH(i,j))**p33
endif
end do
end do
endif
!-----------------------------------------------------------------------
!
! compute velocity scales for tracers
!
!-----------------------------------------------------------------------
if (m_or_s == 2 .or. m_or_s == 3) then
do j=1,ny_block
do i=1,nx_block
if (ZETA(i,j) >= c0) then
WS(i,j) = vonkar*USTAR(i,j)/(c1 + c5*ZETA(i,j))
else if (ZETA(i,j) >= zeta_s) then
WS(i,j) = vonkar*USTAR(i,j)*SQRT(c1 - c16*ZETA(i,j))
else
WS(i,j) = vonkar*(a_s*(USTAR(i,j)**3)-c_s*ZETAH(i,j))**p33
endif
end do
end do
endif
!-----------------------------------------------------------------------
!EOC
end subroutine wscale
!***********************************************************************
!BOP
! !IROUTINE: ddmix
! !INTERFACE:
subroutine ddmix(VDC, TRCR, this_block) 1,2
! !DESCRIPTION:
! $R_\rho$ dependent interior flux parameterization.
! Add double-diffusion diffusivities to Ri-mix values at blending
! interface and below.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TRCR ! tracers at current time
type (block), intent(in) :: &
this_block ! block information for current block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,0:km+1,2),intent(inout) :: &
VDC ! diffusivity for tracer diffusion
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: k,kup,knxt
real (r8), dimension(nx_block,ny_block) :: &
ALPHADT, &! alpha*DT across interfaces
BETADS, &! beta *DS across interfaces
RRHO, &! dd density ratio
DIFFDD, &! dd diffusivity scale
PRANDTL ! prandtl number
real (r8), dimension(nx_block,ny_block,2) :: &
TALPHA, &! temperature expansion coeff
SBETA ! salinity expansion coeff
!-----------------------------------------------------------------------
!
! compute alpha*DT and beta*DS at interfaces. use RRHO and
! PRANDTL for temporary storage for call to state
!
!-----------------------------------------------------------------------
kup = 1
knxt = 2
PRANDTL = merge(-c2,TRCR(:,:,1,1),TRCR(:,:,1,1) < -c2)
call state(1, 1, PRANDTL, TRCR(:,:,1,2), this_block, &
RHOFULL=RRHO, &
DRHODT=TALPHA(:,:,kup), DRHODS=SBETA(:,:,kup))
do k=1,km
if ( k < km ) then
PRANDTL = merge(-c2,TRCR(:,:,k+1,1),TRCR(:,:,k+1,1) < -c2)
call state(k+1, k+1, PRANDTL, TRCR(:,:,k+1,2), &
this_block, &
RHOFULL=RRHO, DRHODT=TALPHA(:,:,knxt), &
DRHODS= SBETA(:,:,knxt))
ALPHADT = -p5*(TALPHA(:,:,kup) + TALPHA(:,:,knxt)) &
*(TRCR(:,:,k,1) - TRCR(:,:,k+1,1))
BETADS = p5*( SBETA(:,:,kup) + SBETA(:,:,knxt)) &
*(TRCR(:,:,k,2) - TRCR(:,:,k+1,2))
kup = knxt
knxt = 3 - kup
else
ALPHADT = c0
BETADS = c0
endif
!-----------------------------------------------------------------------
!
! salt fingering case
!
!-----------------------------------------------------------------------
where ( ALPHADT > BETADS .and. BETADS > c0 )
RRHO = MIN(ALPHADT/BETADS, Rrho0)
DIFFDD = dsfmax*(c1-(RRHO-c1)/(Rrho0-c1))**3
VDC(:,:,k,1) = VDC(:,:,k,1) + 0.7_r8*DIFFDD
VDC(:,:,k,2) = VDC(:,:,k,2) + DIFFDD
endwhere
!-----------------------------------------------------------------------
!
! diffusive convection
!
!-----------------------------------------------------------------------
where ( ALPHADT < c0 .and. BETADS < c0 .and. ALPHADT > BETADS )
RRHO = ALPHADT / BETADS
DIFFDD = 1.5e-2_r8*0.909_r8* &
exp(4.6_r8*exp(-0.54_r8*(c1/RRHO-c1)))
PRANDTL = 0.15_r8*RRHO
elsewhere
RRHO = c0
DIFFDD = c0
PRANDTL = c0
endwhere
where (RRHO > p5) PRANDTL = (1.85_r8 - 0.85_r8/RRHO)*RRHO
VDC(:,:,k,1) = VDC(:,:,k,1) + DIFFDD
VDC(:,:,k,2) = VDC(:,:,k,2) + PRANDTL*DIFFDD
end do
!-----------------------------------------------------------------------
!EOC
end subroutine ddmix
!***********************************************************************
!BOP
! !IROUTINE: buoydiff
! !INTERFACE:
subroutine buoydiff(DBLOC, DBSFC, TRCR, this_block) 1,3
! !DESCRIPTION:
! This routine calculates the buoyancy differences at model levels.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km,nt), intent(in) :: &
TRCR ! tracers at current time
type (block), intent(in) :: &
this_block ! block information for current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,km), intent(out) :: &
DBLOC, &! buoyancy difference between adjacent levels
DBSFC ! buoyancy difference between level and surface
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
k, &! vertical level index
i,j, &! horizontal indices
kprev, klvl, ktmp, &! indices for 2-level TEMPK array
bid ! local block index
real (r8), dimension(nx_block,ny_block) :: &
RHO1, &! density of sfc t,s displaced to k
RHOKM, &! density of t(k-1),s(k-1) displaced to k
RHOK, &! density at level k
TEMPSFC, &! adjusted temperature at surface
TALPHA, &! temperature expansion coefficient
SBETA ! salinity expansion coefficient
real (r8), dimension(nx_block,ny_block,2) :: &
TEMPK ! temp adjusted for freeze at levels k,k-1
!-----------------------------------------------------------------------
!
! calculate density and buoyancy differences at surface
!
!-----------------------------------------------------------------------
TEMPSFC = merge(-c2,TRCR(:,:,1,1),TRCR(:,:,1,1) < -c2)
bid = this_block%local_id
klvl = 2
kprev = 1
TEMPK(:,:,kprev) = TEMPSFC
DBSFC(:,:,1) = c0
!-----------------------------------------------------------------------
!
! calculate DBLOC and DBSFC for all other levels
!
!-----------------------------------------------------------------------
do k = 2,km
TEMPK(:,:,klvl) = merge(-c2,TRCR(:,:,k,1),TRCR(:,:,k,1) < -c2)
call state(k, k, TEMPSFC, TRCR(:,:,1 ,2), &
this_block, RHOFULL=RHO1)
call state(k, k, TEMPK(:,:,kprev), TRCR(:,:,k-1,2), &
this_block, RHOFULL=RHOKM)
call state(k, k, TEMPK(:,:,klvl), TRCR(:,:,k ,2), &
this_block, RHOFULL=RHOK)
do j=1,ny_block
do i=1,nx_block
if (RHOK(i,j) /= c0) then
DBSFC(i,j,k) = grav*(c1 - RHO1 (i,j)/RHOK(i,j))
DBLOC(i,j,k-1) = grav*(c1 - RHOKM(i,j)/RHOK(i,j))
else
DBSFC(i,j,k) = c0
DBLOC(i,j,k-1) = c0
endif
if (k-1 >= KMT(i,j,bid)) DBLOC(i,j,k-1) = c0
end do
end do
ktmp = klvl
klvl = kprev
kprev = ktmp
enddo
DBLOC(:,:,km) = c0
!-----------------------------------------------------------------------
!EOC
end subroutine buoydiff
!***********************************************************************
!BOP
! !IROUTINE: add_kpp_sources
! !INTERFACE:
subroutine add_kpp_sources(SRCARRAY, k, this_block) 1,2
! !DESCRIPTION:
! This routine adds KPP non local mixing term to the tracer source
! tendencies.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k ! vertical level index
type (block), intent(in) :: &
this_block ! block information for current block
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,nt), intent(inout) :: &
SRCARRAY ! array of tracer sources
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (int_kind) :: &
n, &! local dummy index for tracer
bid ! local block address
!-----------------------------------------------------------------------
!
! if KPP not chosen, return
!
!-----------------------------------------------------------------------
if (.not. allocated(KPP_SRC)) return
!-----------------------------------------------------------------------
!
! determine location and add sources
!
!-----------------------------------------------------------------------
bid = this_block%local_id
SRCARRAY = SRCARRAY + KPP_SRC(:,:,k,:,bid)
if (mix_pass /= 1) then
do n=1,nt
if (tavg_requested(tavg_KPP_SRC(n))) &
call accumulate_tavg_field(KPP_SRC(:,:,k,n,bid),tavg_KPP_SRC(n),bid,k)
enddo
endif
!-----------------------------------------------------------------------
!EOC
end subroutine add_kpp_sources
!***********************************************************************
!BOP
! !IROUTINE: smooth_hblt
! !INTERFACE:
subroutine smooth_hblt (overwrite_hblt, use_hmxl, & 2,4
bid, HBLT, KBL, SMOOTH_OUT)
! !DESCRIPTION:
! This subroutine uses a 1-1-4-1-1 Laplacian filter one time
! on HBLT or HMXL to reduce any horizontal two-grid-point noise.
! If HBLT is overwritten, KBL is adjusted after smoothing.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
logical (log_kind), intent(in) :: &
overwrite_hblt, & ! if .true., HBLT is overwritten
! if .false., the result is returned in
! a dummy array
use_hmxl ! if .true., smooth HMXL
! if .false., smooth HBLT
integer (int_kind), intent(in) :: &
bid ! local block address
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), optional, intent(inout) :: &
HBLT ! boundary layer depth
integer (int_kind), dimension(nx_block,ny_block), optional, intent(inout) :: &
KBL ! index of first lvl below hbl
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), optional, intent(out) :: &
SMOOTH_OUT ! optional output array containing the
! smoothened field if overwrite_hblt is false
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
character (char_len) :: &
message
integer (int_kind) :: &
i, j, & ! horizontal loop indices
k ! vertical level index
real (r8), dimension(nx_block,ny_block) :: &
WORK1, WORK2
real (r8) :: &
cc, cw, ce, cn, cs, & ! averaging weights
ztmp ! temp for level depth
!-----------------------------------------------------------------------
!
! consistency checks
!
!-----------------------------------------------------------------------
if ( overwrite_hblt .and. ( .not.present(KBL) .or. &
.not.present(HBLT) ) ) then
message = 'incorrect subroutine arguments for smooth_hblt, error # 1'
call exit_POP (sigAbort, trim(message))
endif
if ( .not.overwrite_hblt .and. .not.present(SMOOTH_OUT) ) then
message = 'incorrect subroutine arguments for smooth_hblt, error # 2'
call exit_POP (sigAbort, trim(message))
endif
if ( use_hmxl .and. .not.present(SMOOTH_OUT) ) then
message = 'incorrect subroutine arguments for smooth_hblt, error # 3'
call exit_POP (sigAbort, trim(message))
endif
if ( overwrite_hblt .and. use_hmxl ) then
message = 'incorrect subroutine arguments for smooth_hblt, error # 4'
call exit_POP (sigAbort, trim(message))
endif
!-----------------------------------------------------------------------
!
! perform one smoothing pass since we cannot do the necessary
! boundary updates for multiple passes.
!
!-----------------------------------------------------------------------
if ( use_hmxl ) then
WORK2 = HMXL(:,:,bid)
else
WORK2 = HBLT
endif
WORK1 = WORK2
do j=2,ny_block-1
do i=2,nx_block-1
if ( KMT(i,j,bid) /= 0 ) then
cw = p125
ce = p125
cn = p125
cs = p125
cc = p5
if ( KMT(i-1,j,bid) == 0 ) then
cc = cc + cw
cw = c0
endif
if ( KMT(i+1,j,bid) == 0 ) then
cc = cc + ce
ce = c0
endif
if ( KMT(i,j-1,bid) == 0 ) then
cc = cc + cs
cs = c0
endif
if ( KMT(i,j+1,bid) == 0 ) then
cc = cc + cn
cn = c0
endif
WORK2(i,j) = cw * WORK1(i-1,j) &
+ ce * WORK1(i+1,j) &
+ cs * WORK1(i,j-1) &
+ cn * WORK1(i,j+1) &
+ cc * WORK1(i,j)
endif
enddo
enddo
do k=1,km
do j=2,ny_block-1
do i=2,nx_block-1
if (partial_bottom_cells) then
ztmp = -zgrid(k-1) + p5*(DZT(i,j,k-1,bid) + &
DZT(i,j,k ,bid))
else
ztmp = -zgrid(k)
endif
if ( k == KMT(i,j,bid) .and. WORK2(i,j) > ztmp ) then
WORK2(i,j) = ztmp
endif
enddo
enddo
enddo
if ( overwrite_hblt .and. .not.use_hmxl ) then
HBLT = WORK2
do k=1,km
do j=2,ny_block-1
do i=2,nx_block-1
if (partial_bottom_cells) then
ztmp = -zgrid(k-1) + p5*(DZT(i,j,k-1,bid) + &
DZT(i,j,k ,bid))
else
ztmp = -zgrid(k)
endif
if ( KMT(i,j,bid) /= 0 .and. &
( HBLT(i,j) > -zgrid(k-1) ) .and. &
( HBLT(i,j) <= ztmp ) ) KBL(i,j) = k
enddo
enddo
enddo
else
SMOOTH_OUT = WORK2
endif
!-----------------------------------------------------------------------
end subroutine smooth_hblt
!***********************************************************************
end module vmix_kpp
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||