[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 |
---|
[38] | 23 | USE dynspg_fsc_atsk |
---|
[3] | 24 | |
---|
| 25 | IMPLICIT NONE |
---|
| 26 | |
---|
| 27 | !!---------------------------------------------------------------------- |
---|
| 28 | !! OPA 9.0 , LODYC-IPSL (2003) |
---|
| 29 | !!---------------------------------------------------------------------- |
---|
| 30 | |
---|
| 31 | CONTAINS |
---|
| 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 l_isl=T) |
---|
[16] | 40 | !! * lk_dynspg_fsc = T : transport divergence system. No specific |
---|
[3] | 41 | !! treatment of islands. |
---|
| 42 | !! |
---|
| 43 | !! ** Method : |
---|
| 44 | !! - Compute the local depth of the water column at u- and v-point |
---|
[16] | 45 | !! (lk_dynspg_fsc = T) or its inverse (lk_dynspg_rl = T). |
---|
[3] | 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 l_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 |
---|
[16] | 62 | !! u- and v-point. (lk_dynspg_rl = T) |
---|
[3] | 63 | !! - hu, hv : masked local depth at u- and v- points |
---|
[16] | 64 | !! (lk_dynspg_fsc = T) |
---|
| 65 | !! - c_solver_pt : nature of the gridpoint at which the |
---|
| 66 | !! solver is applied |
---|
[3] | 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,*) 'ini_sol : 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,*) |
---|
| 113 | WRITE(numout,*) ' type of elliptic solver nsolv = ', nsolv |
---|
| 114 | WRITE(numout,*) ' maximum iterations for solver nmax = ', nmax |
---|
| 115 | WRITE(numout,*) ' absolute precision of solver eps = ', eps |
---|
| 116 | WRITE(numout,*) ' optimal coefficient of sor sor = ', sor |
---|
| 117 | IF(l_isl) WRITE(numout,*) ' absolute precision stream fct epsisl = ', epsisl |
---|
| 118 | IF(l_isl) WRITE(numout,*) ' maximum pcg iterations island nmisl = ', nmisl |
---|
| 119 | WRITE(numout,*) ' free surface parameter rnu = ', rnu |
---|
| 120 | WRITE(numout,*) |
---|
| 121 | ENDIF |
---|
| 122 | |
---|
[41] | 123 | IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN |
---|
[3] | 124 | IF(lwp) WRITE(numout,*) |
---|
| 125 | IF(lwp) WRITE(numout,*) ' *** free surface formulation' |
---|
| 126 | IF( l_isl ) THEN |
---|
| 127 | IF(lwp) WRITE(numout,cform_err) |
---|
| 128 | IF(lwp) WRITE(numout,*) ' key_islands inconsistent with key_dynspg_fsc' |
---|
| 129 | nstop = nstop + 1 |
---|
| 130 | ENDIF |
---|
[16] | 131 | ELSEIF( lk_dynspg_rl ) THEN |
---|
[3] | 132 | IF(lwp) WRITE(numout,*) |
---|
| 133 | IF(lwp) WRITE(numout,*) ' *** Rigid lid formulation' |
---|
[16] | 134 | ELSE |
---|
[3] | 135 | IF(lwp) WRITE(numout,cform_err) |
---|
[16] | 136 | IF(lwp) WRITE(numout,*) ' Chose at least one surface pressure gradient calculation: free surface or rigid-lid' |
---|
| 137 | nstop = nstop + 1 |
---|
| 138 | ENDIF |
---|
| 139 | IF( lk_dynspg_fsc .AND. lk_dynspg_rl ) THEN |
---|
| 140 | IF(lwp) WRITE(numout,cform_err) |
---|
[3] | 141 | IF(lwp) WRITE(numout,*) ' Chose between free surface or rigid-lid, not both' |
---|
| 142 | nstop = nstop + 1 |
---|
[16] | 143 | ENDIF |
---|
[3] | 144 | |
---|
| 145 | SELECT CASE ( nsolv ) |
---|
| 146 | |
---|
| 147 | CASE ( 1 ) ! preconditioned conjugate gradient solver |
---|
| 148 | IF(lwp) WRITE(numout,*) ' use a preconditioned conjugate gradient solver' |
---|
| 149 | |
---|
| 150 | CASE ( 2 ) ! successive-over-relaxation solver |
---|
| 151 | IF(lwp) WRITE(numout,*) ' use a successive-over-relaxation solver' |
---|
| 152 | |
---|
| 153 | CASE ( 3 ) ! FETI solver |
---|
| 154 | IF(lwp) WRITE(numout,*) ' use the FETI solver' |
---|
[16] | 155 | IF( .NOT.lk_mpp ) THEN |
---|
| 156 | IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option' |
---|
[3] | 157 | nstop = nstop + 1 |
---|
[16] | 158 | ELSE |
---|
| 159 | IF( jpnij == 1 ) THEN |
---|
| 160 | IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor' |
---|
| 161 | nstop = nstop + 1 |
---|
| 162 | ENDIF |
---|
[3] | 163 | ENDIF |
---|
| 164 | |
---|
| 165 | CASE DEFAULT |
---|
| 166 | IF(lwp) WRITE(numout,cform_err) |
---|
| 167 | IF(lwp) WRITE(numout,*) ' bad flag value for nsolv = ', nsolv |
---|
| 168 | nstop = nstop + 1 |
---|
| 169 | |
---|
| 170 | END SELECT |
---|
| 171 | |
---|
[16] | 172 | ! Grid-point at which the solver is applied |
---|
| 173 | ! ----------------------------------------- |
---|
[3] | 174 | |
---|
[16] | 175 | IF( lk_dynspg_rl ) THEN ! rigid-lid |
---|
| 176 | IF( lk_mpp ) THEN |
---|
| 177 | c_solver_pt = 'G' ! G= F with special staff ??? which one? |
---|
| 178 | ELSE |
---|
| 179 | c_solver_pt = 'F' |
---|
| 180 | ENDIF |
---|
| 181 | ELSE ! free surface T-point |
---|
| 182 | IF( lk_mpp ) THEN |
---|
| 183 | c_solver_pt = 'S' ! S=T with special staff ??? which one? |
---|
| 184 | ELSE |
---|
| 185 | c_solver_pt = 'T' |
---|
| 186 | ENDIF |
---|
| 187 | ENDIF |
---|
| 188 | |
---|
| 189 | |
---|
[3] | 190 | ! Construction of the elliptic system matrix |
---|
| 191 | ! ------------------------------------------ |
---|
| 192 | |
---|
| 193 | CALL sol_mat |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | IF( l_isl ) THEN |
---|
| 197 | |
---|
| 198 | ! Islands in the domain |
---|
| 199 | ! --------------------- |
---|
| 200 | |
---|
| 201 | IF ( jpisl == 0 ) THEN |
---|
| 202 | IF(lwp)WRITE(numout,cform_err) |
---|
| 203 | IF(lwp)WRITE(numout,*) ' bad islands parameter jpisl =', jpisl |
---|
| 204 | nstop = nstop + 1 |
---|
| 205 | ENDIF |
---|
| 206 | |
---|
| 207 | ! open Island streamfunction statistic file |
---|
| 208 | CALL ctlopn( numisp, 'islands.stat', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', & |
---|
| 209 | & 1 , numout , lwp , 1 ) |
---|
| 210 | |
---|
| 211 | CALL isl_dom ! Island identification |
---|
| 212 | |
---|
| 213 | CALL isl_bsf ! Island barotropic stream function |
---|
| 214 | |
---|
| 215 | CALL isl_mat ! Comput and invert the island matrix |
---|
| 216 | |
---|
| 217 | ! mbathy set to the number of w-level (minimum value 2) |
---|
| 218 | DO jj = 1, jpj |
---|
| 219 | DO ji = 1, jpi |
---|
| 220 | mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1 |
---|
| 221 | END DO |
---|
| 222 | END DO |
---|
| 223 | |
---|
| 224 | ENDIF |
---|
| 225 | |
---|
| 226 | END SUBROUTINE solver_init |
---|
| 227 | |
---|
| 228 | !!====================================================================== |
---|
| 229 | END MODULE solver |
---|