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

Last change on this file since 389 was 389, checked in by opalod, 18 years ago

RB:nemo_v1_update_038: first integration of Agrif :

  • configuration parameters are just integer when agrif is used
  • add call to agrif routines with key_agrif
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.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_oce      ! choice/control of key cpp for surface pressure gradient
22
23   IMPLICIT NONE
24
25   !!----------------------------------------------------------------------
26   !!   OPA 9.0 , LOCEAN-IPSL (2005)
27   !! $Header$
28   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
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_flt = 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_flt = 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_flt = 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      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
78      !!----------------------------------------------------------------------
79      !! * Local declarations
80      INTEGER :: ji, jj   ! dummy loop indices
81      CHARACTER(len=80) :: clname
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      clname = 'solver.stat'
92      CALL ctlopn( numsol, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
93                   1, numout, lwp, 1 )
94
95
96      ! 0. Define the solver parameters
97      !    ----------------------------
98      ! Namelist namsol : elliptic solver / islands / free surface
99      REWIND( numnam )
100      READ  ( numnam, namsol )
101
102#if defined key_feti
103      ! FETI algorithm, we force nsolv at 3
104      nsolv = 3
105#endif
106
107
108      ! 0. Parameter control and print
109      !    ---------------------------
110
111      ! Control print
112      IF(lwp) WRITE(numout,*) '          Namelist namsol : set solver parameters'
113
114      IF(lwp) THEN
115         WRITE(numout,*) '             type of elliptic solver            nsolv    = ', nsolv
116         WRITE(numout,*) '             absolute/relative (0/1) precision  nsol_arp = ', nsol_arp
117         WRITE(numout,*) '             minimum iterations for solver      nmin     = ', nmin
118         WRITE(numout,*) '             maximum iterations for solver      nmax     = ', nmax
119         WRITE(numout,*) '             frequency for test                 nmod     = ', nmod
120         WRITE(numout,*) '             absolute precision of solver       eps      = ', eps
121         WRITE(numout,*) '             absolute precision for SOR solver  resmax   = ', resmax
122         WRITE(numout,*) '             optimal coefficient of sor         sor      = ', sor
123         IF(lk_isl) WRITE(numout,*) '             absolute precision stream fct    epsisl   = ', epsisl
124         IF(lk_isl) WRITE(numout,*) '             maximum pcg iterations island    nmisl    = ', nmisl
125         WRITE(numout,*) '             free surface parameter         rnu    = ', rnu
126         WRITE(numout,*)
127      ENDIF
128
129      IF( lk_dynspg_flt ) THEN
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) '          free surface formulation'
132         IF( lk_isl ) THEN
133            IF(lwp) WRITE(numout,cform_err)
134            IF(lwp) WRITE(numout,*) ' key_islands inconsistent with key_dynspg_flt'
135            nstop = nstop + 1
136         ENDIF
137      ELSEIF( lk_dynspg_rl ) THEN
138         IF(lwp) WRITE(numout,*)
139         IF(lwp) WRITE(numout,*) '          Rigid lid formulation'
140      ELSE
141         IF(lwp) WRITE(numout,cform_err)
142         IF(lwp) WRITE(numout,*) '          Choose only one surface pressure gradient calculation: filtered or rigid-lid'
143         IF(lwp) WRITE(numout,*) '          Should not call this routine if dynspg_exp or dynspg_ts has been chosen'
144         nstop = nstop + 1
145      ENDIF
146      IF( lk_dynspg_flt .AND. lk_dynspg_rl ) THEN
147         IF(lwp) WRITE(numout,cform_err)
148         IF(lwp) WRITE(numout,*) '          Chose between free surface or rigid-lid, not both'
149         nstop = nstop + 1
150      ENDIF
151
152      SELECT CASE ( nsolv )
153
154      CASE ( 1 )                ! preconditioned conjugate gradient solver
155         IF(lwp) WRITE(numout,*) '          a preconditioned conjugate gradient solver is used'
156         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
157            IF(lwp) WRITE(numout,cform_err)
158            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero'
159            nstop = nstop + 1
160         ENDIF
161
162      CASE ( 2 )                ! successive-over-relaxation solver
163         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver is used'
164         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
165            IF(lwp) WRITE(numout,cform_err)
166            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero'
167            nstop = nstop + 1
168         ENDIF
169
170      CASE ( 3 )                ! FETI solver
171         IF(lwp) WRITE(numout,*) '          the FETI solver is used'
172         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
173            IF(lwp) WRITE(numout,cform_err)
174            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero'
175            nstop = nstop + 1
176         ENDIF
177         IF( .NOT.lk_mpp ) THEN
178            IF(lwp) WRITE(numout,cform_err)
179            IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option'
180            nstop = nstop + 1
181         ELSE
182            IF( jpnij == 1 ) THEN
183               IF(lwp) WRITE(numout,cform_err)
184               IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor'
185               nstop = nstop + 1
186            ENDIF
187         ENDIF
188         
189      CASE ( 4 )                ! successive-over-relaxation solver with extra outer halo
190         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver with extra outer halo is used'
191         IF(lwp) WRITE(numout,*) '          with jpr2di =', jpr2di, ' and  jpr2dj =', jpr2dj
192         IF( .NOT. lk_mpp .AND. jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
193            IF(lwp) WRITE(numout,cform_err)
194            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj are not equal to zero'
195            IF(lwp) WRITE(numout,*) ' In this case this algorithm should be used only with the key_mpp_... option'
196            nstop = nstop + 1
197         ELSE
198            IF( ( ( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) .OR. ( jpni /= 1 ) ) &
199              &  .AND. ( jpr2di /= jpr2dj ) ) THEN 
200               IF(lwp) WRITE(numout,cform_err)
201               IF(lwp) WRITE(numout,*) '          jpr2di should be equal to jpr2dj'
202               nstop = nstop + 1
203            ENDIF
204         ENDIF
205
206      CASE DEFAULT
207         IF(lwp) WRITE(numout,cform_err)
208         IF(lwp) WRITE(numout,*) '          bad flag value for nsolv = ', nsolv
209         nstop = nstop + 1
210         
211      END SELECT
212
213      ! Grid-point at which the solver is applied
214      ! -----------------------------------------
215
216      IF( lk_dynspg_rl ) THEN       ! rigid-lid
217         IF( lk_mpp ) THEN
218            c_solver_pt = 'G'   ! G= F with special staff ??? which one?
219         ELSE
220            c_solver_pt = 'F'
221         ENDIF
222      ELSE                          ! free surface T-point
223         IF( lk_mpp ) THEN
224            c_solver_pt = 'S'   ! S=T with special staff ??? which one?
225         ELSE
226            c_solver_pt = 'T'
227         ENDIF
228      ENDIF
229
230
231      ! Construction of the elliptic system matrix
232      ! ------------------------------------------
233
234      CALL sol_mat
235
236
237      IF( lk_isl ) THEN
238     
239         ! Islands in the domain
240         ! ---------------------
241
242         IF ( jpisl == 0 ) THEN
243             IF(lwp)WRITE(numout,cform_err)
244             IF(lwp)WRITE(numout,*) ' bad islands parameter jpisl =', jpisl
245             nstop = nstop + 1
246         ENDIF
247
248         ! open Island streamfunction statistic file
249         CALL ctlopn( numisp, 'islands.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   &
250            &         1     , numout        , lwp      , 1                         )
251   
252         CALL isl_dom       ! Island identification
253
254         CALL isl_bsf       ! Island barotropic stream function
255
256         CALL isl_mat       ! Comput and invert the island matrix
257
258         ! mbathy set to the number of w-level (minimum value 2)
259         DO jj = 1, jpj
260            DO ji = 1, jpi
261               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
262            END DO
263         END DO
264
265      ENDIF
266
267   END SUBROUTINE solver_init
268
269   !!======================================================================
270END MODULE solver
Note: See TracBrowser for help on using the repository browser.