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 @ 111

Last change on this file since 111 was 111, 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.5 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, 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 (1/0) 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,*) '             optimal coefficient of sor         sor      = ', sor
119         IF(lk_isl) WRITE(numout,*) '             absolute precision stream fct    epsisl   = ', epsisl
120         IF(lk_isl) WRITE(numout,*) '             maximum pcg iterations island    nmisl    = ', nmisl
121         WRITE(numout,*) '             free surface parameter         rnu    = ', rnu
122         WRITE(numout,*)
123      ENDIF
124
125      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN
126         IF(lwp) WRITE(numout,*)
127         IF(lwp) WRITE(numout,*) '          free surface formulation'
128         IF( lk_isl ) THEN
129            IF(lwp) WRITE(numout,cform_err)
130            IF(lwp) WRITE(numout,*) ' key_islands inconsistent with key_dynspg_fsc'
131            nstop = nstop + 1
132         ENDIF
133      ELSEIF( lk_dynspg_rl ) THEN
134         IF(lwp) WRITE(numout,*)
135         IF(lwp) WRITE(numout,*) '          Rigid lid formulation'
136      ELSE
137         IF(lwp) WRITE(numout,cform_err)
138         IF(lwp) WRITE(numout,*) '          Chose at least one surface pressure gradient calculation: free surface or rigid-lid'
139         nstop = nstop + 1
140      ENDIF
141      IF( ( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) .AND. lk_dynspg_rl ) THEN
142         IF(lwp) WRITE(numout,cform_err)
143         IF(lwp) WRITE(numout,*) '          Chose between free surface or rigid-lid, not both'
144         nstop = nstop + 1
145      ENDIF
146
147      SELECT CASE ( nsolv )
148
149      CASE ( 1 )                ! preconditioned conjugate gradient solver
150         IF(lwp) WRITE(numout,*) '          a preconditioned conjugate gradient solver is used'
151
152      CASE ( 2 )                ! successive-over-relaxation solver
153         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver is used'
154
155      CASE ( 3 )                ! FETI solver
156         IF(lwp) WRITE(numout,*) '          the FETI solver is used'
157         IF( .NOT.lk_mpp ) THEN
158            IF(lwp) WRITE(numout,cform_err)
159            IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option'
160            nstop = nstop + 1
161         ELSE
162            IF( jpnij == 1 ) THEN
163               IF(lwp) WRITE(numout,cform_err)
164               IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor'
165               nstop = nstop + 1
166            ENDIF
167         ENDIF
168         
169      CASE DEFAULT
170         IF(lwp) WRITE(numout,cform_err)
171         IF(lwp) WRITE(numout,*) '          bad flag value for nsolv = ', nsolv
172         nstop = nstop + 1
173         
174      END SELECT
175
176      ! Grid-point at which the solver is applied
177      ! -----------------------------------------
178
179      IF( lk_dynspg_rl ) THEN       ! rigid-lid
180         IF( lk_mpp ) THEN
181            c_solver_pt = 'G'   ! G= F with special staff ??? which one?
182         ELSE
183            c_solver_pt = 'F'
184         ENDIF
185      ELSE                          ! free surface T-point
186         IF( lk_mpp ) THEN
187            c_solver_pt = 'S'   ! S=T with special staff ??? which one?
188         ELSE
189            c_solver_pt = 'T'
190         ENDIF
191      ENDIF
192
193
194      ! Construction of the elliptic system matrix
195      ! ------------------------------------------
196
197      CALL sol_mat
198
199
200      IF( lk_isl ) THEN
201     
202         ! Islands in the domain
203         ! ---------------------
204
205         IF ( jpisl == 0 ) THEN
206             IF(lwp)WRITE(numout,cform_err)
207             IF(lwp)WRITE(numout,*) ' bad islands parameter jpisl =', jpisl
208             nstop = nstop + 1
209         ENDIF
210
211         ! open Island streamfunction statistic file
212         CALL ctlopn( numisp, 'islands.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
213            &         1     , numout        , lwp      , 1                         )
214   
215         CALL isl_dom       ! Island identification
216
217         CALL isl_bsf       ! Island barotropic stream function
218
219         CALL isl_mat       ! Comput and invert the island matrix
220
221         ! mbathy set to the number of w-level (minimum value 2)
222         DO jj = 1, jpj
223            DO ji = 1, jpi
224               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
225            END DO
226         END DO
227
228      ENDIF
229
230   END SUBROUTINE solver_init
231
232   !!======================================================================
233END MODULE solver
Note: See TracBrowser for help on using the repository browser.