source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d.F90 @ 4596

MODULE ldfc1d
   !!======================================================================
   !!                    ***  MODULE  ldfc1d  ***
   !! Ocean physics:  profile of lateral eddy coefficients 
   !!=====================================================================
   !! History :  3.7  ! 2013-12  (G. Madec)  restructuration/simplification of aht/aeiv specification,
   !!                 !                      add velocity dependent coefficient and optional read in file
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   ldf_c1d       : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m) 
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers
   USE dom_oce        ! ocean space and time domain
   USE phycst         ! physical constants
   !
   USE in_out_manager ! I/O manager
   USE lib_mpp        ! distribued memory computing library
   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)

   IMPLICIT NONE
   PRIVATE

   PUBLIC   ldf_c1d   ! called by ldftra and ldfdyn modules

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

   SUBROUTINE ldf_c1d( cd_type, prat, pahs1, pahs2, pah1, pah2 )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE ldftra_c1d  ***
      !!              
      !! ** Purpose :   1D eddy diffusivity/viscosity coefficients
      !!
      !! ** Method  :   1D eddy diffusivity coefficients F( depth )
      !!                Reduction by prat from surface to bottom 
      !!                hyperbolic tangent profile with inflection point 
      !!                at zh=500m and a width of zw=200m
      !!
      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
      !!             DYN      pah1, pah2 defined at T- and F-points
      !!----------------------------------------------------------------------
      CHARACTER(len=2)                , INTENT(in   ) ::   cd_type        ! DYNamique or TRAcers
      REAL(wp)                        , INTENT(in   ) ::   prat           ! ratio surface/deep values           [-]
      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pahs1, pahs2   ! surface value of eddy coefficient   [m2/s]
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pah1 , pah2    ! eddy coefficient                    [m2/s]
      !
      INTEGER  ::   ji, jj, jk      ! dummy loop indices
      REAL(wp) ::   zh, zc, zdep1   ! local scalars
      REAL(wp) ::   zw    , zdep2   !   -      -
      !!----------------------------------------------------------------------

      ! initialization of the profile
      zh =  500._wp              ! depth    of the inflection point [m]
      zw =  1._wp / 200._wp      ! width^-1     -        -      -   [1/m]
      !                          ! associated coefficient           [-]
      zc = ( 1._wp - prat ) / ( 1._wp + TANH( zh * zw) )
      !
      !
      SELECT CASE( cd_type )        ! point of calculation
      !
      CASE( 'DYN' )                     ! T- and F-points
         DO jk = 1, jpk                      ! pah1 at T-point
            pah1(:,:,jk) = pahs1(:,:) * (  prat + zc * ( 1._wp + TANH( - ( fsdept(:,:,jk) - zh ) * zw) )  ) * tmask(:,:,jk)
         END DO
         DO jk = 1, jpk                      ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
            DO jj = 1, jpjm1
               DO ji = 1, fs_jpim1
                  zdep2 = (  fsdept(ji,jj+1,jk) + fsdept(ji+1,jj+1,jk)   &
                     &     + fsdept(ji,jj  ,jk) + fsdept(ji+1,jj  ,jk)  ) * 0.25_wp
                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * fmask(ji,jj,jk)
               END DO
            END DO
         END DO
         CALL lbc_lnk( pah2, 'F', 1. )   ! Lateral boundary conditions
         !
      CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.)
         DO jk = 1, jpk
            DO jj = 1, jpjm1
               DO ji = 1, fs_jpim1
                  zdep1 = (  fsdept(ji,jj,jk) + fsdept(ji+1,jj,jk)  ) * 0.5_wp
                  zdep2 = (  fsdept(ji,jj,jk) + fsdept(ji,jj+1,jk)  ) * 0.5_wp
                  pah1(ji,jj,jk) = pahs1(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  ) * umask(ji,jj,jk)
                  pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * vmask(ji,jj,jk)
               END DO
            END DO
         END DO
         CALL lbc_lnk( pah1, 'U', 1. )   ! Lateral boundary conditions
         CALL lbc_lnk( pah2, 'V', 1. )   
         !
      END SELECT
      !
   END SUBROUTINE ldf_c1d

   !!======================================================================
END MODULE ldfc1d