source: branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2872

MODULE dynzdf_imp
   !!======================================================================
   !!                    ***  MODULE  dynzdf_imp  ***
   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
   !!======================================================================
   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dyn_zdf_imp  : update the momentum trend with the vertical diffusion using a implicit time-stepping
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE sbc_oce         ! surface boundary condition: ocean
   USE zdf_oce         ! ocean vertical physics
   USE phycst          ! physical constants
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! MPP library
   USE zdfbfr          ! bottom friction setup

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dyn_zdf_imp   ! called by step.F90

   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id$
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE dyn_zdf_imp( kt, p2dt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_zdf_imp  ***
      !!                   
      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
      !!      and the surface forcing, and add it to the general trend of 
      !!      the momentum equations.
      !!
      !! ** Method  :   The vertical momentum mixing trend is given by :
      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
      !!      backward time stepping
      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
      !!      Add this trend to the general trend ua :
      !!         ua = ua + dz( avmu dz(u) )
      !!
      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
      !!---------------------------------------------------------------------
      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
      USE oce     , ONLY:  zwd  => ta       , zws   => sa   ! (ta,sa) used as 3D workspace
      USE wrk_nemo, ONLY:   zwi => wrk_3d_3                 ! 3D workspace
      !!
      INTEGER , INTENT(in) ::   kt    ! ocean time-step index
      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
      !!
      INTEGER  ::   ji, jj, jk     ! dummy loop indices
      INTEGER  ::   ikbum1, ikbvm1 ! local variable
      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs   ! local scalars

      !! * Local variables for implicit bottom friction.    H. Liu
      REAL(wp) ::   zbfru, zbfrv 
      REAL(wp) ::   bfr_imp = 1._wp
      !!----------------------------------------------------------------------
      !!----------------------------------------------------------------------

      IF( wrk_in_use(3, 3) ) THEN
         CALL ctl_stop('dyn_zdf_imp: requested workspace array unavailable')   ;   RETURN
      END IF

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
         IF(.NOT.ln_bfrimp) bfr_imp = 0._wp
      ENDIF

      ! 0. Local constant initialization
      ! --------------------------------
      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep

      ! 1. Vertical diffusion on u

      ! Vertical diffusion on u&v
      ! ---------------------------
      ! Matrix and second member construction
      !! bottom boundary condition: both zwi and zws must be masked as avmu can take
      !! non zero value at the ocean bottom depending on the bottom friction
      !! used but the bottom velocities have already been updated with the bottom
      !! friction velocity in dyn_bfr using values from the previous timestep. There
      !! is no need to include these in the implicit calculation.

      ! The code has been modified here to implicitly implement bottom
      ! friction: u(v)mask is not necessary here anymore. 
      ! H. Liu, April 2010.

      ! 1. Vertical diffusion on u
      DO jj = 2, jpjm1 
         DO ji = fs_2, fs_jpim1   ! vector opt.
            ikbum1 = mbku(ji,jj)
            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
!               zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbum1)*zinv )
               zbfru = bfrua(ji,jj)

            DO jk = 1, ikbum1
               zcoef = - p2dt / fse3u(ji,jj,jk)
               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
            END DO

      ! Surface boudary conditions
            zwi(ji,jj,1) = 0._wp
            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)

      ! Bottom boudary conditions  ! H. Liu, May, 2010
!            !commented out to be consisent with v3.2, h.liu
!            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0 * bfr_imp
            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * bfr_imp
            zws(ji,jj,ikbum1) = 0._wp
            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 

      ! Matrix inversion starting from the first level
      !-----------------------------------------------------------------------
      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
      !
      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
      !        (        ...               )( ...  ) ( ...  )
      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
      !
      !   m is decomposed in the product of an upper and a lower triangular matrix
      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
      !   The solution (the after velocity) is in ua
      !-----------------------------------------------------------------------

      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
            DO jk = 2, ikbum1
               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
            END DO

      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 )
            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj))
           DO jk = 2, ikbum1   
               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
            END DO


      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1)
            DO jk = ikbum1-1, 1, -1
               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
            END DO
         END DO
      END DO

      ! Normalization to obtain the general momentum trend ua
      DO jk = 1, jpkm1
         DO jj = 2, jpjm1   
            DO ji = fs_2, fs_jpim1   ! vector opt.
               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
            END DO
         END DO
      END DO


      ! 2. Vertical diffusion on v
      ! ---------------------------

      DO ji = fs_2, fs_jpim1   ! vector opt.
         DO jj = 2, jpjm1   
            ikbvm1 = mbkv(ji,jj)
            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
!               zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbvm1)*zinv )
               zbfrv = bfrva(ji,jj)

            DO jk = 1, ikbvm1
               zcoef = -p2dt / fse3v(ji,jj,jk)
               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
               zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
            END DO

      ! Surface boudary conditions
            zwi(ji,jj,1) = 0._wp
            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)

      ! Bottom boudary conditions  
      ! Bottom boudary conditions  ! H. Liu, May, 2010
!            !commented out to be consisent with v3.2, h.liu
!            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0 * bfr_imp
            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * bfr_imp
            zws(ji,jj,ikbvm1) = 0._wp
            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 

      ! Matrix inversion
      !-----------------------------------------------------------------------
      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
      !
      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
      !        (        ...               )( ...  ) ( ...  )
      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
      !
      !   m is decomposed in the product of an upper and lower triangular
      !   matrix
      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
      !   The solution (after velocity) is in 2d array va
      !-----------------------------------------------------------------------

      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
            DO jk = 2, ikbvm1
               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
            END DO

      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 )
            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj))
            DO jk = 2, ikbvm1
               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
            END DO

      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1)

            DO jk = ikbvm1-1, 1, -1
               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
            END DO

         END DO
      END DO

      ! Normalization to obtain the general momentum trend va
      DO jk = 1, jpkm1
         DO jj = 2, jpjm1   
            DO ji = fs_2, fs_jpim1   ! vector opt.
               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
            END DO
         END DO
      END DO
      !
      IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array')
      !

   END SUBROUTINE dyn_zdf_imp

   !!==============================================================================
END MODULE dynzdf_imp