vdif_q.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	vdif_q

Function/Subroutine Documentation

vdif_q()

subroutine vdif_q

(

)

Definition at line 46 of file vdif_q.f90.

!------------------------------------------------------------------------------|
   USE all_vars
   USE mod_utils
   USE mod_wd
   USE mod_par



   IMPLICIT NONE
   INTEGER :: I,J,K,KI
   REAL(SP)  :: CONST1,COEF1,COEF2,COEF3,COEF4,COEF5,LMAX
   REAL(SP), DIMENSION(0:MT,KB) :: GH,SM,SH,PROD,DTEF,KN,BOYGR,A,C,VHP,VH,STF
   REAL(SP), DIMENSION(0:MT) :: L0
   REAL(SP)  :: UTAU2
   REAL(SP)  :: WUSURF_NODE,WVSURF_NODE,WUBOT_NODE,WVBOT_NODE
   REAL(SP)  :: UU_NODE_K,VV_NODE_K,UU_NODE_KM1,VV_NODE_KM1

!--Stability Function Coefficients-------------
   REAL(SP), PARAMETER :: A1    = 0.92_sp
   REAL(SP), PARAMETER :: B1    = 16.6_sp
   REAL(SP), PARAMETER :: A2    = 0.74_sp
   REAL(SP), PARAMETER :: B2    = 10.1_sp
   REAL(SP), PARAMETER :: C1    = 0.08_sp

!--Source/Sink Term Coefficients---------------
   REAL(SP), PARAMETER :: E1    = 1.80_sp
   REAL(SP), PARAMETER :: E2    = 1.33_sp
   REAL(SP), PARAMETER :: SEF   = 1.00_sp

!--Von Karmans Constant------------------------
   REAL(SP), PARAMETER :: KAPPA = 0.40_sp

!--Gravitational Constant----------------------
   REAL(SP), PARAMETER :: GEE   = 9.806_sp

!--Limiting Values of Q2/Q2L-------------------
   REAL(SP), PARAMETER :: SMALL = 1.e-8_sp

!--Coefficient of buoyancy length scale--------
   REAL(SP), PARAMETER :: CB    = 1.e-8_sp

!--Denominator small number threshold----------
   REAL(SP), PARAMETER :: DEN_S = 1.e-10_sp

!--Upper bound of GH for unstable flow---------
   REAL(SP), PARAMETER :: GH_MAX = .0233_sp

!--Lower bound for GH for stable flow----------
   REAL(SP), PARAMETER :: GH_MIN = -.281_sp



!----------------------------------------------
   REAL(SP), PARAMETER :: CBCNST = 100.0_sp
   REAL(SP), PARAMETER :: SURFL = 2.e+5_sp
   REAL(SP), PARAMETER :: SHIW = 0.0_sp
   REAL(SP), PARAMETER :: GHC = -6.0_sp

   REAL(SP)  :: KQ_TMP,KM_TMP,KH_TMP

   IF(dbg_set(dbg_sbr)) write(ipt,*) "Start: vdif_q"
!------------------------------------------------------------------------------|

   ! INITIALIZE MEMORY
   boygr = 0.0_sp
   gh    = 0.0_sp
   sm    = 0.0_sp
   sh    = 0.0_sp
   prod  = 0.0_sp
   dtef  = 0.0_sp
   kn    = 0.0_sp
   a     = 0.0_sp
   c     = 0.0_sp
   vhp   = 0.0_sp
   vh    = 0.0_sp
   stf   = 0.0_sp
   l0    = 0.0_sp

!------------------------------------------------------------------------------|
!  place a threshold on the turbulence quantities following advection          |
!------------------------------------------------------------------------------|

   DO k = 2, kbm1
     DO i = 1, m
       IF (q2f(i,k) <= small .OR. q2lf(i,k) <= small) THEN
!         Q2F(I,K)  = 0.0_SP
!         Q2LF(I,K) = 0.0_SP
         q2f(i,k)  = small
         q2lf(i,k) = small
       END IF
     END DO
   END DO

