source: trunk/NEMO/OPA_SRC/TRA/trabbc.F90 @ 719

MODULE trabbc
   !!==============================================================================
   !!                       ***  MODULE  trabbc  ***
   !! Ocean active tracers:  bottom boundary condition
   !!==============================================================================
   !! History :  8.1  ! 99-10 (G. Madec)  original code
   !!            8.5  ! 02-08 (G. Madec)  free form + modules
   !!            8.5  ! 02-11 (A. Bozec)  tra_bbc_init: original code
   !!----------------------------------------------------------------------
#if   defined key_trabbc   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_trabbc'                                  geothermal heat flux
   !!----------------------------------------------------------------------
   !!   tra_bbc      : update the tracer trend at ocean bottom 
   !!   tra_bbc_init : initialization of geothermal heat flux trend
   !!----------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and active tracers
   USE dom_oce         ! ocean space and time domain
   USE phycst          ! physical constants
   USE trdmod          ! ocean trends 
   USE trdmod_oce      ! ocean variables trends
   USE in_out_manager  ! I/O manager
   USE prtctl          ! Print control

   IMPLICIT NONE
   PRIVATE

   PUBLIC tra_bbc          ! routine called by step.F90

   !! to be transfert in the namelist ???!   
   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbc = .TRUE.   !: bbc flag

   !!* Namelist nambbc: bottom boundary condition
   INTEGER  ::   ngeo_flux       = 1            ! Geothermal flux (0:no flux, 1:constant flux, 2:read in file )
   REAL(wp) ::   ngeo_flux_const = 86.4e-3      ! Constant value of geothermal heat flux

   INTEGER , DIMENSION(jpi,jpj) ::   nbotlevt   ! ocean bottom level index at T-pt
   REAL(wp), DIMENSION(jpi,jpj) ::   qgh_trd0   ! geothermal heating trend
 
   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !!  OPA 9.0 , LOCEAN-IPSL (2006) 
   !! $Header$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE tra_bbc( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE tra_bbc  ***
      !!
      !! ** Purpose :   Compute the bottom boundary contition on temperature 
      !!      associated with geothermal heating and add it to the general
      !!      trend of temperature equations.
      !!
      !! ** Method  :   The geothermal heat flux set to its constant value of 
      !!       86.4 mW/m2 (Stein and Stein 1992, Huang 1999).
      !!       The temperature trend associated to this heat flux through the
      !!       ocean bottom can be computed once and is added to the temperature
      !!       trend juste above the bottom at each time step:
      !!            ta = ta + Qsf / (rau0 rcp e3T) for k= mbathy -1
      !!       Where Qsf is the geothermal heat flux.
      !!
      !! ** Action  : - update the temperature trends (ta) with the trend of
      !!                the ocean bottom boundary condition
      !!
      !! References : Stein, C. A., and S. Stein, 1992, Nature, 359, 123-129.
      !!----------------------------------------------------------------------
      USE oce, ONLY :   ztrdt => ua   ! use ua as 3D workspace   
      USE oce, ONLY :   ztrds => va   ! use va as 3D workspace   
      !!
      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
      !!
#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
      INTEGER ::   ji       ! dummy loop indices
#else
      INTEGER ::   ji, jj   ! dummy loop indices
#endif
      REAL(wp) ::   zqgh_trd  ! geothermal heat flux trend
      !!----------------------------------------------------------------------

      IF( kt == nit000 )   CALL tra_bbc_init      ! Initialization

      IF( l_trdtra )   THEN         ! Save ta and sa trends
         ztrdt(:,:,:) = ta(:,:,:) 
         ztrds(:,:,:) = 0.e0
      ENDIF

      ! Add the geothermal heat flux trend on temperature

      SELECT CASE ( ngeo_flux )
      !
      CASE ( 1:2 )                !  geothermal heat flux
#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
            zqgh_trd = ro0cpr * qgh_trd0(ji,1) / fse3t(ji,1,nbotlevt(ji,1) )
            ta(ji,1,nbotlevt(ji,1)) = ta(ji,1,nbotlevt(ji,1)) + zqgh_trd
         END DO
#else
         DO jj = 2, jpjm1
            DO ji = 2, jpim1
               zqgh_trd = ro0cpr * qgh_trd0(ji,jj) / fse3t(ji,jj,nbotlevt(ji,jj))
               ta(ji,jj,nbotlevt(ji,jj)) = ta(ji,jj,nbotlevt(ji,jj)) + zqgh_trd
            END DO
         END DO
#endif
      END SELECT

      IF( l_trdtra ) THEN        ! Save the geothermal heat flux trend for diagnostics
         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
         CALL trd_mod( ztrdt, ztrds, jptra_trd_bbc, 'TRA', kt )
      ENDIF
      !
      IF(ln_ctl)   CALL prt_ctl(tab3d_1=ta, clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta')
      !
   END SUBROUTINE tra_bbc


   SUBROUTINE tra_bbc_init
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE tra_bbc_init  ***
      !!
      !! ** Purpose :   Compute once for all the trend associated with geo-
      !!      thermal heating that will be applied at each time step at the
      !!      bottom ocean level
      !!
      !! ** Method  :   Read the nambbc namelist and check the parameters.
      !!      called at the first time step (nit000)
      !!
      !! ** Input   : - Namlist nambbc
      !!              - NetCDF file  : geothermal_heating.nc ( if necessary )
      !!
      !! ** Action  : - read/fix the geothermal heat qgh_trd0
      !!              - compute the bottom ocean level nbotlevt
      !!----------------------------------------------------------------------
      USE iom
      !!
      INTEGER  ::   ji, jj              ! dummy loop indices
      INTEGER  ::   inum                ! temporary logical unit

      NAMELIST/nambbc/ngeo_flux, ngeo_flux_const 
      !!----------------------------------------------------------------------

      REWIND ( numnam )              ! Read Namelist nambbc : bottom momentum boundary condition
      READ   ( numnam, nambbc )

      !                              ! Control print
      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'tra_bbc : tempearture Bottom Boundary Condition (bbc)'
      IF(lwp) WRITE(numout,*) '~~~~~~~   Geothermal heatflux'
      IF(lwp) WRITE(numout,*) '          Namelist nambbc : set bbc parameters'
      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) '             Geothermal flux           ngeo_flux       = ', ngeo_flux
      IF(lwp) WRITE(numout,*) '             Constant geothermal flux  ngeo_flux_const = ', ngeo_flux_const
      IF(lwp) WRITE(numout,*)

      !                              ! level of the ocean bottom at T-point
      DO jj = 1, jpj
         DO ji = 1, jpi
            nbotlevt(ji,jj) = MAX( mbathy(ji,jj)-1, 1 )
         END DO
      END DO

      SELECT CASE ( ngeo_flux )      ! initialization of geothermal heat flux
      !
      CASE ( 0 )                ! no geothermal heat flux
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '             *** no geothermal heat flux'
         !
      CASE ( 1 )                ! constant flux
         IF(lwp) WRITE(numout,*) '             *** constant heat flux  =   ', ngeo_flux_const
         qgh_trd0(:,:) = ngeo_flux_const
         !
      CASE ( 2 )                ! variable geothermal heat flux
         ! read the geothermal fluxes in mW/m2
         !
         IF(lwp) WRITE(numout,*) '             *** variable geothermal heat flux'
         CALL iom_open ( 'geothermal_heating.nc', inum )
         CALL iom_get ( inum, jpdom_data, 'heatflow', qgh_trd0 )
         CALL iom_close (inum)
         !
         qgh_trd0(:,:) = qgh_trd0(:,:) * 1.e-3 ! conversion in W/m2
         !
      CASE DEFAULT
         WRITE(ctmp1,*) '     bad flag value for ngeo_flux = ', ngeo_flux
         CALL ctl_stop( ctmp1 )
         !
      END SELECT


   END SUBROUTINE tra_bbc_init

#else
   !!----------------------------------------------------------------------
   !!   Default option                                         Empty module
   !!----------------------------------------------------------------------
   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbc = .FALSE.  !: bbc flag
CONTAINS
   SUBROUTINE tra_bbc( kt )           ! Empty routine
      WRITE(*,*) 'tra_bbc: You should not have seen this print! error?', kt
   END SUBROUTINE tra_bbc
#endif

   !!======================================================================
END MODULE trabbc