New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
solver.F90 in trunk/NEMO/OPA_SRC/SOL – NEMO

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

Last change on this file since 100 was 79, checked in by opalod, 20 years ago

CT : UPDATE053 : Use logical key "lk_isl" instead of "l_isl"
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.2 KB
RevLine 
[3]1MODULE solver
2   !!======================================================================
3   !!                     ***  MODULE  solver  ***
4   !! Ocean solver :  initialization of ocean solver
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   solver_init: solver initialization
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE oce             ! ocean dynamics and tracers variables
12   USE dom_oce         ! ocean space and time domain variables
13   USE zdf_oce         ! ocean vertical physics variables
14   USE sol_oce         ! solver variables
15   USE solmat          ! ???
16   USE solisl          ! ???
17   USE obc_oce         ! Lateral open boundary condition
18   USE in_out_manager  ! I/O manager
19   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[16]20   USE lib_mpp
21   USE dynspg_rl       
22   USE dynspg_fsc     
[79]23   USE dynspg_fsc_atsk     
[3]24
25   IMPLICIT NONE
26
27   !!----------------------------------------------------------------------
28   !!   OPA 9.0 , LODYC-IPSL  (2003)
29   !!----------------------------------------------------------------------
30
31CONTAINS
32
33   SUBROUTINE solver_init
34      !!----------------------------------------------------------------------
35      !!                  ***  ROUTINE solver_init  ***
36      !!                   
37      !! ** Purpose :   Initialization for the solver of the elliptic equation:
38      !!       * default option: barotropic stream function system
[79]39      !!         and islands initialization (if lk_isl=T)
[16]40      !!       * lk_dynspg_fsc = T : transport divergence system. No specific
[3]41      !!         treatment of islands.
42      !!     
43      !! ** Method :
44      !!       - Compute the local depth of the water column at u- and v-point
[16]45      !!      (lk_dynspg_fsc = T) or its inverse (lk_dynspg_rl = T).
[3]46      !!      The local depth of the water column is computed by summing
47      !!      the vertical scale factors. For its inverse, the thickness of
48      !!      the first model level is imposed as lower bound. The inverse of
49      !!      this depth is THEN taken and masked, so that the inverse of the
50      !!      local depth is zero when the local depth is zero.
51      !!       - Construct the matrix of the elliptic system by a call to
52      !!      solmat.F routine.
[79]53      !!       - island (if lk_isl=T)
[3]54      !!            isl_dom: find islands from the bathymetry file
55      !!            isl_bsf: compute the island barotropic stream function
56      !!            isl_mat: compute the inverse island matrix
57      !!            set mbathy to the number of non-zero w-levels of a water
58      !!            column (the minimum value of mbathy is 2):
59      !!                  mbathy = min( mbathy, 1 ) + 1
60      !!
61      !! ** Action : - hur, hvr : masked inverse of the local depth at
[16]62      !!                                u- and v-point. (lk_dynspg_rl = T)
[3]63      !!             - hu, hv   : masked local depth at u- and v- points
[16]64      !!                                (lk_dynspg_fsc = T)
65      !!             - c_solver_pt : nature of the gridpoint at which the
66      !!                                solver is applied
[3]67      !! References :
68      !!      Jensen, 1986: adv. phys. oceanogr. num. mod.,ed. o brien,87-110.
69      !!      Madec & Marti, 1990: internal rep. LODYC, 90/03., 29pp.
70      !!
71      !! History :
72      !!        !  90-10  (G. Madec)  Original code           
73      !!        !  93-02  (O. Marti)                         
74      !!        !  97-02  (G. Madec)  local depth inverse computation
75      !!        !  98-10  (G. Roullet, G. Madec)  free surface
76      !!   9.0  !  03-07  (G. Madec)  free form, F90
77      !!----------------------------------------------------------------------
78      !! * Local declarations
79      INTEGER :: ji, jj   ! dummy loop indices
80
81      NAMELIST/namsol/ nsolv, nmax, eps, sor, epsisl, nmisl, rnu
82      !!----------------------------------------------------------------------
83
84      IF(lwp) WRITE(numout,*)
[79]85      IF(lwp) WRITE(numout,*) 'solver_init : solver to compute the surface pressure gradient'
86      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[3]87
88      ! open elliptic solver statistics file
89      CALL ctlopn( numsol, 'solver.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
90                   1, numout, lwp, 1 )
91
92
93      ! 0. Define the solver parameters
94      !    ----------------------------
95      ! Namelist namsol : elliptic solver / islands / free surface
96      REWIND( numnam )
97      READ  ( numnam, namsol )
98
99#if defined key_feti
100      ! FETI algorithm, we force nsolv at 3
101      nsolv = 3
102#endif
103
104
105      ! 0. Parameter control and print
106      !    ---------------------------
107
108      ! Control print
109      IF(lwp) WRITE(numout,*) '          Namelist namsol : set solver parameters'
110
111      IF(lwp) THEN
112         WRITE(numout,*) '             type of elliptic solver        nsolv  = ', nsolv
113         WRITE(numout,*) '             maximum iterations for solver  nmax   = ', nmax
114         WRITE(numout,*) '             absolute precision of solver   eps    = ', eps
115         WRITE(numout,*) '             optimal coefficient of sor     sor    = ', sor
[79]116         IF(lk_isl) WRITE(numout,*) '             absolute precision stream fct  epsisl = ', epsisl
117         IF(lk_isl) WRITE(numout,*) '             maximum pcg iterations island  nmisl  = ', nmisl
[3]118         WRITE(numout,*) '             free surface parameter         rnu    = ', rnu
119         WRITE(numout,*)
120      ENDIF
121
[41]122      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN
[3]123         IF(lwp) WRITE(numout,*)
[79]124         IF(lwp) WRITE(numout,*) '          free surface formulation'
125         IF( lk_isl ) THEN
[3]126            IF(lwp) WRITE(numout,cform_err)
127            IF(lwp) WRITE(numout,*) ' key_islands inconsistent with key_dynspg_fsc'
128            nstop = nstop + 1
129         ENDIF
[16]130      ELSEIF( lk_dynspg_rl ) THEN
[3]131         IF(lwp) WRITE(numout,*)
[79]132         IF(lwp) WRITE(numout,*) '          Rigid lid formulation'
[16]133      ELSE
[3]134         IF(lwp) WRITE(numout,cform_err)
[16]135         IF(lwp) WRITE(numout,*) '          Chose at least one surface pressure gradient calculation: free surface or rigid-lid'
136         nstop = nstop + 1
137      ENDIF
[79]138      IF( ( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) .AND. lk_dynspg_rl ) THEN
[16]139         IF(lwp) WRITE(numout,cform_err)
[3]140         IF(lwp) WRITE(numout,*) '          Chose between free surface or rigid-lid, not both'
141         nstop = nstop + 1
[16]142      ENDIF
[3]143
144      SELECT CASE ( nsolv )
145
146      CASE ( 1 )                ! preconditioned conjugate gradient solver
[79]147         IF(lwp) WRITE(numout,*) '          a preconditioned conjugate gradient solver is used'
[3]148
149      CASE ( 2 )                ! successive-over-relaxation solver
[79]150         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver is used'
[3]151
152      CASE ( 3 )                ! FETI solver
[79]153         IF(lwp) WRITE(numout,*) '          the FETI solver is used'
[16]154         IF( .NOT.lk_mpp ) THEN
[79]155            IF(lwp) WRITE(numout,cform_err)
[16]156            IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option'
[3]157            nstop = nstop + 1
[16]158         ELSE
159            IF( jpnij == 1 ) THEN
[79]160               IF(lwp) WRITE(numout,cform_err)
[16]161               IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor'
162               nstop = nstop + 1
163            ENDIF
[3]164         ENDIF
165         
166      CASE DEFAULT
167         IF(lwp) WRITE(numout,cform_err)
168         IF(lwp) WRITE(numout,*) '          bad flag value for nsolv = ', nsolv
169         nstop = nstop + 1
170         
171      END SELECT
172
[16]173      ! Grid-point at which the solver is applied
174      ! -----------------------------------------
[3]175
[16]176      IF( lk_dynspg_rl ) THEN       ! rigid-lid
177         IF( lk_mpp ) THEN
178            c_solver_pt = 'G'   ! G= F with special staff ??? which one?
179         ELSE
180            c_solver_pt = 'F'
181         ENDIF
182      ELSE                          ! free surface T-point
183         IF( lk_mpp ) THEN
184            c_solver_pt = 'S'   ! S=T with special staff ??? which one?
185         ELSE
186            c_solver_pt = 'T'
187         ENDIF
188      ENDIF
189
190
[3]191      ! Construction of the elliptic system matrix
192      ! ------------------------------------------
193
194      CALL sol_mat
195
196
[79]197      IF( lk_isl ) THEN
[3]198     
199         ! Islands in the domain
200         ! ---------------------
201
202         IF ( jpisl == 0 ) THEN
203             IF(lwp)WRITE(numout,cform_err)
204             IF(lwp)WRITE(numout,*) ' bad islands parameter jpisl =', jpisl
205             nstop = nstop + 1
206         ENDIF
207
208         ! open Island streamfunction statistic file
209         CALL ctlopn( numisp, 'islands.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
210            &         1     , numout        , lwp      , 1                         )
211   
212         CALL isl_dom       ! Island identification
213
214         CALL isl_bsf       ! Island barotropic stream function
215
216         CALL isl_mat       ! Comput and invert the island matrix
217
218         ! mbathy set to the number of w-level (minimum value 2)
219         DO jj = 1, jpj
220            DO ji = 1, jpi
221               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
222            END DO
223         END DO
224
225      ENDIF
226
227   END SUBROUTINE solver_init
228
229   !!======================================================================
230END MODULE solver
Note: See TracBrowser for help on using the repository browser.