source: trunk/NEMO/OPA_SRC/SOL/solver.F90 @ 1596

MODULE solver
   !!======================================================================
   !!                     ***  MODULE  solver  ***
   !! Ocean solver :  initialization of ocean solver
   !!=====================================================================

   !!----------------------------------------------------------------------
   !!   solver_init: solver initialization
   !!----------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and tracers variables
   USE dom_oce         ! ocean space and time domain variables 
   USE zdf_oce         ! ocean vertical physics variables
   USE sol_oce         ! solver variables
   USE solmat          ! ???
   USE obc_oce         ! Lateral open boundary condition
   USE in_out_manager  ! I/O manager
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE lib_mpp
   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient

   IMPLICIT NONE

   !!----------------------------------------------------------------------
   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
   !! $Id$ 
   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE solver_init( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE solver_init  ***
      !!                   
      !! ** Purpose :   Initialization for the solver of the elliptic equation:
      !!       * lk_dynspg_flt = T : transport divergence system. 
      !!      
      !! ** Method :
      !!       - Compute the local depth of the water column at u- and v-point
      !!      The local depth of the water column is computed by summing 
      !!      the vertical scale factors. For its inverse, the thickness of
      !!      the first model level is imposed as lower bound. The inverse of
      !!      this depth is THEN taken and masked, so that the inverse of the
      !!      local depth is zero when the local depth is zero.
      !!
      !! ** Action : - hur, hvr : masked inverse of the local depth at
      !!                                u- and v-point.
      !!             - hu, hv   : masked local depth at u- and v- points
      !!             - c_solver_pt : nature of the gridpoint at which the
      !!                                solver is applied
      !! References :
      !!      Jensen, 1986: adv. phys. oceanogr. num. mod.,ed. o brien,87-110.
      !!      Madec & Marti, 1990: internal rep. LODYC, 90/03., 29pp.
      !!
      !! History :
      !!        !  90-10  (G. Madec)  Original code           
      !!        !  93-02  (O. Marti)                         
      !!        !  97-02  (G. Madec)  local depth inverse computation
      !!        !  98-10  (G. Roullet, G. Madec)  free surface 
      !!   9.0  !  03-07  (G. Madec)  free form, F90
      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT(in) :: kt

      NAMELIST/namsol/ nsolv, nsol_arp, nmin, nmax, nmod, eps, resmax, sor, rnu
      !!----------------------------------------------------------------------

      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'solver_init : solver to compute the surface pressure gradient'
         WRITE(numout,*) '~~~~~~~~~~~'
         
         ! open elliptic solver statistics file
         CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
      ENDIF


      ! 0. Define the solver parameters
      !    ----------------------------
      ! Namelist namsol : elliptic solver / free surface
      REWIND( numnam )
      READ  ( numnam, namsol )

      ! 0. Parameter control and print
      !    ---------------------------

      ! Control print
      IF(lwp) WRITE(numout,*) '          Namelist namsol : set solver parameters'

      IF(lwp) THEN
         WRITE(numout,*) '             type of elliptic solver            nsolv    = ', nsolv
         WRITE(numout,*) '             absolute/relative (0/1) precision  nsol_arp = ', nsol_arp
         WRITE(numout,*) '             minimum iterations for solver      nmin     = ', nmin
         WRITE(numout,*) '             maximum iterations for solver      nmax     = ', nmax
         WRITE(numout,*) '             frequency for test                 nmod     = ', nmod
         WRITE(numout,*) '             absolute precision of solver       eps      = ', eps
         WRITE(numout,*) '             absolute precision for SOR solver  resmax   = ', resmax
         WRITE(numout,*) '             optimal coefficient of sor         sor      = ', sor
         WRITE(numout,*) '             free surface parameter         rnu    = ', rnu
         WRITE(numout,*)
      ENDIF

      IF( lk_dynspg_flt ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '          free surface formulation'
      ELSE
         CALL ctl_stop( '          Choose only one surface pressure gradient calculation: filtered ',   &
              &         '          Should not call this routine if dynspg_exp or dynspg_ts has been chosen' )
      ENDIF

      SELECT CASE ( nsolv )

      CASE ( 1 )                ! preconditioned conjugate gradient solver
         IF(lwp) WRITE(numout,*) '          a preconditioned conjugate gradient solver is used'
         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) &
            CALL ctl_stop( ' jpr2di and jpr2dj should be equal to zero' )

      CASE ( 2 )                ! successive-over-relaxation solver
         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver with extra outer halo is used'
         IF(lwp) WRITE(numout,*) '          with jpr2di =', jpr2di, ' and  jpr2dj =', jpr2dj
         IF( .NOT. lk_mpp .AND. jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
             CALL ctl_stop( ' jpr2di and jpr2dj are not equal to zero',   &
             &              ' In this case this algorithm should be used only with the key_mpp_... option' )
         ELSE
            IF( ( ( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) .OR. ( jpni /= 1 ) ) &
              &  .AND. ( jpr2di /= jpr2dj ) ) CALL ctl_stop( '          jpr2di should be equal to jpr2dj' )
         ENDIF

      CASE DEFAULT
         WRITE(ctmp1,*) '          bad flag value for nsolv = ', nsolv
         CALL ctl_stop( ctmp1 )
         
      END SELECT

      IF( nbit_cmp == 1 ) THEN
         IF( nsolv /= 2 ) THEN
            CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require the SOR solver: nsolv = 2' )
         ELSE IF( MAX( jpr2di, jpr2dj ) > 0 ) THEN
            CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require jpr2di = jpr2dj = 0' )
         END IF 
      END IF

      ! Grid-point at which the solver is applied
      ! -----------------------------------------

      IF( lk_mpp ) THEN
         c_solver_pt = 'S'   ! S=T with special staff ??? which one?
      ELSE
         c_solver_pt = 'T'
      ENDIF

      ! Construction of the elliptic system matrix
      ! ------------------------------------------

      CALL sol_mat( kt )
      !
   END SUBROUTINE solver_init

   !!======================================================================
END MODULE solver