Changeset 86 for trunk/NEMO
- Timestamp:
- 2004-04-22T15:37:51+02:00 (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SOL/solsor.F90
r16 r86 1 1 MODULE solsor 2 2 !!====================================================================== 3 !! *** MODULE solsor 3 !! *** MODULE solsor *** 4 4 !! Ocean solver : Successive Over-Relaxation solver 5 5 !!===================================================================== 6 6 7 7 !!---------------------------------------------------------------------- 8 !! sol_sor : Successive Over-Relaxation solver8 !! sol_sor : Red-Black Successive Over-Relaxation solver 9 9 !!---------------------------------------------------------------------- 10 10 !! * Modules used … … 35 35 !! ** Purpose : Solve the ellipic equation for the barotropic stream 36 36 !! function system (lk_dynspg_rl=T) or the transport divergence 37 !! system (lk_dynspg_fsc=T) using a successive-over-relaxation38 !! method.37 !! system (lk_dynspg_fsc=T) using a red-black successive-over- 38 !! relaxation method. 39 39 !! In the former case, the barotropic stream function trend has a 40 40 !! zero boundary condition along all coastlines (i.e. continent … … 55 55 !! ! 96-05 (G. Madec) merge sor and pcg formulations 56 56 !! ! 96-11 (A. Weaver) correction to preconditioning 57 !! 8.5 ! 02-08 (G. Madec) F90: Free form 58 !! ! 03-02 (C. Deltel) Red-Black SOR <== Not done yet!!! 59 !! ************* 57 !! 9.0 ! 03-04 (C. Deltel, G. Madec) Red-Black SOR in free form 60 58 !!---------------------------------------------------------------------- 61 59 !! * Arguments … … 64 62 ! ! stopped in step 65 63 ! ! set to zero before the call of solsor 66 67 64 !! * Local declarations 68 65 INTEGER :: ji, jj, jn ! dummy loop indices 66 INTEGER :: ishift 67 REAL(wp) :: ztmp, zres, zres2 68 69 INTEGER :: ijmppodd, ijmppeven 69 70 !!---------------------------------------------------------------------- 70 71 72 ijmppeven = MOD(njmpp ,2) 73 ijmppodd = MOD(njmpp+1,2) 71 74 ! ! ============== 72 75 DO jn = 1, nmax ! Iterative loop … … 77 80 ! Residus 78 81 ! ------- 79 DO jj = 2, jpjm1 80 DO ji = 2, jpim1 81 gcr(ji,jj) = gcb(ji,jj ) - gcx(ji ,jj ) & 82 - gcp(ji,jj,1) * gcx(ji ,jj-1) & 83 - gcp(ji,jj,2) * gcx(ji-1,jj ) & 84 - gcp(ji,jj,3) * gcx(ji+1,jj ) & 85 - gcp(ji,jj,4) * gcx(ji ,jj+1) 82 83 zres2 = 0.e0 84 85 ! Guess black update 86 DO jj = 2, jpjm1 87 ishift = MOD( jj-ijmppodd, 2 ) 88 DO ji = 2+ishift, jpim1, 2 89 ztmp = gcb(ji ,jj ) & 90 & - gcp(ji,jj,1) * gcx(ji ,jj-1) & 91 & - gcp(ji,jj,2) * gcx(ji-1,jj ) & 92 & - gcp(ji,jj,3) * gcx(ji+1,jj ) & 93 & - gcp(ji,jj,4) * gcx(ji ,jj+1) 94 ! Estimate of the residual 95 zres = ztmp - gcx(ji,jj) 96 zres2 = zres2 + zres * gcdmat(ji,jj) * zres 97 gcr(ji,jj) = zres * gcdmat(ji,jj) * zres 98 ! Guess update 99 gcx(ji,jj) = sor * ztmp + (1-sor) * gcx(ji,jj) 86 100 END DO 87 101 END DO 88 CALL lbc_lnk( gcr, c_solver_pt, 1. ) ! applied the lateral boubary conditions89 102 90 ! Successive over relaxation 91 DO jj = 2, jpj 92 DO ji = 1, jpi 93 gcr(ji,jj) = gcr(ji,jj) - sor * gcp(ji,jj,1) * gcr(ji,jj-1) 94 END DO 95 DO ji = 2, jpi 96 gcr(ji,jj) = gcr(ji,jj) - sor * gcp(ji,jj,2) * gcr(ji-1,jj) 103 CALL lbc_lnk( gcx, c_solver_pt, 1. ) ! applied the lateral boubary conditions 104 105 ! Guess red update 106 DO jj = 2, jpjm1 107 ishift = MOD( jj-ijmppeven, 2 ) 108 DO ji = 2+ishift, jpim1, 2 109 ztmp = gcb(ji ,jj ) & 110 & - gcp(ji,jj,1) * gcx(ji ,jj-1) & 111 & - gcp(ji,jj,2) * gcx(ji-1,jj ) & 112 & - gcp(ji,jj,3) * gcx(ji+1,jj ) & 113 & - gcp(ji,jj,4) * gcx(ji ,jj+1) 114 ! Estimate of the residual 115 zres = ztmp - gcx(ji,jj) 116 zres2 = zres2 + zres * gcdmat(ji,jj) * zres 117 gcr(ji,jj) = zres * gcdmat(ji,jj) * zres 118 ! Guess update 119 gcx(ji,jj) = sor * ztmp + (1-sor) * gcx(ji,jj) 97 120 END DO 98 121 END DO 99 100 ! gcx guess 101 DO jj = 2, jpjm1 102 DO ji = 1, jpi 103 gcx(ji,jj) = ( gcx(ji,jj) + sor * gcr(ji,jj) ) * bmask(ji,jj) 104 END DO 105 END DO 106 CALL lbc_lnk( gcx, c_solver_pt, 1. ) 107 122 123 CALL lbc_lnk( gcx, c_solver_pt, 1. ) ! applied the lateral boubary conditions 124 108 125 ! relative precision 109 rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) ) 110 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 111 126 rnorme = zres2 127 !!!! IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 128 !i 129 zres2 = MAXVAL( gcr(2:jpim1,2:jpjm1) ) 130 131 IF( lk_mpp ) CALL mpp_max( zres2 ) ! sum over the global domain 132 !i 112 133 ! test of convergence 113 134 ! old test (either res<resmax or jn=nmax) 114 135 ! IF( res < resmax .OR. jn == nmax ) THEN 136 IF( zres2 < 1.e-10 .OR. jn == nmax ) THEN 115 137 ! relative precision 116 IF( rnorme < epsr .OR. jn == nmax ) THEN 117 res = SQRT( rnorme ) 138 ! IF( rnorme < epsr .OR. jn == nmax ) THEN 139 ! res = SQRT( rnorme ) 140 res = SQRT( zres2 ) 118 141 niter = jn 119 142 ncut = 999 143 write(numout,*) 'solsor : res max = ', zres2, 'niter= ', niter 120 144 ENDIF 121 145
Note: See TracChangeset
for help on using the changeset viewer.