@ -20,6 +20,7 @@
! KB use get_field , only : get_3d_field
use variables_2d , only : fwf_int
use variables_3d , only : rk , S , hn , kmin , sss
use variables_3d , only : saltbins , saltbinmin , saltbinmax , Si_s
use halo_zones , only : update_3d_halo , wait_halo , D_TAG , H_TAG
IMPLICIT NONE
@ -61,6 +62,7 @@
integer :: salt_check = 0
REALTYPE :: min_salt = _ ZERO_ , max_salt = 40 * _ ONE_
integer , public :: nonnegsalt_method = 0
integer , allocatable :: binidx_kij ( : , : , : )
!
! ! REVISION HISTORY :
! Original author ( s ) : Karsten Bolding & Hans Burchard
@ -77,7 +79,7 @@
! \ label { sec - init - salinity }
!
! ! INTERFACE :
subroutine init_salinity ( hotstart )
subroutine init_salinity ( runtype , hotstart )
!
! ! DESCRIPTION :
!
@ -104,10 +106,11 @@
IMPLICIT NONE
!
! ! INPUT PARAMETERS :
integer , intent ( in ) :: runtype
logical , intent ( in ) :: hotstart
!
! ! LOCAL VARIABLES :
integer :: i , j , k , n , rc
integer :: i , j , k , l , n , rc
integer :: status
NAMELIST / salt / &
salt_method , salt_const , salt_file , &
@ -115,6 +118,7 @@
salt_adv_split , salt_adv_hor , salt_adv_ver , &
avmols , salt_AH_method , salt_AH_const , salt_AH_Prt , &
salt_AH_stirr_const , sss_nudging_time , &
saltbins , saltbinmin , saltbinmax , &
salt_check , min_salt , max_salt , &
nonnegsalt_method
! EOP
@ -136,10 +140,6 @@
LEVEL3 'avmols = ' , real ( avmols )
end if
if ( . not . hotstart ) then
call init_salinity_field ( )
end if
! Sanity checks for advection specifications
LEVEL3 'Advection of salinity'
if ( salt_adv_hor . eq . J7 ) stop 'init_salinity: J7 not implemented yet'
@ -188,6 +188,22 @@
allocate ( sss ( E2DFIELD ) , stat = rc )
end if
if ( runtype . eq . 4 . and . saltbins . gt . 0 ) then
LEVEL3 'salinity space for TEF:' , real ( saltbinmin ) , real ( saltbinmax ) , '(' , saltbins , 'bins)'
allocate ( Si_s ( 0 : saltbins ) , stat = rc )
if ( rc / = 0 ) stop 'postinit_3d: Error allocating memory (Si_s)'
do l = 0 , saltbins
Si_s ( l ) = saltbinmin + ( saltbinmax - saltbinmin ) / saltbins * l
end do
allocate ( binidx_kij ( _ KRANGE_ , I2DFIELD ) , stat = rc )
if ( rc / = 0 ) stop 'postinit_3d: Error allocating memory (binidx_kij)'
binidx_kij = - 1
end if
if ( . not . hotstart ) then
call init_salinity_field ( )
end if
LEVEL3 'salt_check=' , salt_check
if ( salt_check . ne . 0 ) then
LEVEL4 'doing sanity check on salinity'
@ -298,6 +314,14 @@
call wait_halo ( D_TAG )
call mirror_bdy_3d ( S , D_TAG )
if ( allocated ( binidx_kij ) ) then
do j = jmin - HALO , jmax + HALO
do i = imin - HALO , imax + HALO
if ( az ( i , j ) . ne . 0 ) call rebin ( kmax , S ( i , j , : ) , saltbins , Si_s , binidx_kij ( : , i , j ) )
end do
end do
end if
# ifdef DEBUG
write ( debug , * ) 'Leaving init_salinity_field()'
write ( debug , * )
@ -312,7 +336,7 @@
! ! IROUTINE : do_salinity - salinity equation \ label { sec - do - salinity }
!
! ! INTERFACE :
subroutine do_salinity ( n )
subroutine do_salinity ( loop )
!
! ! DESCRIPTION :
!
@ -333,24 +357,40 @@
!
! ! USES :
use advection_3d , only : do_advection_3d
use advection_3d , only : ffluxu2_3d , ffluxv2_3d
use variables_3d , only : dt , cnpar , hn , ho , nuh , uu , vv , ww , hun , hvn
use domain , only : imin , imax , jmin , jmax , kmax , az
use getm_timers , only : tic , toc , TIM_SALT , TIM_SALTH , TIM_MIXANALYSIS
use domain , only : imin , imax , jmin , jmax , kmax , az , au , av
use domain , only : dyu , dxv
use domain , only : nriverl , ri , rj
use getm_timers , only : tic , toc , TIM_SALT , TIM_SALTH , TIM_MIXANALYSIS , TIM_TEF
use variables_3d , only : Sfluxu , Sfluxv , Sfluxw
use variables_3d , only : nummix_S , phymix_S
use variables_3d , only : Sfluxu2 , Sfluxv2
use variables_3d , only : counts_s , flags_s , h_s , hS_s , hS2_s , hpmS_s , hnmS_s
use variables_3d , only : hu_s , uu_s , Sfluxu_s , S2fluxu_s
use variables_3d , only : hv_s , vv_s , Sfluxv_s , S2fluxv_s
use variables_3d , only : fwf_s , fwfS2_s , phymix_S_fwf_s , phymix_S_riv_s
use variables_3d , only : nummix_S , phymix_S , phymix_S_fwf
! $ use omp_lib
IMPLICIT NONE
!
! ! INPUT PARAMETERS :
integer , intent ( in ) :: n
integer , intent ( in ) :: loop
!
! ! LOCAL VARIABLES :
integer :: i , j , k , rc
integer :: i , j , k , l , n , rc
REALTYPE , POINTER :: Res ( : )
REALTYPE , POINTER :: auxn ( : ) , auxo ( : )
REALTYPE , POINTER :: a1 ( : ) , a2 ( : ) , a3 ( : ) , a4 ( : )
REALTYPE , pointer :: fluxw ( : )
REALTYPE :: sss_old ( I2DFIELD )
REALTYPE , dimension ( I3DFIELD ) :: Sold , Sfluxu_tmp , Sfluxv_tmp
REALTYPE , dimension ( I2DFIELD ) :: wrk2d
REALTYPE :: rivmix ( 0 : kmax , 1 : nriverl )
REALTYPE :: deltaS , hu_k , uu_k
REALTYPE :: S_k ( 0 : kmax )
integer :: binidx ( 0 : kmax )
integer :: status
logical :: calc_phymix_fwf
! EOP
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! BOC
@ -361,6 +401,15 @@
# endif
call tic ( TIM_SALT )
if ( allocated ( phymix_S_riv_s ) ) then
do n = 1 , nriverl
rivmix ( : , n ) = phymix_S ( ri ( n ) , rj ( n ) , : )
end do
end if
calc_phymix_fwf = ( associated ( phymix_S ) . or . allocated ( phymix_S_fwf ) )
if ( calc_phymix_fwf ) sss_old = S ( : , : , kmax )
do j = jmin - HALO , jmax + HALO
do i = imin - HALO , imax + HALO
if ( az ( i , j ) . eq . 1 ) then
@ -369,28 +418,202 @@
end do
end do
if ( calc_phymix_fwf ) then
! Note ( KK ) : division by hn because we already refer to new volume
! where diffusion will be added
where ( az . eq . 1 ) wrk2d = sss_old * S ( : , : , kmax ) * fwf_int / dt
if ( associated ( phymix_S ) ) then
where ( az . eq . 1 ) phymix_S ( : , : , kmax ) = phymix_S ( : , : , kmax ) + wrk2d
end if
if ( allocated ( phymix_S_fwf ) ) then
where ( az . eq . 1 ) phymix_S_fwf = wrk2d
end if
end if
if ( allocated ( fwf_s ) . or . allocated ( fwfS2_s ) ) then
call tic ( TIM_TEF )
if ( allocated ( fwf_s ) ) fwf_s = _ ZERO_
if ( allocated ( fwfS2_s ) ) fwfS2_s = _ ZERO_
do j = jmin , jmax
do i = imin , imax
if ( az ( i , j ) . eq . 1 ) then
l = binidx_kij ( kmax , i , j ) ! binning based on old salinity
if ( l . ge . 0 ) then
if ( allocated ( fwf_s ) ) fwf_s ( i , j , l ) = fwf_int ( i , j )
if ( allocated ( fwfS2_s ) ) fwfS2_s ( i , j , l ) = fwf_int ( i , j ) * S ( i , j , kmax ) * S ( i , j , kmax )
end if
end if
end do
end do
call toc ( TIM_TEF )
end if
call do_advection_3d ( dt , S , uu , vv , ww , hun , hvn , ho , hn , &
salt_adv_split , salt_adv_hor , salt_adv_ver , _ ZERO_ , H_TAG , &
ffluxu = Sfluxu , ffluxv = Sfluxv , ffluxw = Sfluxw , nvd = nummix_S )
if ( allocated ( Sfluxu2 ) ) Sfluxu2 = ffluxu2_3d
if ( allocated ( Sfluxv2 ) ) Sfluxv2 = ffluxv2_3d
! KK - TODO : do we have to treat salt fluxes from different
! directional - split steps separately ?
if ( allocated ( hu_s ) . or . allocated ( uu_s ) . or . allocated ( Sfluxu_s ) . or . allocated ( S2fluxu_s ) ) then
call tic ( TIM_TEF )
if ( allocated ( hu_s ) ) hu_s = _ ZERO_
if ( allocated ( uu_s ) ) uu_s = _ ZERO_
if ( allocated ( Sfluxu_s ) ) Sfluxu_s = _ ZERO_
if ( allocated ( S2fluxu_s ) ) S2fluxu_s = _ ZERO_
do j = jmin , jmax
do i = imin , imax
if ( au ( i , j ) . eq . 1 . or . au ( i , j ) . eq . 2 ) then
where ( uu ( i , j , 1 : kmax ) . ne . _ ZERO_ )
S_k ( 1 : kmax ) = Sfluxu ( i , j , 1 : kmax ) / uu ( i , j , 1 : kmax ) / DYU
else where
S_k ( 1 : kmax ) = saltbinmax + 999 9.0_rk
end where
call rebin ( kmax , S_k , saltbins , Si_s , binidx )
do k = 1 , kmax
l = binidx ( k )
if ( l . lt . 0 ) cycle
if ( allocated ( hu_s ) ) hu_s ( i , j , l ) = hu_s ( i , j , l ) + hun ( i , j , k )
if ( allocated ( uu_s ) ) uu_s ( i , j , l ) = uu_s ( i , j , l ) + uu ( i , j , k )
if ( allocated ( Sfluxu_s ) ) Sfluxu_s ( i , j , l ) = Sfluxu_s ( i , j , l ) + Sfluxu ( i , j , k )
if ( allocated ( S2fluxu_s ) ) S2fluxu_s ( i , j , l ) = S2fluxu_s ( i , j , l ) + Sfluxu2 ( i , j , k )
end do
end if
end do
end do
call toc ( TIM_TEF )
end if
if ( allocated ( hv_s ) . or . allocated ( vv_s ) . or . allocated ( Sfluxv_s ) . or . allocated ( S2fluxv_s ) ) then
call tic ( TIM_TEF )
if ( allocated ( hv_s ) ) hv_s = _ ZERO_
if ( allocated ( vv_s ) ) vv_s = _ ZERO_
if ( allocated ( Sfluxv_s ) ) Sfluxv_s = _ ZERO_
if ( allocated ( S2fluxv_s ) ) S2fluxv_s = _ ZERO_
do j = jmin , jmax
do i = imin , imax
if ( av ( i , j ) . eq . 1 . or . av ( i , j ) . eq . 2 ) then
where ( vv ( i , j , 1 : kmax ) . ne . _ ZERO_ )
S_k ( 1 : kmax ) = Sfluxv ( i , j , 1 : kmax ) / vv ( i , j , 1 : kmax ) / DXV
else where
S_k ( 1 : kmax ) = saltbinmax + 999 9.0_rk
end where
call rebin ( kmax , S_k , saltbins , Si_s , binidx )
do k = 1 , kmax
l = binidx ( k )
if ( l . lt . 0 ) cycle
if ( allocated ( hv_s ) ) hv_s ( i , j , l ) = hv_s ( i , j , l ) + hvn ( i , j , k )
if ( allocated ( vv_s ) ) vv_s ( i , j , l ) = vv_s ( i , j , l ) + vv ( i , j , k )
if ( allocated ( Sfluxv_s ) ) Sfluxv_s ( i , j , l ) = Sfluxv_s ( i , j , l ) + Sfluxv ( i , j , k )
if ( allocated ( S2fluxv_s ) ) S2fluxv_s ( i , j , l ) = S2fluxv_s ( i , j , l ) + Sfluxv2 ( i , j , k )
end do
end if
end do
end do
call toc ( TIM_TEF )
end if
if ( salt_AH_method . gt . 0 ) then
! S is not halo updated after advection
call tic ( TIM_SALTH )
call update_3d_halo ( S , S , az , imin , jmin , imax , jmax , kmax , D_TAG )
call wait_halo ( D_TAG )
call toc ( TIM_SALTH )
if ( allocated ( hu_s ) . or . allocated ( uu_s ) . or . allocated ( Sfluxu_s ) . or . allocated ( S2fluxu_s ) &
. or . allocated ( hv_s ) . or . allocated ( vv_s ) . or . allocated ( Sfluxv_s ) . or . allocated ( S2fluxv_s ) ) then
call tic ( TIM_TEF )
Sold = S
do j = jmin , jmax + 1
do i = imin , imax + 1
if ( az ( i , j ) . ne . 0 ) call rebin ( kmax , S ( i , j , : ) , saltbins , Si_s , binidx_kij ( : , i , j ) )
end do
end do
call toc ( TIM_TEF )
end if
if ( allocated ( hu_s ) . or . allocated ( uu_s ) . or . allocated ( Sfluxu_s ) ) Sfluxu_tmp = Sfluxu
if ( allocated ( hv_s ) . or . allocated ( vv_s ) . or . allocated ( Sfluxv_s ) ) Sfluxv_tmp = Sfluxv
call tracer_diffusion ( S , hn , salt_AH_method , salt_AH_const , salt_AH_Prt , salt_AH_stirr_const , &
ffluxu = Sfluxu , ffluxv = Sfluxv , &
phymix = phymix_S )
end if
if ( allocated ( Sfluxu2 ) ) Sfluxu2 = Sfluxu2 + ffluxu2_3d
if ( allocated ( Sfluxv2 ) ) Sfluxv2 = Sfluxv2 + ffluxv2_3d
if ( allocated ( hu_s ) . or . allocated ( uu_s ) . or . allocated ( Sfluxu_s ) . or . allocated ( S2fluxu_s ) ) then
call tic ( TIM_TEF )
Sfluxu_tmp = Sfluxu - Sfluxu_tmp
do j = jmin , jmax
do i = imin , imax
if ( au ( i , j ) . eq . 1 . or . au ( i , j ) . eq . 2 ) then
do k = 1 , kmax
if ( Sfluxu_tmp ( i , j , k ) . ne . _ ZERO_ ) then
deltaS = Sold ( i + 1 , j , k ) - Sold ( i , j , k )
! if ( deltaS . ne . _ ZERO_ ) then
hu_k = _ HALF_ * ( hn ( i , j , k ) + hn ( i + 1 , j , k ) )
uu_k = Sfluxu_tmp ( i , j , k ) / deltaS / DYU
! fluxes with S ( i )
l = binidx_kij ( k , i , j )
if ( l . ge . 0 ) then
if ( allocated ( hu_s ) ) hu_s ( i , j , l ) = hu_s ( i , j , l ) + hu_k
if ( allocated ( uu_s ) ) uu_s ( i , j , l ) = uu_s ( i , j , l ) - uu_k
if ( allocated ( Sfluxu_s ) ) Sfluxu_s ( i , j , l ) = Sfluxu_s ( i , j , l ) - uu_k * DYU * Sold ( i , j , k )
if ( allocated ( S2fluxu_s ) ) S2fluxu_s ( i , j , l ) = S2fluxu_s ( i , j , l ) - uu_k * DYU * Sold ( i , j , k ) * Sold ( i , j , k )
end if
! fluxes with S ( i + 1 )
l = binidx_kij ( k , i + 1 , j )
if ( l . ge . 0 ) then
if ( allocated ( hu_s ) ) hu_s ( i , j , l ) = hu_s ( i , j , l ) + hu_k
if ( allocated ( uu_s ) ) uu_s ( i , j , l ) = uu_s ( i , j , l ) + uu_k
if ( allocated ( Sfluxu_s ) ) Sfluxu_s ( i , j , l ) = Sfluxu_s ( i , j , l ) + uu_k * DYU * Sold ( i + 1 , j , k )
if ( allocated ( S2fluxu_s ) ) S2fluxu_s ( i , j , l ) = S2fluxu_s ( i , j , l ) + uu_k * DYU * Sold ( i + 1 , j , k ) * Sold ( i + 1 , j , k )
end if
end if
end do
end if
end do
end do
call toc ( TIM_TEF )
end if
if ( allocated ( hv_s ) . or . allocated ( vv_s ) . or . allocated ( Sfluxv_s ) . or . allocated ( S2fluxv_s ) ) then
call tic ( TIM_TEF )
Sfluxv_tmp = Sfluxv - Sfluxv_tmp
do j = jmin , jmax
do i = imin , imax
if ( av ( i , j ) . eq . 1 . or . av ( i , j ) . eq . 2 ) then
do k = 1 , kmax
if ( Sfluxv_tmp ( i , j , k ) . ne . _ ZERO_ ) then
deltaS = Sold ( i , j + 1 , k ) - Sold ( i , j , k )
! if ( deltaS . ne . _ ZERO_ ) then
hu_k = _ HALF_ * ( hn ( i , j , k ) + hn ( i , j + 1 , k ) )
uu_k = Sfluxv_tmp ( i , j , k ) / deltaS / DXV
! fluxes with S ( j )
l = binidx_kij ( k , i , j )
if ( l . ge . 0 ) then
if ( allocated ( hv_s ) ) hv_s ( i , j , l ) = hv_s ( i , j , l ) + hu_k
if ( allocated ( vv_s ) ) vv_s ( i , j , l ) = vv_s ( i , j , l ) - uu_k
if ( allocated ( Sfluxv_s ) ) Sfluxv_s ( i , j , l ) = Sfluxv_s ( i , j , l ) - uu_k * DXV * Sold ( i , j , k )
if ( allocated ( S2fluxv_s ) ) S2fluxv_s ( i , j , l ) = S2fluxv_s ( i , j , l ) - uu_k * DXV * Sold ( i , j , k ) * Sold ( i , j , k )
end if
! fluxes with S ( j + 1 )
l = binidx_kij ( k , i , j + 1 )
if ( l . ge . 0 ) then
if ( allocated ( hv_s ) ) hv_s ( i , j , l ) = hv_s ( i , j , l ) + hu_k
if ( allocated ( vv_s ) ) vv_s ( i , j , l ) = vv_s ( i , j , l ) + uu_k
if ( allocated ( Sfluxv_s ) ) Sfluxv_s ( i , j , l ) = Sfluxv_s ( i , j , l ) + uu_k * DXV * Sold ( i , j + 1 , k )
if ( allocated ( S2fluxv_s ) ) S2fluxv_s ( i , j , l ) = S2fluxv_s ( i , j , l ) + uu_k * DXV * Sold ( i , j + 1 , k ) * Sold ( i , j + 1 , k )
end if
end if
end do
end if
end do
end do
call toc ( TIM_TEF )
end if
if ( associated ( phymix_S ) ) then
call toc ( TIM_SALT )
call tic ( TIM_MIXANALYSIS )
call physical_mixing ( S , avmols , phymix_S , salt_AH_method )
call toc ( TIM_MIXANALYSIS )
call tic ( TIM_SALT )
end if
@ -417,7 +640,7 @@
if ( rc / = 0 ) stop 'do_salinity: Error allocating memory (auxo)'
fluxw = > null ( )
if ( associated ( Sfluxw ) ) then
if ( associated ( Sfluxw ) . or . associated ( phymix_S ) ) then
allocate ( fluxw ( 0 : kmax ) , stat = rc ) ! work array
if ( rc / = 0 ) stop 'do_salinity: Error allocating memory (fluxw)'
end if
@ -486,6 +709,28 @@
if ( associated ( Sfluxw ) ) then
Sfluxw ( i , j , 1 : kmax - 1 ) = Sfluxw ( i , j , 1 : kmax - 1 ) + fluxw ( 1 : kmax - 1 )
end if
if ( associated ( phymix_S ) ) then
! KK - TODO : only consider changes of S2 due to diffusion !
! ( no sss_nudging or radiation for temperature )
do k = 1 , kmax - 1
# if 0
fluxw ( k ) = fluxw ( k ) * &
( cnpar * ( Res ( k ) + Res ( k + 1 ) ) &
+ ( _ ONE_ - cnpar ) * ( S ( i , j , k ) + S ( i , j , k + 1 ) ) )
# else
fluxw ( k ) = - ( &
auxo ( k ) * ( S ( i , j , k + 1 ) - S ( i , j , k ) ) * ( S ( i , j , k ) + S ( i , j , k + 1 ) ) &
+ auxn ( k ) * ( Res ( k + 1 ) - Res ( k ) ) * ( Res ( k ) + Res ( k + 1 ) ) &
) / dt
# endif
end do
do k = 1 , kmax
phymix_S ( i , j , k ) = phymix_S ( i , j , k ) &
- ( hn ( i , j , k ) * ( Res ( k ) * Res ( k ) - S ( i , j , k ) * S ( i , j , k ) ) / dt &
+ ( fluxw ( k ) - fluxw ( k - 1 ) ) )
end do
end if
do k = 1 , kmax
S ( i , j , k ) = Res ( k )
@ -520,8 +765,55 @@
! $ OMP END PARALLEL
if ( allocated ( counts_s ) . or . allocated ( flags_s ) . or . allocated ( h_s ) . or . allocated ( hS_s ) . or . allocated ( hS2_s ) . or . allocated ( hpmS_s ) . or . allocated ( hnmS_s ) . or . allocated ( phymix_S_fwf_s ) . or . allocated ( phymix_S_riv_s ) ) then
call tic ( TIM_TEF )
if ( allocated ( counts_s ) ) counts_s = 0
if ( allocated ( flags_s ) ) flags_s = 0
if ( allocated ( h_s ) ) h_s = _ ZERO_
if ( allocated ( hS_s ) ) hS_s = _ ZERO_
if ( allocated ( hS2_s ) ) hS2_s = _ ZERO_
if ( allocated ( hpmS_s ) ) hpmS_s = _ ZERO_
if ( allocated ( hnmS_s ) ) hnmS_s = _ ZERO_
if ( allocated ( phymix_S_fwf_s ) ) phymix_S_fwf_s = _ ZERO_
if ( allocated ( phymix_S_riv_s ) ) phymix_S_riv_s = _ ZERO_
do j = jmin , jmax
do i = imin , imax
if ( az ( i , j ) . eq . 1 ) then
call rebin ( kmax , S ( i , j , : ) , saltbins , Si_s , binidx )
binidx_kij ( : , i , j ) = binidx ( : )
do k = 1 , kmax
l = binidx ( k )
if ( l . lt . 0 ) cycle
if ( allocated ( counts_s ) ) counts_s ( i , j , l ) = counts_s ( i , j , l ) + 1
if ( allocated ( flags_s ) ) flags_s ( i , j , l ) = 1
if ( allocated ( h_s ) ) h_s ( i , j , l ) = h_s ( i , j , l ) + hn ( i , j , k )
if ( allocated ( hS_s ) ) hS_s ( i , j , l ) = hS_s ( i , j , l ) + hn ( i , j , k ) * S ( i , j , k )
if ( allocated ( hS2_s ) ) hS2_s ( i , j , l ) = hS2_s ( i , j , l ) + hn ( i , j , k ) * S ( i , j , k ) * S ( i , j , k )
if ( allocated ( hpmS_s ) ) hpmS_s ( i , j , l ) = hpmS_s ( i , j , l ) + phymix_S ( i , j , k )
if ( allocated ( hnmS_s ) ) hnmS_s ( i , j , l ) = hnmS_s ( i , j , l ) + nummix_S ( i , j , k )
end do
l = binidx ( kmax )
if ( l . ge . 0 ) then
if ( allocated ( phymix_S_fwf_s ) ) phymix_S_fwf_s ( i , j , l ) = phymix_S_fwf ( i , j )
end if
end if
end do
end do
if ( allocated ( phymix_S_riv_s ) ) then
do n = 1 , nriverl
j = rj ( n )
i = ri ( n )
do k = 1 , kmax
l = binidx_kij ( k , i , j )
if ( l . lt . 0 ) cycle
phymix_S_riv_s ( i , j , l ) = phymix_S_riv_s ( i , j , l ) + rivmix ( k , n )
end do
end do
end if
call toc ( TIM_TEF )
end if
if ( salt_check . ne . 0 . and . mod ( n , abs ( salt_check ) ) . eq . 0 ) then
if ( salt_check . ne . 0 . and . mod ( loop , abs ( salt_check ) ) . eq . 0 ) then
call check_3d_fields ( imin , jmin , imax , jmax , kmin , kmax , az , &
S , min_salt , max_salt , status )
if ( status . gt . 0 ) then