ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/LDF/ldfdyn_smag.F90

MODULE ldfdyn_smag
   !!======================================================================
   !!                     ***  MODULE  ldftrasmag  ***
   !! Ocean physics:  variable eddy induced velocity coefficients
   !!======================================================================
#if   defined key_dynldf_smag   &&   defined key_dynldf_c3d
   !!----------------------------------------------------------------------
   !!   'key_dynldf_smag'      and           smagorinsky  diffusivity
   !!   'key_dynldf_c3d'                    3D tracer lateral  mixing coef.
   !!----------------------------------------------------------------------
   !!   ldf_eiv      : compute the eddy induced velocity coefficients
   !!----------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE sbc_oce         ! surface boundary condition: ocean
   USE sbcrnf          ! river runoffs
   USE ldfdyn_oce      ! ocean tracer   lateral physics
   USE phycst          ! physical constants
   USE ldfslp          ! iso-neutral slopes
   USE in_out_manager  ! I/O manager
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE prtctl          ! Print control
   USE iom
   USE wrk_nemo
   IMPLICIT NONE
   PRIVATE

   !! * Routine accessibility
   PUBLIC ldf_dyn_smag               ! routine called by step.F90
   !!----------------------------------------------------------------------
   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
   !! $Id: ldfdyn_smag.F90 2355 2015-05-20 07:11:50Z ufla $
   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
   !!----------------------------------------------------------------------
   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------

