source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90 @ 7698

MODULE trazdf
   !!==============================================================================
   !!                 ***  MODULE  trazdf  ***
   !! Ocean active tracers:  vertical component of the tracer mixing trend
   !!==============================================================================
   !! History :  1.0  ! 2005-11  (G. Madec)  Original code
   !!            3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   tra_zdf       : Update the tracer trend with the vertical diffusion
   !!   tra_zdf_init  : initialisation of the computation
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers variables
   USE dom_oce        ! ocean space and time domain variables 
   USE domvvl         ! variable volume
   USE phycst         ! physical constant
   USE zdf_oce        ! ocean vertical physics variables
   USE sbc_oce        ! surface boundary condition: ocean
   USE ldftra         ! lateral diffusion: eddy diffusivity
   USE ldfslp         ! lateral diffusion: iso-neutral slope 
   USE trazdf_exp     ! vertical diffusion: explicit (tra_zdf_exp routine)
   USE trazdf_imp     ! vertical diffusion: implicit (tra_zdf_imp routine)
   USE trd_oce        ! trends: ocean variables
   USE trdtra         ! trends: tracer trend manager
   !
   USE in_out_manager ! I/O manager
   USE prtctl         ! Print control
   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
   USE lib_mpp        ! MPP library
   USE wrk_nemo       ! Memory allocation
   USE timing         ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC   tra_zdf        ! routine called by step.F90
   PUBLIC   tra_zdf_init   ! routine called by nemogcm.F90

   INTEGER ::   nzdf = 0   ! type vertical diffusion algorithm used (defined from ln_zdf...  namlist logicals)

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

   SUBROUTINE tra_zdf( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE tra_zdf  ***
      !!
      !! ** Purpose :   compute the vertical ocean tracer physics.
      !!---------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
      !
      INTEGER  ::   jk, jj, ji           ! Dummy loop indices
      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
      !!---------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('tra_zdf')
      !
      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
         r2dt =  rdt                          ! = rdt (restarting with Euler time stepping)
      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
         r2dt = 2. * rdt                      ! = 2 rdt (leapfrog)
      ENDIF
      !
      IF( l_trdtra )   THEN                    !* Save ta and sa trends
         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
         DO jk = 1, jpk
            DO jj = 1, jpj
               DO ji = 1, jpi
                  ztrdt(ji,jj,jk) = tsa(ji,jj,jk,jp_tem)
                  ztrds(ji,jj,jk) = tsa(ji,jj,jk,jp_sal)
               END DO
            END DO
         END DO
      ENDIF
      !
      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme 
      END SELECT
!!gm WHY here !   and I don't like that !
      ! DRAKKAR SSS control {
      ! JMM avoid negative salinities near river outlet ! Ugly fix
      ! JMM : restore negative salinities to small salinities:
!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               IF( tsa(ji,jj,jk,jp_sal) < 0._wp )   tsa(ji,jj,jk,jp_sal) = 0.1_wp
            END DO
         END DO
      END DO
!!gm

      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
         DO jk = 1, jpkm1
            DO jj = 1, jpj
               DO ji = 1, jpi
                  ztrdt(ji,jj,jk) = ( ( tsa(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) / r2dt ) - ztrdt(ji,jj,jk)
                  ztrds(ji,jj,jk) = ( ( tsa(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) / r2dt ) - ztrds(ji,jj,jk)
               END DO
            END DO
         END DO
!!gm this should be moved in trdtra.F90 and done on all trends
         CALL lbc_lnk( ztrdt, 'T', 1. )
         CALL lbc_lnk( ztrds, 'T', 1. )
!!gm
         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
      ENDIF
      !                                          ! print mean trends (used for debugging)
      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
      !
      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf')
      !
   END SUBROUTINE tra_zdf


   SUBROUTINE tra_zdf_init
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE tra_zdf_init  ***
      !!
      !! ** Purpose :   Choose the vertical mixing scheme
      !!
      !! ** Method  :   Set nzdf from ln_zdfexp
      !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T)
      !!           = 1   implicit (euler backward) scheme (ln_zdfexp=F)
      !!      NB: rotation of lateral mixing operator or TKE & GLS schemes,
      !!          an implicit scheme is required.
      !!----------------------------------------------------------------------
      USE zdftke
      USE zdfgls
      !!----------------------------------------------------------------------
      !
      ! Choice from ln_zdfexp already read in namelist in zdfini module
      IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
      ELSE                   ;   nzdf = 1           ! use implicit scheme
      ENDIF
      !
      ! Force implicit schemes
      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE, or GLS physics
      IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
      IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
      IF( ln_zdfexp .AND. nzdf == 1 )   CALL ctl_stop( 'tra_zdf : If using the rotation of lateral mixing operator',   &
            &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' )
            !
      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme'
         WRITE(numout,*) '~~~~~~~~~~~'
         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
      ENDIF
      !
   END SUBROUTINE tra_zdf_init

   !!==============================================================================
END MODULE trazdf