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 tags/nemo_dev_x6/NEMO/OPA_SRC/SOL – NEMO

source: tags/nemo_dev_x6/NEMO/OPA_SRC/SOL/solver.F90 @ 3532

Last change on this file since 3532 was 120, checked in by opalod, 20 years ago

CT : UPDATE070 : Optimisation of the Red-Black SOR algorithm convergence
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 KB
Line 
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)
20   USE lib_mpp
21   USE dynspg_rl       
22   USE dynspg_fsc     
23   USE dynspg_fsc_atsk     
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
39      !!         and islands initialization (if lk_isl=T)
40      !!       * lk_dynspg_fsc = T : transport divergence system. No specific
41      !!         treatment of islands.
42      !!     
43      !! ** Method :
44      !!       - Compute the local depth of the water column at u- and v-point
45      !!      (lk_dynspg_fsc = T) or its inverse (lk_dynspg_rl = T).
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.
53      !!       - island (if lk_isl=T)
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
62      !!                                u- and v-point. (lk_dynspg_rl = T)
63      !!             - hu, hv   : masked local depth at u- and v- points
64      !!                                (lk_dynspg_fsc = T)
65      !!             - c_solver_pt : nature of the gridpoint at which the
66      !!                                solver is applied
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, nsol_arp, nmin, nmax, nmod, eps, resmax, sor, epsisl, nmisl, rnu
82      !!----------------------------------------------------------------------
83
84      IF(lwp) WRITE(numout,*)
85      IF(lwp) WRITE(numout,*) 'solver_init : solver to compute the surface pressure gradient'
86      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
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,*) '             absolute/relative (0/1) precision  nsol_arp = ', nsol_arp
114         WRITE(numout,*) '             minimum iterations for solver      nmin     = ', nmin
115         WRITE(numout,*) '             maximum iterations for solver      nmax     = ', nmax
116         WRITE(numout,*) '             frequency for test                 nmod     = ', nmod
117         WRITE(numout,*) '             absolute precision of solver       eps      = ', eps
118         WRITE(numout,*) '             absolute precision for SOR solver  resmax   = ', resmax
119         WRITE(numout,*) '             optimal coefficient of sor         sor      = ', sor
120         IF(lk_isl) WRITE(numout,*) '             absolute precision stream fct    epsisl   = ', epsisl
121         IF(lk_isl) WRITE(numout,*) '             maximum pcg iterations island    nmisl    = ', nmisl
122         WRITE(numout,*) '             free surface parameter         rnu    = ', rnu
123         WRITE(numout,*)
124      ENDIF
125
126      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN
127         IF(lwp) WRITE(numout,*)
128         IF(lwp) WRITE(numout,*) '          free surface formulation'
129         IF( lk_isl ) THEN
130            IF(lwp) WRITE(numout,cform_err)
131            IF(lwp) WRITE(numout,*) ' key_islands inconsistent with key_dynspg_fsc'
132            nstop = nstop + 1
133         ENDIF
134      ELSEIF( lk_dynspg_rl ) THEN
135         IF(lwp) WRITE(numout,*)
136         IF(lwp) WRITE(numout,*) '          Rigid lid formulation'
137      ELSE
138         IF(lwp) WRITE(numout,cform_err)
139         IF(lwp) WRITE(numout,*) '          Chose at least one surface pressure gradient calculation: free surface or rigid-lid'
140         nstop = nstop + 1
141      ENDIF
142      IF( ( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) .AND. lk_dynspg_rl ) THEN
143         IF(lwp) WRITE(numout,cform_err)
144         IF(lwp) WRITE(numout,*) '          Chose between free surface or rigid-lid, not both'
145         nstop = nstop + 1
146      ENDIF
147
148      SELECT CASE ( nsolv )
149
150      CASE ( 1 )                ! preconditioned conjugate gradient solver
151         IF(lwp) WRITE(numout,*) '          a preconditioned conjugate gradient solver is used'
152
153      CASE ( 2 )                ! successive-over-relaxation solver
154         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver is used'
155
156      CASE ( 3 )                ! FETI solver
157         IF(lwp) WRITE(numout,*) '          the FETI solver is used'
158         IF( .NOT.lk_mpp ) THEN
159            IF(lwp) WRITE(numout,cform_err)
160            IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option'
161            nstop = nstop + 1
162         ELSE
163            IF( jpnij == 1 ) THEN
164               IF(lwp) WRITE(numout,cform_err)
165               IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor'
166               nstop = nstop + 1
167            ENDIF
168         ENDIF
169         
170      CASE DEFAULT
171         IF(lwp) WRITE(numout,cform_err)
172         IF(lwp) WRITE(numout,*) '          bad flag value for nsolv = ', nsolv
173         nstop = nstop + 1
174         
175      END SELECT
176
177      ! Grid-point at which the solver is applied
178      ! -----------------------------------------
179
180      IF( lk_dynspg_rl ) THEN       ! rigid-lid
181         IF( lk_mpp ) THEN
182            c_solver_pt = 'G'   ! G= F with special staff ??? which one?
183         ELSE
184            c_solver_pt = 'F'
185         ENDIF
186      ELSE                          ! free surface T-point
187         IF( lk_mpp ) THEN
188            c_solver_pt = 'S'   ! S=T with special staff ??? which one?
189         ELSE
190            c_solver_pt = 'T'
191         ENDIF
192      ENDIF
193
194
195      ! Construction of the elliptic system matrix
196      ! ------------------------------------------
197
198      CALL sol_mat
199
200
201      IF( lk_isl ) THEN
202     
203         ! Islands in the domain
204         ! ---------------------
205
206         IF ( jpisl == 0 ) THEN
207             IF(lwp)WRITE(numout,cform_err)
208             IF(lwp)WRITE(numout,*) ' bad islands parameter jpisl =', jpisl
209             nstop = nstop + 1
210         ENDIF
211
212         ! open Island streamfunction statistic file
213         CALL ctlopn( numisp, 'islands.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
214            &         1     , numout        , lwp      , 1                         )
215   
216         CALL isl_dom       ! Island identification
217
218         CALL isl_bsf       ! Island barotropic stream function
219
220         CALL isl_mat       ! Comput and invert the island matrix
221
222         ! mbathy set to the number of w-level (minimum value 2)
223         DO jj = 1, jpj
224            DO ji = 1, jpi
225               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
226            END DO
227         END DO
228
229      ENDIF
230
231   END SUBROUTINE solver_init
232
233   !!======================================================================
234END MODULE solver
Note: See TracBrowser for help on using the repository browser.