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
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, nmax, 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,*) '             maximum iterations for solver  nmax   = ', nmax
114         WRITE(numout,*) '             absolute precision of solver   eps    = ', eps
115         WRITE(numout,*) '             optimal coefficient of sor     sor    = ', sor
116         IF(lk_isl) WRITE(numout,*) '             absolute precision stream fct  epsisl = ', epsisl
117         IF(lk_isl) WRITE(numout,*) '             maximum pcg iterations island  nmisl  = ', nmisl
118         WRITE(numout,*) '             free surface parameter         rnu    = ', rnu
119         WRITE(numout,*)
120      ENDIF
121
122      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN
123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) '          free surface formulation'
125         IF( lk_isl ) THEN
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
130      ELSEIF( lk_dynspg_rl ) THEN
131         IF(lwp) WRITE(numout,*)
132         IF(lwp) WRITE(numout,*) '          Rigid lid formulation'
133      ELSE
134         IF(lwp) WRITE(numout,cform_err)
135         IF(lwp) WRITE(numout,*) '          Chose at least one surface pressure gradient calculation: free surface or rigid-lid'
136         nstop = nstop + 1
137      ENDIF
138      IF( ( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) .AND. lk_dynspg_rl ) THEN
139         IF(lwp) WRITE(numout,cform_err)
140         IF(lwp) WRITE(numout,*) '          Chose between free surface or rigid-lid, not both'
141         nstop = nstop + 1
142      ENDIF
143
144      SELECT CASE ( nsolv )
145
146      CASE ( 1 )                ! preconditioned conjugate gradient solver
147         IF(lwp) WRITE(numout,*) '          a preconditioned conjugate gradient solver is used'
148
149      CASE ( 2 )                ! successive-over-relaxation solver
150         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver is used'
151
152      CASE ( 3 )                ! FETI solver
153         IF(lwp) WRITE(numout,*) '          the FETI solver is used'
154         IF( .NOT.lk_mpp ) THEN
155            IF(lwp) WRITE(numout,cform_err)
156            IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option'
157            nstop = nstop + 1
158         ELSE
159            IF( jpnij == 1 ) THEN
160               IF(lwp) WRITE(numout,cform_err)
161               IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor'
162               nstop = nstop + 1
163            ENDIF
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
173      ! Grid-point at which the solver is applied
174      ! -----------------------------------------
175
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
191      ! Construction of the elliptic system matrix
192      ! ------------------------------------------
193
194      CALL sol_mat
195
196
197      IF( lk_isl ) THEN
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.