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

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 solisl          ! ???
   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) 
   !! $Header$ 
   !! 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:
      !!       * default option: barotropic stream function system
      !!         and islands initialization (if lk_isl=T)
      !!       * lk_dynspg_flt = T : transport divergence system. No specific
      !!         treatment of islands.
      !!      
      !! ** Method :
      !!       - Compute the local depth of the water column at u- and v-point
      !!      (lk_dynspg_flt = T) or its inverse (lk_dynspg_rl = T).
      !!      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.
      !!       - Construct the matrix of the elliptic system by a call to
      !!      solmat.F routine.
      !!       - island (if lk_isl=T)
      !!            isl_dom: find islands from the bathymetry file
      !!            isl_bsf: compute the island barotropic stream function
      !!            isl_mat: compute the inverse island matrix
      !!            set mbathy to the number of non-zero w-levels of a water
      !!            column (the minimum value of mbathy is 2):
      !!                  mbathy = min( mbathy, 1 ) + 1
      !!
      !! ** Action : - hur, hvr : masked inverse of the local depth at
      !!                                u- and v-point. (lk_dynspg_rl = T)
      !!             - hu, hv   : masked local depth at u- and v- points
      !!                                (lk_dynspg_flt = T)
      !!             - 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

      !! * Local declarations
      INTEGER :: ji, jj   ! dummy loop indices
      CHARACTER(len=80) :: clname

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

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'solver_init : solver to compute the surface pressure gradient'
      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'

      ! open elliptic solver statistics file
      clname = 'solver.stat'
      CALL ctlopn( numsol, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
                   1, numout, lwp, 1 )


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

#if defined key_feti
      ! FETI algorithm, we force nsolv at 3
      nsolv = 3
#endif

      ! 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
         IF(lk_isl) WRITE(numout,*) '             absolute precision stream fct    epsisl   = ', epsisl
         IF(lk_isl) WRITE(numout,*) '             maximum pcg iterations island    nmisl    = ', nmisl
         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'
         IF( lk_isl ) CALL ctl_stop( ' key_islands inconsistent with key_dynspg_flt' )
    
      ELSEIF( lk_dynspg_rl ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '          Rigid lid formulation'
      ELSE
         CALL ctl_stop( '          Choose only one surface pressure gradient calculation: filtered or rigid-lid',   &
              &         '          Should not call this routine if dynspg_exp or dynspg_ts has been chosen' )
      ENDIF
      IF( lk_dynspg_flt .AND. lk_dynspg_rl ) &
         CALL ctl_stop( '          Chose between free surface or rigid-lid, not both' )

      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 is used'
         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) &
             CALL ctl_stop( ' jpr2di and jpr2dj should be equal to zero' )

      CASE ( 3 )                ! FETI solver
         IF(lwp) WRITE(numout,*) '          the FETI solver is used'
         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) &
              CALL ctl_stop( ' jpr2di and jpr2dj should be equal to zero' )
           
         IF( .NOT.lk_mpp ) THEN
            CALL ctl_stop( ' The FETI algorithm is used only with the key_mpp_... option' )
         ELSE
            IF( jpnij == 1 ) THEN
               CALL ctl_stop( ' The FETI algorithm needs more than one processor' )
            ENDIF
         ENDIF
         
      CASE ( 4 )                ! successive-over-relaxation solver with extra outer halo
         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 .AND. nsolv /= 2 ) THEN
         CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require the SOR solver: nsolv = 2' )
      END IF

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

      IF( lk_dynspg_rl ) THEN       ! rigid-lid
         IF( lk_mpp ) THEN
            c_solver_pt = 'G'   ! G= F with special staff ??? which one?
         ELSE
            c_solver_pt = 'F'
         ENDIF
      ELSE                          ! free surface T-point
         IF( lk_mpp ) THEN
            c_solver_pt = 'S'   ! S=T with special staff ??? which one?
         ELSE
            c_solver_pt = 'T'
         ENDIF
      ENDIF


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

      CALL sol_mat( kt )


      IF( lk_isl ) THEN
      
         ! Islands in the domain
         ! ---------------------

         IF ( jpisl == 0 ) THEN
             WRITE(ctmp1,*) ' bad islands parameter jpisl =', jpisl
             CALL ctl_stop( ctmp1 )
         ENDIF

         ! open Island streamfunction statistic file
         CALL ctlopn( numisp, 'islands.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
            &         1     , numout        , lwp      , 1                         )
   
         CALL isl_dom       ! Island identification

         CALL isl_bsf       ! Island barotropic stream function

         CALL isl_mat       ! Comput and invert the island matrix

         ! mbathy set to the number of w-level (minimum value 2)
         DO jj = 1, jpj
            DO ji = 1, jpi
               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
            END DO
         END DO

      ENDIF

   END SUBROUTINE solver_init

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