CONTAINS





   !!----------------------------------------------------------------------
   !!                        ***  ldfdyn_smag.F90  ***
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
   !! $Id: ldfdyn_smag.F90 2355 2015-05-20 07:11:50Z ufla $ 
   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   'key_dynldf_smag'             3D lateral eddy viscosity coefficients
   !!----------------------------------------------------------------------

   SUBROUTINE ldf_dyn_smag( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE ldf_dyn_smag  ***
      !!                   
      !! ** Purpose :   initializations of the horizontal ocean physics
      !!
      !! ** Method  :   3D eddy viscosity coef. 
      !!    M.Griffies, R.Hallberg AMS, 2000
      !! for laplacian:
      !!   Asmag=(C/pi)^2*dx*dy sqrt(D^2), C=3-4
      !! for bilaplacian:
      !!   Bsmag=Asmag*dx*dy/8
      !!   D^2=(du/dx-dv/dy)^2+(dv/dx+du/dy)^2 for Cartesian coordinates
      !!  in general case du/dx ==> e2 d(u/e2)/dx;  du/dy ==> e1 d(u/e1)/dy; 
      !!                  dv/dx ==> e2 d(v/e2)/dx;  dv/dy ==> e1 d(v/e1)/dy
      !!
      !!       laplacian operator   : ahm1, ahm2 defined at T- and F-points
      !!                              ahm3, ahm4 never used
      !!       bilaplacian operator : ahm1, ahm2 never used
      !!                           :  ahm3, ahm4 defined at U- and V-points
      !!       explanation of the default is missingi
      !!  last modified : Maria Luneva, September 2011
      !!----------------------------------------------------------------------
      !! * Modules used
      !! ahm0 here is a background viscosity

      !! * Arguments

      !! * local variables

      INTEGER              :: kt                   ! timestep

      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
      REAL (wp):: zdeltat,zdeltaf,zdeltau,zdeltav  ! temporary scalars
      REAL (wp), POINTER, DIMENSION (:,:) ::   zux, zuy , zvx ,zvy, zue1, zue2, zve1, zve2 
      REAL (wp)::  zcmsmag_1, zcmsmag_2 , zcmsh


      !!----------------------------------------------------------------------

      CALL wrk_alloc( jpi,jpj,zux,zuy,zvx,zvy )
      CALL wrk_alloc( jpi,jpj,zue1,zue2,zve1,zve2 )


      IF(  kt == nit000 ) THEN


      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'ldf_dyn_smag : 3D lateral eddy viscosity coefficient'
      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
     
      ENDIF
      
      zcmsmag_1 = rn_cmsmag_1
      zcmsmag_2 = rn_cmsmag_2
      zcmsh     = rn_cmsh

   


      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operators
      ! =================                       whatever its orientation is)
      IF( ln_dynldf_lap ) THEN
         ! define ahm1 and ahm2 at the right grid point position
         ! (USER: modify ahm1 and ahm2 following your desiderata)
         
         DO jk=1,jpk
           zue2(:,:)=un(:,:,jk)/e2u(:,:)
           zve1(:,:)=vn(:,:,jk)/e1v(:,:)
           zue1(:,:)=un(:,:,jk)/e1u(:,:)
           zve2(:,:)=vn(:,:,jk)/e2v(:,:)

 
           DO jj=2,jpj
            DO ji=2,jpi
            zux(ji,jj)=(zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk) * zcmsh
            zvy(ji,jj)=(zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk) * zcmsh
            ENDDO
           ENDDO

           DO jj=1,jpjm1
            DO ji=1,jpim1
            zuy(ji,jj)=(zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
            zvx(ji,jj)=(zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
            ENDDO
           ENDDO
             
          DO jj=2,jpjm1
           DO ji=2,jpim1

            zdeltat=2._wp /(e1t(ji,jj)**(-2)+e2t(ji,jj)**(-2))
            zdeltaf=2._wp /(e1f(ji,jj)**(-2)+e2f(ji,jj)**(-2))   
            ahm1(ji,jj,jk)=(zcmsmag_1/rpi)**2*zdeltat*                                         &
                            sqrt( (zux(ji,jj)-zvy(ji,jj))**2+                                    &
                     0.0625_wp*(zuy(ji,jj)+zuy(ji,jj-1)+zuy(ji-1,jj)+zuy(ji-1,jj-1)+             &
                            zvx(ji,jj)+zvx(ji,jj-1)+zvx(ji-1,jj)+zvx(ji-1,jj-1))**2)

            ahm2(ji,jj,jk)=(zcmsmag_1/rpi)**2*zdeltaf*                                         &
                            sqrt( (zuy(ji,jj)+zvx(ji,jj))**2+                                    &
                     0.0625_wp*(zux(ji,jj)+zux(ji,jj+1)+zux(ji+1,jj)+zux(ji+1,jj+1)-             &
                             zvy(ji,jj)-zvy(ji,jj+1)-zvy(ji+1,jj)-zvy(ji+1,jj+1))**2)

            ahm1(ji,jj,jk)=MAX(ahm1(ji,jj,jk),rn_ahm_0_lap)
            ahm2(ji,jj,jk)=MAX(ahm2(ji,jj,jk),rn_ahm_0_lap)

! stability criteria or upper limit set from namelist
            ahm1(ji,jj,jk)=MIN(ahm1(ji,jj,jk),zdeltat / (16_wp*rdt),rn_ahm_m_lap)
            ahm2(ji,jj,jk)=MIN(ahm2(ji,jj,jk),zdeltaf / (16_wp*rdt),rn_ahm_m_lap)

            ENDDO
           ENDDO

         ENDDO ! jpk
            ahm1(:,:,jpk) = ahm1(:,:,jpkm1)
            ahm2(:,:,jpk) = ahm2(:,:,jpkm1)
            IF(lwp.and.kt==nit000) WRITE(numout,'(36x," ahm ", 7x)')
            DO jk = 1, jpk

               IF(lwp.and.kt==nit000) WRITE(numout,'(30x,E10.2,8x,i3)') ahm1(jpi/2,jpj/2,jk), jk
            END DO
      CALL lbc_lnk( ahm1, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
      CALL lbc_lnk( ahm2, 'F', 1. )   ! Lateral boundary conditions on ( ahtt )

      ENDIF    ! ln_dynldf_lap
      


      ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
      ! ================================  whatever its orientation is)
      ! Here: ahm is proportional to the cube of the maximum of the grid spacing
      !       in the to horizontal direction

      IF( ln_dynldf_bilap ) THEN
         DO jk=1,jpk
           zue2(:,:) = un(:,:,jk)/e2u(:,:)
           zve1(:,:) = vn(:,:,jk)/e1v(:,:)
           zue1(:,:) = un(:,:,jk)/e1u(:,:)
           zve2(:,:) = vn(:,:,jk)/e2v(:,:)


           DO jj=2,jpj
            DO ji=2,jpi
            zux(ji,jj) = (zue2(ji,jj)-zue2(ji-1,jj))/e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
            zvy(ji,jj) = (zve1(ji,jj)-zve1(ji,jj-1))/e2t(ji,jj)*e1t(ji,jj)*tmask(ji,jj,jk)
            ENDDO
           ENDDO

           DO jj=1,jpjm1
            DO ji=1,jpim1
            zuy(ji,jj) = (zue1(ji,jj+1)-zue1(ji,jj))/e2f(ji,jj)*e1f(ji,jj)*fmask(ji,jj,jk)
            zvx(ji,jj) = (zve2(ji+1,jj)-zve2(ji,jj))/e1f(ji,jj)*e2f(ji,jj)*fmask(ji,jj,jk)
            ENDDO
           ENDDO


          DO jj=2,jpjm1
           DO ji=2,jpim1
            zdeltau = 2._wp/(e1u(ji,jj)**(-2)+e2u(ji,jj)**(-2))
            zdeltav = 2._wp/(e1v(ji,jj)**(-2)+e2v(ji,jj)**(-2))

            ahm3(ji,jj,jk) = -(zcmsmag_2/rpi)**2/8.0_wp*zdeltau**2*                          &

                         sqrt(0.25_wp*(zux(ji,jj)+zux(ji+1,jj)-zvy(ji,jj)-zvy(ji+1,jj))**2+    &
                              0.25_wp*(zuy(ji,jj)+zuy(ji,jj-1)+zvx(ji,jj)+zvx(ji,jj-1))**2)

            ahm4(ji,jj,jk) = -(zcmsmag_2/rpi)**2/8.0_wp*zdeltav**2*                           &

                         sqrt(0.25_wp*(zux(ji,jj)+zux(ji,jj+1)-zvy(ji,jj)-zvy(ji,jj+1))**2+    &
                              0.25_wp*(zuy(ji,jj)+zuy(ji-1,jj)+zvx(ji-1,jj)+zvx(ji,jj))**2)

            ahm3(ji,jj,jk) = MIN (rn_ahm_0_blp , ahm3(ji,jj,jk) )
            ahm4(ji,jj,jk) = MIN (rn_ahm_0_blp , ahm4(ji,jj,jk) ) 

! stability criteria or upper limit set in namelist

            ahm3(ji,jj,jk) = MAX( ahm3(ji,jj,jk),-zdeltau**2/( 128._wp*rdt ),rn_ahm_m_blp )
            ahm4(ji,jj,jk) = MAX( ahm4(ji,jj,jk),-zdeltav**2/( 128._wp*rdt ),rn_ahm_m_blp )
            

            ENDDO
           ENDDO

         ENDDO
            ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
            ahm4(:,:,jpk) = ahm4(:,:,jpkm1)

      DO jk = 1, jpk
      IF(  kt == nit000 ) THEN

               IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(jpi/2,jpj/2,jk), jk
      ENDIF   
      END DO
      CALL lbc_lnk( ahm3, 'U', 1. )   ! Lateral boundary conditions
      CALL lbc_lnk( ahm4, 'V', 1. )
   ENDIF

      CALL wrk_dealloc( jpi,jpj,zux,zuy,zvx,zvy )
      CALL wrk_dealloc( jpi,jpj,zue1,zue2,zve1,zve2 )
   !    zumax = MAXVAL( ABS( ahm3(:,:,:) ) )                ! slower than the following loop on NEC SX5
      zdeltat = 0._wp
   If(ln_dynldf_lap)THEN
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               zdeltat = MAX(zdeltat,ABS(ahm1(ji,jj,jk)),ABS(ahm2(ji,jj,jk)) )
          END DO
        END DO
      END DO
      IF( lk_mpp )   CALL mpp_max( zdeltat )                 ! max over the global domain
      !
      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,'dynlap:  abs(ahm) max: ', zdeltat
    ENDIF
    If(ln_dynldf_bilap)THEN
     zdeltat = 0._wp
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               zdeltat = MAX(zdeltat,ABS(ahm3(ji,jj,jk)),ABS(ahm3(ji,jj,jk)) )
          END DO
        END DO
      END DO
      IF( lk_mpp )   CALL mpp_max( zdeltat )                 ! max over the global domain
      !
      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,'dyn_bilap abs(ahm) max: ', zdeltat
      !
   ENDIF
      !

END SUBROUTINE ldf_dyn_smag
#else
   !!----------------------------------------------------------------------
   !!   Default option                                         Dummy module
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE ldf_dyn_smag( kt )       ! Empty routine
      WRITE(*,*) 'ldf_dyn_smag: You should not have seen this print! error? check keys ldf:c3d+smag', kt
   END SUBROUTINE ldf_dyn_smag
#endif

   END MODULE ldfdyn_smag