[3] | 1 | MODULE 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) |
---|
[16] | 20 | USE lib_mpp |
---|
| 21 | USE dynspg_rl |
---|
| 22 | USE dynspg_fsc |
---|
[79] | 23 | USE dynspg_fsc_atsk |
---|
[3] | 24 | |
---|
| 25 | IMPLICIT NONE |
---|
| 26 | |
---|
| 27 | !!---------------------------------------------------------------------- |
---|
[247] | 28 | !! OPA 9.0 , LOCEAN-IPSL (2005) |
---|
| 29 | !! $Header$ |
---|
| 30 | !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt |
---|
[3] | 31 | !!---------------------------------------------------------------------- |
---|
| 32 | |
---|
| 33 | CONTAINS |
---|
| 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 |
---|
[79] | 41 | !! and islands initialization (if lk_isl=T) |
---|
[16] | 42 | !! * lk_dynspg_fsc = T : transport divergence system. No specific |
---|
[3] | 43 | !! treatment of islands. |
---|
| 44 | !! |
---|
| 45 | !! ** Method : |
---|
| 46 | !! - Compute the local depth of the water column at u- and v-point |
---|
[16] | 47 | !! (lk_dynspg_fsc = T) or its inverse (lk_dynspg_rl = T). |
---|
[3] | 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. |
---|
[79] | 55 | !! - island (if lk_isl=T) |
---|
[3] | 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 |
---|
[16] | 64 | !! u- and v-point. (lk_dynspg_rl = T) |
---|
[3] | 65 | !! - hu, hv : masked local depth at u- and v- points |
---|
[16] | 66 | !! (lk_dynspg_fsc = T) |
---|
| 67 | !! - c_solver_pt : nature of the gridpoint at which the |
---|
| 68 | !! solver is applied |
---|
[3] | 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 | |
---|
[120] | 83 | NAMELIST/namsol/ nsolv, nsol_arp, nmin, nmax, nmod, eps, resmax, sor, epsisl, nmisl, rnu |
---|
[3] | 84 | !!---------------------------------------------------------------------- |
---|
| 85 | |
---|
| 86 | IF(lwp) WRITE(numout,*) |
---|
[79] | 87 | IF(lwp) WRITE(numout,*) 'solver_init : solver to compute the surface pressure gradient' |
---|
| 88 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' |
---|
[3] | 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 |
---|
[111] | 114 | WRITE(numout,*) ' type of elliptic solver nsolv = ', nsolv |
---|
[120] | 115 | WRITE(numout,*) ' absolute/relative (0/1) precision nsol_arp = ', nsol_arp |
---|
[111] | 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] | 120 | WRITE(numout,*) ' absolute precision for SOR solver resmax = ', resmax |
---|
[111] | 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 |
---|
[3] | 124 | WRITE(numout,*) ' free surface parameter rnu = ', rnu |
---|
| 125 | WRITE(numout,*) |
---|
| 126 | ENDIF |
---|
| 127 | |
---|
[41] | 128 | IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN |
---|
[3] | 129 | IF(lwp) WRITE(numout,*) |
---|
[79] | 130 | IF(lwp) WRITE(numout,*) ' free surface formulation' |
---|
| 131 | IF( lk_isl ) THEN |
---|
[3] | 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 |
---|
[16] | 136 | ELSEIF( lk_dynspg_rl ) THEN |
---|
[3] | 137 | IF(lwp) WRITE(numout,*) |
---|
[79] | 138 | IF(lwp) WRITE(numout,*) ' Rigid lid formulation' |
---|
[16] | 139 | ELSE |
---|
[3] | 140 | IF(lwp) WRITE(numout,cform_err) |
---|
[16] | 141 | IF(lwp) WRITE(numout,*) ' Chose at least one surface pressure gradient calculation: free surface or rigid-lid' |
---|
| 142 | nstop = nstop + 1 |
---|
| 143 | ENDIF |
---|
[79] | 144 | IF( ( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) .AND. lk_dynspg_rl ) THEN |
---|
[16] | 145 | IF(lwp) WRITE(numout,cform_err) |
---|
[3] | 146 | IF(lwp) WRITE(numout,*) ' Chose between free surface or rigid-lid, not both' |
---|
| 147 | nstop = nstop + 1 |
---|
[16] | 148 | ENDIF |
---|
[3] | 149 | |
---|
| 150 | SELECT CASE ( nsolv ) |
---|
| 151 | |
---|
| 152 | CASE ( 1 ) ! preconditioned conjugate gradient solver |
---|
[79] | 153 | IF(lwp) WRITE(numout,*) ' a preconditioned conjugate gradient solver is used' |
---|
[312] | 154 | IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN |
---|
| 155 | IF(lwp) WRITE(numout,cform_err) |
---|
| 156 | IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero' |
---|
| 157 | nstop = nstop + 1 |
---|
| 158 | ENDIF |
---|
[3] | 159 | |
---|
| 160 | CASE ( 2 ) ! successive-over-relaxation solver |
---|
[79] | 161 | IF(lwp) WRITE(numout,*) ' a successive-over-relaxation solver is used' |
---|
[312] | 162 | IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN |
---|
| 163 | IF(lwp) WRITE(numout,cform_err) |
---|
| 164 | IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero' |
---|
| 165 | nstop = nstop + 1 |
---|
| 166 | ENDIF |
---|
[3] | 167 | |
---|
| 168 | CASE ( 3 ) ! FETI solver |
---|
[79] | 169 | IF(lwp) WRITE(numout,*) ' the FETI solver is used' |
---|
[312] | 170 | IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN |
---|
| 171 | IF(lwp) WRITE(numout,cform_err) |
---|
| 172 | IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero' |
---|
| 173 | nstop = nstop + 1 |
---|
| 174 | ENDIF |
---|
[16] | 175 | IF( .NOT.lk_mpp ) THEN |
---|
[79] | 176 | IF(lwp) WRITE(numout,cform_err) |
---|
[16] | 177 | IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option' |
---|
[3] | 178 | nstop = nstop + 1 |
---|
[16] | 179 | ELSE |
---|
| 180 | IF( jpnij == 1 ) THEN |
---|
[79] | 181 | IF(lwp) WRITE(numout,cform_err) |
---|
[16] | 182 | IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor' |
---|
| 183 | nstop = nstop + 1 |
---|
| 184 | ENDIF |
---|
[3] | 185 | ENDIF |
---|
| 186 | |
---|
[312] | 187 | CASE ( 4 ) ! successive-over-relaxation solver with extra outer halo |
---|
| 188 | IF(lwp) WRITE(numout,*) ' a successive-over-relaxation solver with extra outer halo is used' |
---|
| 189 | IF(lwp) WRITE(numout,*) ' with jpr2di =', jpr2di, ' and jpr2dj =', jpr2dj |
---|
| 190 | IF( .NOT. lk_mpp .AND. jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN |
---|
| 191 | IF(lwp) WRITE(numout,cform_err) |
---|
| 192 | IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj are not equal to zero' |
---|
| 193 | IF(lwp) WRITE(numout,*) ' In this case this algorithm should be used only with the key_mpp_... option' |
---|
| 194 | nstop = nstop + 1 |
---|
| 195 | ELSE |
---|
| 196 | IF( ( ( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) .OR. ( jpni /= 1 ) ) & |
---|
| 197 | & .AND. ( jpr2di /= jpr2dj ) ) THEN |
---|
| 198 | IF(lwp) WRITE(numout,cform_err) |
---|
| 199 | IF(lwp) WRITE(numout,*) ' jpr2di should be equal to jpr2dj' |
---|
| 200 | nstop = nstop + 1 |
---|
| 201 | ENDIF |
---|
| 202 | ENDIF |
---|
| 203 | |
---|
[3] | 204 | CASE DEFAULT |
---|
| 205 | IF(lwp) WRITE(numout,cform_err) |
---|
| 206 | IF(lwp) WRITE(numout,*) ' bad flag value for nsolv = ', nsolv |
---|
| 207 | nstop = nstop + 1 |
---|
| 208 | |
---|
| 209 | END SELECT |
---|
| 210 | |
---|
[16] | 211 | ! Grid-point at which the solver is applied |
---|
| 212 | ! ----------------------------------------- |
---|
[3] | 213 | |
---|
[16] | 214 | IF( lk_dynspg_rl ) THEN ! rigid-lid |
---|
| 215 | IF( lk_mpp ) THEN |
---|
| 216 | c_solver_pt = 'G' ! G= F with special staff ??? which one? |
---|
| 217 | ELSE |
---|
| 218 | c_solver_pt = 'F' |
---|
| 219 | ENDIF |
---|
| 220 | ELSE ! free surface T-point |
---|
| 221 | IF( lk_mpp ) THEN |
---|
| 222 | c_solver_pt = 'S' ! S=T with special staff ??? which one? |
---|
| 223 | ELSE |
---|
| 224 | c_solver_pt = 'T' |
---|
| 225 | ENDIF |
---|
| 226 | ENDIF |
---|
| 227 | |
---|
| 228 | |
---|
[3] | 229 | ! Construction of the elliptic system matrix |
---|
| 230 | ! ------------------------------------------ |
---|
| 231 | |
---|
| 232 | CALL sol_mat |
---|
| 233 | |
---|
| 234 | |
---|
[79] | 235 | IF( lk_isl ) THEN |
---|
[3] | 236 | |
---|
| 237 | ! Islands in the domain |
---|
| 238 | ! --------------------- |
---|
| 239 | |
---|
| 240 | IF ( jpisl == 0 ) THEN |
---|
| 241 | IF(lwp)WRITE(numout,cform_err) |
---|
| 242 | IF(lwp)WRITE(numout,*) ' bad islands parameter jpisl =', jpisl |
---|
| 243 | nstop = nstop + 1 |
---|
| 244 | ENDIF |
---|
| 245 | |
---|
| 246 | ! open Island streamfunction statistic file |
---|
| 247 | CALL ctlopn( numisp, 'islands.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', & |
---|
| 248 | & 1 , numout , lwp , 1 ) |
---|
| 249 | |
---|
| 250 | CALL isl_dom ! Island identification |
---|
| 251 | |
---|
| 252 | CALL isl_bsf ! Island barotropic stream function |
---|
| 253 | |
---|
| 254 | CALL isl_mat ! Comput and invert the island matrix |
---|
| 255 | |
---|
| 256 | ! mbathy set to the number of w-level (minimum value 2) |
---|
| 257 | DO jj = 1, jpj |
---|
| 258 | DO ji = 1, jpi |
---|
| 259 | mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1 |
---|
| 260 | END DO |
---|
| 261 | END DO |
---|
| 262 | |
---|
| 263 | ENDIF |
---|
| 264 | |
---|
| 265 | END SUBROUTINE solver_init |
---|
| 266 | |
---|
| 267 | !!====================================================================== |
---|
| 268 | END MODULE solver |
---|