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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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