Changeset 16 for trunk/NEMO/OPA_SRC/SOL/solpcg.F90
- Timestamp:
- 2004-02-17T09:06:15+01:00 (20 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/SOL/solpcg.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SOL/solpcg.F90
r3 r16 14 14 USE lib_mpp ! distributed memory computing 15 15 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 16 USE in_out_manager ! I/O manager 16 17 17 18 IMPLICIT NONE … … 32 33 !! 33 34 !! ** Purpose : Solve the ellipic equation for the barotropic stream 34 !! function system ( default option) or the transport divergence35 !! system ( "key_dynspg_fsc") using a diagonal preconditionned35 !! function system (lk_dynspg_rl=T) or the transport divergence 36 !! system (lk_dynspg_fsc=T) using a diagonal preconditionned 36 37 !! conjugate gradient method. 37 38 !! In the former case, the barotropic stream function trend has a … … 93 94 ! !================ 94 95 95 !,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 96 97 IF( jn == 1 ) THEN 98 99 ! 1.0 Initialization of the algorithm 100 ! ----------------------------------- 101 102 #if defined key_dynspg_fsc 103 # if defined key_mpp 104 ! Mpp: export boundary values to neighbouring processors 105 CALL lbc_lnk( gcx, 'S', 1. ) 106 # else 107 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 108 CALL lbc_lnk( gcx, 'T', 1. ) 109 # endif 110 #else 111 # if defined key_mpp 112 ! ... Mpp: export boundary values to neighbouring processors 113 CALL lbc_lnk( gcx, 'G', 1. ) 114 # else 115 ! ... mono- or macro-tasking: F-point, >0, 2D array, no slab 116 CALL lbc_lnk( gcx, 'F', 1. ) 117 # endif 118 #endif 96 IF( jn == 1 ) THEN ! Initialization of the algorithm 119 97 120 !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,121 98 CALL lbc_lnk( gcx, c_solver_pt, 1. ) ! lateral boundary condition 99 122 100 ! gcr = gcb-a.gcx 123 101 ! gcdes = gsr 124 125 102 DO jj = 2, jpjm1 126 103 DO ji = fs_2, fs_jpim1 ! vector opt. 127 zgcad = bmask(ji,jj)*( & 128 gcb(ji,jj ) - gcx(ji ,jj ) & 129 - gcp(ji,jj,1)*gcx(ji ,jj-1) & 130 - gcp(ji,jj,2)*gcx(ji-1,jj ) & 131 - gcp(ji,jj,3)*gcx(ji+1,jj ) & 132 - gcp(ji,jj,4)*gcx(ji ,jj+1) ) 104 zgcad = bmask(ji,jj) * ( gcb(ji,jj ) - gcx(ji ,jj ) & 105 & - gcp(ji,jj,1) * gcx(ji ,jj-1) & 106 & - gcp(ji,jj,2) * gcx(ji-1,jj ) & 107 & - gcp(ji,jj,3) * gcx(ji+1,jj ) & 108 & - gcp(ji,jj,4) * gcx(ji ,jj+1) ) 133 109 gcr (ji,jj) = zgcad 134 110 gcdes(ji,jj) = zgcad … … 136 112 END DO 137 113 138 !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,139 140 114 rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) ) 141 142 #if defined key_mpp 143 ! Mpp: sum over all the global domain 144 CALL mpp_sum( rnorme ) 145 #endif 115 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 146 116 rr = rnorme 147 117 148 ENDIF 149 !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,, 118 ENDIF 150 119 120 ! ! Algorithm 151 121 152 ! 1.1 Algorithm 153 ! ------------- 122 CALL lbc_lnk( gcdes, c_solver_pt, 1. ) ! lateral boundary condition 154 123 155 ! boundary condition on gcdes (only cyclic bc are required) 156 #if defined key_dynspg_fsc 157 # if defined key_mpp 158 ! Mpp: export boundary values to neighbouring processors 159 CALL lbc_lnk( gcdes, 'S', 1. ) 160 # else 161 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 162 CALL lbc_lnk( gcdes, 'T', 1. ) 163 # endif 164 #else 165 # if defined key_mpp 166 ! Mpp: export boundary values to neighbouring processors 167 CALL lbc_lnk( gcdes, 'G', 1. ) 168 # else 169 ! mono- or macro-tasking: F-point, >0, 2D array, no slab 170 CALL lbc_lnk( gcdes, 'F', 1. ) 171 # endif 172 #endif 124 ! ... gccd = matrix . gcdes 125 DO jj = 2, jpjm1 126 DO ji = fs_2, fs_jpim1 ! vector opt. 127 gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj) & 128 & +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj) & 129 & +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj) ) 130 END DO 131 END DO 132 133 ! alph = (gcr,gcr)/(gcdes,gccd) 134 radd = SUM( gcdes(:,:) * gcdmat(:,:) * gccd(:,:) ) 135 IF( lk_mpp ) CALL mpp_sum( radd ) ! sum over the global domain 136 alph = rr / radd 137 138 ! gcx = gcx + alph * gcdes 139 ! gcr = gcr - alph * gccd 140 DO jj = 2, jpjm1 141 DO ji = fs_2, fs_jpim1 ! vector opt. 142 gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) ) 143 gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) ) 144 END DO 145 END DO 173 146 174 !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,, 147 ! rnorme = (gcr,gcr) 148 rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) ) 149 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 175 150 176 ! ... gccd = matrix . gcdes177 DO jj = 2, jpjm1178 DO ji = fs_2, fs_jpim1 ! vector opt.179 gccd(ji,jj) = bmask(ji,jj)*( &180 gcdes(ji,jj) &181 +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj) &182 +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj) &183 184 END DO185 END DO151 ! test of convergence 152 IF( rnorme < epsr .OR. jn == nmax ) THEN 153 res = SQRT( rnorme ) 154 niter = jn 155 ncut = 999 156 ENDIF 157 158 ! beta = (rk+1,rk+1)/(rk,rk) 159 beta = rnorme / rr 160 rr = rnorme 186 161 187 !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,, 162 ! indicator of non-convergence or explosion 163 IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2 164 IF( ncut == 999 ) GOTO 999 165 166 ! gcdes = gcr + beta * gcdes 167 DO jj = 2, jpjm1 168 DO ji = fs_2, fs_jpim1 ! vector opt. 169 gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) ) 170 END DO 171 END DO 188 172 189 ! alph = (gcr,gcr)/(gcdes,gccd) 190 191 radd = SUM( gcdes(:,:) * gcdmat(:,:) * gccd(:,:) ) 192 193 #if defined key_mpp 194 ! Mpp: sum over all the global domain 195 CALL mpp_sum( radd ) 196 #endif 197 alph = rr / radd 198 199 !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,, 200 201 ! gcx = gcx + alph * gcdes 202 ! gcr = gcr - alph * gccd 203 DO jj = 2, jpjm1 204 DO ji = fs_2, fs_jpim1 ! vector opt. 205 gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) ) 206 gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) ) 207 END DO 208 END DO 209 210 !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,, 211 212 ! rnorme = (gcr,gcr) 213 214 rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) ) 215 216 #if defined key_mpp 217 ! Mpp: sum over all the global domain 218 CALL mpp_sum( rnorme ) 219 #endif 220 221 ! test of convergence 222 IF ( rnorme < epsr .OR. jn == nmax ) THEN 223 res = SQRT( rnorme ) 224 niter = jn 225 ncut = 999 226 ENDIF 227 228 ! beta = (rk+1,rk+1)/(rk,rk) 229 beta = rnorme / rr 230 rr = rnorme 231 232 !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,, 233 234 ! indicator of non-convergence or explosion 235 IF( jn == nmax .OR. sqrt(epsr)/eps > 1.e+20 ) kindic = -2 236 IF( ncut == 999 ) GOTO 999 237 238 ! gcdes = gcr + beta * gcdes 239 DO jj = 2, jpjm1 240 DO ji = fs_2, fs_jpim1 ! vector opt. 241 gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) ) 242 END DO 243 END DO 244 245 ! !================ 246 END DO ! End Loop 247 ! !================ 173 ! !================ 174 END DO ! End Loop 175 ! !================ 248 176 249 999 CONTINUE177 999 CONTINUE 250 178 251 179 252 ! 2.Output in gcx with lateral b.c. applied253 ! ------------------------------------------180 ! Output in gcx with lateral b.c. applied 181 ! --------------------------------------- 254 182 255 ! boundary conditions !!bug ??? 256 #if defined key_mpp 257 ! Mpp: export boundary values to neighbouring processors 258 # if defined key_dynspg_fsc 259 CALL lbc_lnk( gcx, 'S', 1. ) 260 # else 261 CALL lbc_lnk( gcx, 'G', 1. ) 262 # endif 263 #else 264 IF ( nperio /= 0 ) THEN 265 # if defined key_dynspg_fsc 266 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 267 CALL lbc_lnk( gcx, 'T', 1. ) 268 # else 269 ! mono- or macro-tasking: F-point, >0, 2D array, no slab 270 CALL lbc_lnk( gcx, 'F', 1. ) 271 # endif 272 ENDIF 273 #endif 183 CALL lbc_lnk( gcx, c_solver_pt, 1. ) 274 184 275 185 END SUBROUTINE sol_pcg
Note: See TracChangeset
for help on using the changeset viewer.