!------------------------------------------------------------------------------!
!  set up coefficients for implicit calculation of vertical diffusion          !
!------------------------------------------------------------------------------!

   DO i=1,m
       IF(iswetn(i) == 1)THEN
     DO k=2,kbm1
       a(i,k) = -0.5_sp* dti * (kq(i,k+1)+kq(i,k)+2.0_sp*umol)/ &
             (dzz(i,k-1)*dz(i,k)*d(i)*d(i))
       c(i,k) = -0.5_sp* dti*(kq(i,k-1)+kq(i,k)+2.0_sp*umol)/ &
             (dzz(i,k-1)*dz(i,k-1)*d(i)*d(i))
     END DO
     END IF
   END DO

!------------------------------------------------------------------------------!
!    the following section solves the equation                                 !
!    dti*(kq*q2')' - q2*(2.*dti*dtef+1.) = -q2b                                !
!------------------------------------------------------------------------------!

   const1 = 16.6_sp ** .6666667_sp * sef
   DO i = 1, m
!       VHP(I,1) = SQRT(WUSURF(I)**2+WVSURF(I)**2) * CONST1
!      VH(I,1)  = 0.0_SP

     wusurf_node = 0.0_sp
     wvsurf_node = 0.0_sp
     DO j=1,ntve(i)
       wusurf_node = wusurf_node + wusurf(nbve(i,j))
       wvsurf_node = wvsurf_node + wvsurf(nbve(i,j))
     END DO
     wusurf_node = wusurf_node/float(ntve(i))
     wvsurf_node = wvsurf_node/float(ntve(i))

     utau2 = sqrt(wusurf_node**2+wvsurf_node**2)
     vhp(i,1) = (15.8_sp*cbcnst)**0.6666667_sp*utau2
     vh(i,1)  = 0.0_sp
!     L0(I)    = SURFL*UTAU2/GEE
     l0(i)    = surfl*utau2/grav_n(i)
!     HLC   = 2.4_SP          
!      L0(I) = HLC*HSC1(I)            


     wubot_node = 0.0_sp
     wvbot_node = 0.0_sp
     DO j=1,ntve(i)
       wubot_node = wubot_node + wubot(nbve(i,j))
       wvbot_node = wvbot_node + wvbot(nbve(i,j))
     END DO
     wubot_node = wubot_node/float(ntve(i))
     wvbot_node = wvbot_node/float(ntve(i))

     q2f(i,kb) = sqrt(wubot_node**2+wvbot_node**2) * const1
   END DO

   q2  = abs(q2+small)
   q2l = abs(q2l)
 
!------------------------------------------------------------------------------!
!  calculate boygr = -Brunt Vaisala Frequency^2                                !
!  calculate internal wave shear energy contribution: prod = 200*(b.v.freq)**3 !
!------------------------------------------------------------------------------!

   DO k=2,kbm1
     DO i=1,m
      IF(iswetn(i) == 1)THEN
!       BOYGR(I,K) = GEE * (RHO1(I,K-1)-RHO1(I,K))/(DZZ(I,K-1)*D(I))
       boygr(i,k) = grav_n(i) * (rho1(i,k-1)-rho1(i,k))/(dzz(i,k-1)*d(i))
       prod(i,k)  = km(i,k) * 0.0_sp * (.5_sp*(-boygr(i,k)+abs(boygr(i,k)))) ** 1.5_sp
      END IF
     END DO
   END DO

!------------------------------------------------------------------------------!
!  calculate shear production source term = prod                               !
!  calculate buoyancy production source term = kh*boygr                        !
!------------------------------------------------------------------------------!
   DO  k = 2, kbm1
     DO  i = 1, m
       IF(iswetn(i) == 1)THEN

         uu_node_k   = 0.0_sp
         vv_node_k   = 0.0_sp
         uu_node_km1 = 0.0_sp
         vv_node_km1 = 0.0_sp
         DO j=1,ntve(i)
           uu_node_k   = uu_node_k + u(nbve(i,j),k)
           vv_node_k   = vv_node_k + v(nbve(i,j),k)
           uu_node_km1 = uu_node_km1 + u(nbve(i,j),k-1)
           vv_node_km1 = vv_node_km1 + v(nbve(i,j),k-1)
         END DO
         uu_node_k   = uu_node_k/float(ntve(i))
         vv_node_k   = vv_node_k/float(ntve(i))
         uu_node_km1 = uu_node_km1/float(ntve(i))
         vv_node_km1 = vv_node_km1/float(ntve(i))

         prod(i,k) = prod(i,k) + km(i,k) * sef * &
                     ((uu_node_k-uu_node_km1)**2+(vv_node_k-vv_node_km1)**2)/ &
                     (dzz(i,k-1)*d(i))**2
         prod(i,k) = prod(i,k) + kh(i,k) * boygr(i,k)
       END IF
     END DO
   END DO
!------------------------------------------------------------------------------!
!  solve for turbulent length scale l = q2l/q2                                 !
!------------------------------------------------------------------------------!

!   L = Q2L/Q2
    DO k=2,kbm1
      DO i=1,m
        l(i,k) = q2l(i,k)/q2(i,k)
      END DO
    END DO

!------------------------------------------------------------------------------!
!  in stably stratified regions, length scale is limited by buoyancy length    !
!  scale, l_b = c_b*(k^.5/N).  length scale limitation also puts an effective  !
!  lower bound of -.281 on gh in stably stratified regions                     !
!  length scale near surface is modified with wind mixed length scale          !
!------------------------------------------------------------------------------!

   DO k=2,kbm1
     DO i=1,m
       IF(z(i,k) > -0.5_sp)l(i,k) = max(l(i,k),kappa*l0(i)) 
       IF(boygr(i,k) < 0.0_sp)THEN
         lmax      = sqrt(abs(gh_min) * q2(i,k) / (-boygr(i,k) + small) )
         l(i,k)    = min(l(i,k),lmax)
         q2l(i,k) = q2(i,k)*l(i,k)
       END IF
     END DO
   END DO

!------------------------------------------------------------------------------!
!  solve for gh = (-l^2*N^2/q^2)                                               !
!------------------------------------------------------------------------------!

!!   GH = (L**2/Q2)*BOYGR
!   GH(1:MT,1:KBM1) = (L(1:MT,1:KBM1)**2/Q2(1:MT,1:KBM1))*BOYGR(1:MT,1:KBM1)
   gh(1:m,1:kbm1) = (l(1:m,1:kbm1)**2/q2(1:m,1:kbm1))*boygr(1:m,1:kbm1)

!------------------------------------------------------------------------------!
!  limit maximum of GH in unstable regions to reasonable physical value        !
!  limiting GH also constrains stability functions SH/SM to finite values      !
!------------------------------------------------------------------------------!

   gh = min(gh,gh_max)

   l(:,1)   = kappa*l0(:)         !0.0_SP
   l(:,kb)  = 0.0_sp
   gh(:,1)  = 0.0_sp
   gh(:,kb) = 0.0_sp


!------------------------------------------------------------------------------!
!  calculate eddy kinetic energy dissipation rate: dtef                        !
!------------------------------------------------------------------------------!

   stf = 1.0_sp

   IF(.not.surface_wave_mixing)THEN
   DO k = 1,kb
     DO i = 1,m
       IF(gh(i,k) < 0.0_sp) stf(i,k)=1.0-0.9*(gh(i,k)/ghc)**1.5
       IF(gh(i,k) < ghc) stf(i,k) = 0.1
     END DO
   END DO
   END IF

!   DTEF = Q2*SQRT(Q2)/(B1*Q2L + small)*STF
   dtef = sqrt(q2)/(b1*l + small)*stf

   DO k = 2, kbm1
     DO i = 1, m
       vhp(i,k) = 1. / (a(i,k)+c(i,k)*(1.-vh(i,k-1))-(2.*dti*dtef(i,k)+1.))
       vh(i,k) = a(i,k) * vhp(i,k)
       vhp(i,k) = (-2.*dti*prod(i,k)+c(i,k)*vhp(i,k-1)-q2f(i,k))*vhp(i,k)
     END DO
   END DO

   DO k=1,kbm1
     ki = kb-k
     DO i = 1, m
     IF(iswetn(i) == 1)THEN
       q2f(i,ki) = vh(i,ki) * q2f(i,ki+1) + vhp(i,ki)
     END IF
     END DO
   END DO

!------------------------------------------------------------------------------!
!      the following section solves the equation                               !
!      dti*(kq*q2l')' - q2l*(dti*dtef+1.) = -q2lb                              !
!------------------------------------------------------------------------------!

   vh(:,1)  = 0.0_sp
   vhp(:,1) = 0.0_sp

   DO  k = 2, kbm1
     DO  i = 1, m
       IF(iswetn(i) == 1)THEN
         dtef(i,k) = dtef(i,k) * (1.+e2*((1./abs(z(i,k)-z(i,1))+1./ &
                     abs(z(i,k)-z(i,kb)))*l(i,k)/(d(i)*kappa))**2)
         vhp(i,k)  = 1. / (a(i,k)+c(i,k)*(1.-vh(i,k-1))- &
                     (dti*dtef(i,k)+1.))
         vh(i,k)   = a(i,k) * vhp(i,k)
         vhp(i,k)  = (dti*(-prod(i,k)*l(i,k)*e1)+c(i,k)* &
                     vhp(i,k-1)-q2lf(i,k)) * vhp(i,k)
       END IF
     END DO
   END DO


   DO k=1,kbm1
     ki = kb - k
     DO i = 1, m
     IF(iswetn(i) == 1)THEN
       q2lf(i,ki) = vh(i,ki) * q2lf(i,ki+1) + vhp(i,ki)
     END IF
     END DO
   END DO

!------------------------------------------------------------------------------|
!  place a threshold on the turbulence quantities following vertical diffusion |
!------------------------------------------------------------------------------|


   DO k = 2, kbm1
     DO i = 1, m
       IF (q2f(i,k) <= small .OR. q2lf(i,k) <= small) THEN
         q2f(i,k)  = small
         q2lf(i,k) = small
       END IF
     END DO
   END DO


!===========THRESHOLD ON LENGTH SCALE ALT STYLE (Active in Barotropic Case)====
!   DO K=2,KBM1
!     DO I=1,N
!       LMAX      = SQRT(ABS(GH_MIN) * Q2F(I,K) / MAX(0.0_SP,-BOYGR(I,K)+SMALL))
!       L(I,K)    = MIN(L(I,K),LMAX)
!       Q2LF(I,K) = Q2F(I,K)*L(I,K)
!     END DO
!   END DO
!------------------------------------------------------------------------------!

!   L(:,1)   = 0.0_SP
!   L(:,KB)  = 0.0_SP
!   GH(:,1)  = 0.0_SP
!   GH(:,KB) = 0.0_SP


!------------------------------------------------------------------------------!
!  calculate stability functions sh + sm using Galperin 1988 formulation       !
!------------------------------------------------------------------------------!
   coef4 = 18.0_sp * a1 * a1 + 9.0_sp * a1 * a2
   coef5 = 9.0_sp * a1 * a2

   DO k = 1,kb
     DO i = 1,m
       coef1 = a2 * (1.0_sp-6.0_sp*a1/b1*stf(i,k))
       coef2 = 3.0_sp * a2 * b2/stf(i,k) + 18.0_sp * a1 * a2
       coef3 = a1 * (1.0_sp-3.0_sp*c1-6.0_sp*a1/b1*stf(i,k))

       sh(i,k) = coef1/(1.0_sp-coef2*gh(i,k))
       sm(i,k) = (coef3+sh(i,k)*coef4*gh(i,k))/(1.-coef5*gh(i,k))
     END DO
   END DO

!------------------------------------------------------------------------------!
!  calculate turbulent eddy viscosities/diffusivities kq,km,kh                 !
!------------------------------------------------------------------------------!

   kn = l*sqrt(q2)
   kq = 0.5_sp*(kn*kappa*sm+kq)
   km = 0.5_sp*(kn*sm+km)
   kh = 0.5_sp*(kn*sh*vprnu+kh)

    DO i=1,m
     IF(h(i) <= static_ssh_adj)THEN
      kq_tmp = 0.0_sp
      km_tmp = 0.0_sp
      kh_tmp = 0.0_sp
      DO k=2,kbm1
        kq_tmp = kq_tmp + kq(i,k)
        km_tmp = km_tmp + km(i,k)
        kh_tmp = kh_tmp + kh(i,k)
      END DO
        kq_tmp = kq_tmp/(kbm1-1)
        km_tmp = km_tmp/(kbm1-1)
        kh_tmp = kh_tmp/(kbm1-1)
      DO k=2,kbm1
        kq(i,k) = kq_tmp
        km(i,k) = km_tmp
        kh(i,k) = kh_tmp
      END DO
     END IF
    END DO




   IF(dbg_set(dbg_sbr)) write(ipt,*) "End: vdif_q"

Here is the call graph for this function:

Here is the caller graph for this function: