Changeset 16 for trunk/NEMO/OPA_SRC/SOL
- Timestamp:
- 2004-02-17T09:06:15+01:00 (20 years ago)
- Location:
- trunk/NEMO/OPA_SRC/SOL
- Files:
-
- 8 edited
-
sol_oce.F90 (modified) (4 diffs)
-
solfet.F90 (modified) (2 diffs)
-
solisl.F90 (modified) (22 diffs)
-
solisl_fdir.h90 (modified) (4 diffs)
-
solmat.F90 (modified) (14 diffs)
-
solpcg.F90 (modified) (4 diffs)
-
solsor.F90 (modified) (5 diffs)
-
solver.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SOL/sol_oce.F90
r3 r16 16 16 17 17 IMPLICIT NONE 18 PRIVATE 18 19 19 !!----------------------------------- -----------------------------------20 !!----------------------------------- 20 21 !! elliptic solver: SOR, PCG or FETI 21 !! ---------------------------------- -----------------22 INTEGER :: & !!!namsol elliptic solver / island / free surface23 nsolv = 1 , & ! = 1/2/3 type of elliptic solver24 nmax = 800 , & ! maximum of iterations for the solver25 nmisl = 4000 ! maximum pcg iterations for island22 !! ---------------------------------- 23 INTEGER , PUBLIC :: & !!: namsol elliptic solver / island / free surface 24 nsolv = 1 , & !: = 1/2/3 type of elliptic solver 25 nmax = 800 , & !: maximum of iterations for the solver 26 nmisl = 4000 !: maximum pcg iterations for island 26 27 27 REAL(wp) :: & !!!namsol elliptic solver / island / free surface28 eps = 1.e-6_wp , & ! absolute precision of the solver29 sor = 1.76_wp , & ! optimal coefficient for sor solver30 epsisl = 1.e-10_wp, & ! absolute precision on stream function solver31 rnu = 1.0_wp ! strength of the additional force used in free surface28 REAL(wp), PUBLIC :: & !!: namsol elliptic solver / island / free surface 29 eps = 1.e-6_wp , & !: absolute precision of the solver 30 sor = 1.76_wp , & !: optimal coefficient for sor solver 31 epsisl = 1.e-10_wp, & !: absolute precision on stream function solver 32 rnu = 1.0_wp !: strength of the additional force used in free surface 32 33 33 INTEGER :: &34 ncut, & ! indicator of solver convergence35 niter ! number of iteration done by the solver34 CHARACTER(len=1), PUBLIC :: & !: 35 c_solver_pt = 'T' !: nature of grid-points T (S) for free surface case 36 ! ! F (G) for rigid-lid case 36 37 37 REAL(wp) :: & 38 epsr, & ! relative precision for SOR & PCG solvers 39 epsilo, & ! precision for the FETI solver 40 rnorme, res, & ! intermediate modulus, solver residu 41 alph, & ! coefficient =(gcr,gcr)/(gcx,gccd) 42 beta, & ! coefficient =(rn+1,rn+1)/(rn,rn) 43 radd, & ! coefficient =(gccd,gcdes) 44 rr ! coefficient =(rn,rn) 38 INTEGER , PUBLIC :: & !: 39 ncut, & !: indicator of solver convergence 40 niter !: number of iteration done by the solver 45 41 46 REAL(wp), DIMENSION(jpi,jpj,4) :: & 47 gcp ! barotropic matrix extra-diagonal elements 42 REAL(wp), PUBLIC :: & !: 43 epsr, & !: relative precision for SOR & PCG solvers 44 epsilo, & !: precision for the FETI solver 45 rnorme, res, & !: intermediate modulus, solver residu 46 alph, & !: coefficient =(gcr,gcr)/(gcx,gccd) 47 beta, & !: coefficient =(rn+1,rn+1)/(rn,rn) 48 radd, & !: coefficient =(gccd,gcdes) 49 rr !: coefficient =(rn,rn) 48 50 49 REAL(wp), DIMENSION(jpi,jpj) :: & 50 gcx, gcxb, & ! now, before solution of the elliptic equation 51 gcdprc, & ! inverse diagonal preconditioning matrix 52 gcdmat, & ! diagonal preconditioning matrix 53 gcb, & ! second member of the barotropic linear system 54 gcr, & ! residu =b-a.x 55 gcdes, & ! vector descente 56 gccd ! vector such that ca.gccd=a.d (ca-1=gcdprc) 51 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,4) :: & !: 52 gcp !: barotropic matrix extra-diagonal elements 53 54 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 55 gcx, gcxb, & !: now, before solution of the elliptic equation 56 gcdprc, & !: inverse diagonal preconditioning matrix 57 gcdmat, & !: diagonal preconditioning matrix 58 gcb, & !: second member of the barotropic linear system 59 gcr, & !: residu =b-a.x 60 gcdes, & !: vector descente 61 gccd !: vector such that ca.gccd=a.d (ca-1=gcdprc) 57 62 58 63 #if defined key_feti … … 67 72 !! malistin() : concatened list of interface nodes 68 73 69 INTEGER :: nim,nxm, &74 INTEGER, PUBLIC :: nim,nxm, & 70 75 malxm,malim,malxmax,malimax, & 71 76 nifmat,njfmat,nelem,npe,matopo, & … … 85 90 madwork 86 91 87 INTEGER :: mfet(jpi*jpj+2*jpi+2*jpj+51)92 INTEGER, PUBLIC :: mfet(jpi*jpj+2*jpi+2*jpj+51) 88 93 89 REAL(wp) :: wfeti(jpj*jpi*jpi+13*jpi*jpj+19*(jpi+jpj) & 94 REAL(wp), PUBLIC :: & !: 95 wfeti(jpj*jpi*jpi+13*jpi*jpj+19*(jpi+jpj) & 90 96 +4*jpnij+33 & 91 97 +2*(jpi+jpj)*(jpnij-jpni)*jpi & … … 94 100 +3*(jpnij-jpnj+jperio)*jpj) 95 101 96 REAL(wp) :: res2, rcompt102 REAL(wp), PUBLIC :: res2, rcompt 97 103 98 104 #endif -
trunk/NEMO/OPA_SRC/SOL/solfet.F90
r3 r16 33 33 !! Solve the ellipic equation for the barotropic stream function 34 34 !! system (default option) or the transport divergence system 35 !! ( "key_dynspg_fsc") using a Finite Elements Tearing &35 !! (lk_dynspg_fsc=T) using a Finite Elements Tearing and 36 36 !! Interconnecting (FETI) approach. 37 37 !! In the former case, the barotropic stream function trend has a … … 142 142 CALL feti_vmov( noeuds, wfeti(miax), gcx ) 143 143 144 ! boundary conditions !!bug ??? check arguments... 145 # if defined key_dynspg_fsc 146 # if defined key_mpp 147 ! Mpp: export boundary values to neighbouring processors 148 CALL lbc_lnk( gcx, 'S', 1. ) 149 # else 150 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 151 IF( nperio /= 0 ) THEN 152 CALL lbc_lnk( gcx, 'T', 1. ) 153 ENDIF 154 # endif 155 # else 156 # if defined key_mpp 157 ! Mpp: export boundary values to neighbouring processors 158 CALL lbc_lnk( gcx, 'G', 1. ) 159 # else 160 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 161 IF( nperio /= 0 ) THEN 162 CALL lbc_lnk( gcx, 'F', 1. ) 163 ENDIF 164 # endif 165 # endif 144 CALL lbc_lnk( gcx, c_solver_pt, 1. ) ! lateral boundary condition 166 145 167 146 END SUBROUTINE sol_fet -
trunk/NEMO/OPA_SRC/SOL/solisl.F90
r3 r16 36 36 37 37 !! * Shared module variables 38 LOGICAL, PUBLIC :: & 39 l_isl = .TRUE. ! 'key_islands' flag 38 LOGICAL, PUBLIC, PARAMETER :: l_isl = .TRUE. !: 'key_islands' flag 40 39 41 40 !! * module variable … … 157 156 zwb(jpi,:) = 0.e0 158 157 ENDIF 159 # if defined key_mpp 160 ! Mpp: export boundary values to neighboring processors 161 CALL lbc_lnk( zwb, 'G', 1. ) 162 # endif 158 IF( lk_mpp ) CALL lbc_lnk( zwb, 'G', 1. ) 163 159 164 160 165 161 ! 1. Initialization for the search of island grid-points 166 162 ! ------------------------------------------------------ 167 # if defined key_mpp 168 169 ! Mpp : The overlap region are not taken into account 170 ! (islands bondaries are searched over subdomain only) 171 iista = 1 + jpreci 172 iiend = nlci - jpreci 173 ijsta = 1 + jprecj 174 ijend = nlcj - jprecj 175 ijstm1= 1 + jprecj 176 ijenm1= nlcj - jprecj 177 IF( nbondi == -1 .OR. nbondi == 2 ) THEN 163 164 IF( lk_mpp ) THEN 165 166 ! Mpp : The overlap region are not taken into account 167 ! (islands bondaries are searched over subdomain only) 168 iista = 1 + jpreci 169 iiend = nlci - jpreci 170 ijsta = 1 + jprecj 171 ijend = nlcj - jprecj 172 ijstm1= 1 + jprecj 173 ijenm1= nlcj - jprecj 174 IF( nbondi == -1 .OR. nbondi == 2 ) THEN 175 iista = 1 176 ENDIF 177 IF( nbondi == 1 .OR. nbondi == 2 ) THEN 178 iiend = nlci 179 ENDIF 180 IF( nbondj == -1 .OR. nbondj == 2 ) THEN 181 ijsta = 1 182 ijstm1 = 2 183 ENDIF 184 IF( nbondj == 1 .OR. nbondj == 2 ) THEN 185 ijend = nlcj 186 ijenm1 = nlcj-1 187 ENDIF 188 IF( npolj == 3 .OR. npolj == 4 ) THEN 189 ijend = nlcj-2 190 ijenm1 = nlcj-2 191 ENDIF 192 ELSE 193 ! mono- or macro-tasking environnement: full domain scan 178 194 iista = 1 179 ENDIF 180 IF( nbondi == 1 .OR. nbondi == 2 ) THEN 181 iiend = nlci 182 ENDIF 183 IF( nbondj == -1 .OR. nbondj == 2 ) THEN 195 iiend = jpi 184 196 ijsta = 1 185 197 ijstm1 = 2 186 ENDIF 187 IF( nbondj == 1 .OR. nbondj == 2 ) THEN 188 ijend = nlcj 189 ijenm1 = nlcj-1 190 ENDIF 191 IF( npolj == 3 .OR. npolj == 4 ) THEN 192 ijend = nlcj-2 193 ijenm1 = nlcj-2 194 ENDIF 195 # else 196 197 ! mono- or macro-tasking environnement: full domain scan 198 iista = 1 199 iiend = jpi 200 ijsta = 1 201 ijstm1 = 2 202 IF( nperio == 3 .OR. nperio == 4 ) THEN 203 ijend = jpj-2 204 ijenm1 = jpj-2 205 ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN 206 ijend = jpj-1 207 ijenm1 = jpj-1 208 ELSE 209 ijend = jpj 210 ijenm1 = jpj-1 211 ENDIF 212 # endif 198 IF( nperio == 3 .OR. nperio == 4 ) THEN 199 ijend = jpj-2 200 ijenm1 = jpj-2 201 ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN 202 ijend = jpj-1 203 ijenm1 = jpj-1 204 ELSE 205 ijend = jpj 206 ijenm1 = jpj-1 207 ENDIF 208 ENDIF 213 209 214 210 … … 247 243 inilt = inilt + indil(jj) 248 244 END DO 249 # if defined key_mpp 250 CALL mpp_sum( inilt ) 251 # endif 245 IF( lk_mpp ) CALL mpp_sum( inilt ) ! sum over the global domain 246 252 247 IF( inilt == 0 ) THEN 253 248 IF(lwp) THEN … … 255 250 WRITE(numout,*) ' change parameter.h' 256 251 ENDIF 257 STOP 'isldom' 252 STOP 'isldom' !cr replace by nstop 258 253 ENDIF 259 254 … … 381 376 ! Take account of redundant points 382 377 383 # if defined key_mpp 384 CALL mpp_sum( ip ) 385 # endif 378 IF( lk_mpp ) CALL mpp_sum( ip ) ! sum over the global domain 386 379 387 380 IF( ip > jpnisl ) THEN … … 391 384 WRITE(numout,*) ' change parameter.h' 392 385 ENDIF 393 STOP 'isldom' 386 STOP 'isldom' !cr => nstop 394 387 ENDIF 395 388 … … 409 402 410 403 inilt = isrchne( jpij, zwb(1,1), 1, 0. ) 411 # if defined key_mpp 412 CALL mpp_min( inilt ) 413 # endif 404 IF( lk_mpp ) CALL mpp_min( inilt ) ! min over the global domain 414 405 415 406 IF( inilt /= jpij+1 ) THEN … … 426 417 ! ---------------------------------------- 427 418 428 CALL isl pri419 CALL isl_pri 429 420 430 421 … … 432 423 ! ------------------------------------------------------- 433 424 434 CALL isl pth425 CALL isl_pth 435 426 436 427 END SUBROUTINE isl_dom … … 466 457 ipe = mnisl(3,jni) 467 458 ipw = mnisl(4,jni) 468 # if defined key_mpp 469 CALL mpp_sum( ip )470 CALL mpp_sum( ipn )471 CALL mpp_sum( ips )472 CALL mpp_sum( ipe )473 CALL mpp_sum( ipw )474 # endif 459 IF( lk_mpp ) THEN 460 CALL mpp_sum( ip ) ! sums over the global domain 461 CALL mpp_sum( ipn ) 462 CALL mpp_sum( ips ) 463 CALL mpp_sum( ipe ) 464 CALL mpp_sum( ipw ) 465 ENDIF 475 466 IF(lwp) THEN 476 467 WRITE(numout,9000) jni … … 484 475 END DO 485 476 486 ! FORMAT 487 477 ! FORMAT !!cr => no more format 488 478 9000 FORMAT(/, /, 'island number= ', i2 ) 489 479 9010 FORMAT(/, 'npil=',i4,' npn=',i3,' nps=',i3,' npe=',i3,' npw=',i3 ) … … 514 504 !!---------------------------------------------------------------------- 515 505 !! * Local declarations 516 INTEGER :: j i, jj, jni, jii, jnp ! dummy loop indices517 INTEGER :: ii, ij 506 INTEGER :: jni, jii, jnp ! dummy loop indices 507 INTEGER :: ii, ij ! temporary integers 518 508 !!---------------------------------------------------------------------- 519 509 … … 587 577 REAL(wp), DIMENSION(jpi,jpj) :: zlamt, zphit 588 578 REAL(wp), DIMENSION(jpi,jpj,2) :: zwx 589 # if defined key_mpp590 579 REAL(wp), DIMENSION(jpisl*jpisl) :: ztab 591 # endif592 580 !!---------------------------------------------------------------------- 593 581 … … 674 662 675 663 END DO 676 # if defined key_mpp 677 DO jnj=1,jpisl 678 DO jni=1,jpisl 679 ztab(jni+(jnj-1)*jpisl) = aisl(jni,jnj) 680 END DO 681 END DO 682 683 CALL mpp_sum( ztab, jpisl*jpisl ) 684 !! CALL mpp_sum( aisl, jpisl*jpisl ) 685 # endif 664 IF( lk_mpp ) THEN 665 DO jnj = 1, jpisl 666 DO jni = 1, jpisl 667 ztab(jni+(jnj-1)*jpisl) = aisl(jni,jnj) 668 END DO 669 END DO 670 CALL mpp_sum( ztab, jpisl*jpisl ) ! sum over the global domain 671 !! CALL mpp_sum( aisl, jpisl*jpisl ) 672 ENDIF 686 673 687 674 ! 1.3 Control print … … 775 762 REAL(wp) :: zep(jpisl), zlamt(jpi,jpj), zphit(jpi,jpj), zdept(1), zprec(4) 776 763 REAL(wp) :: zdate0, zdt 777 # if defined key_mpp 764 REAL(wp) :: t2p1(jpi,1,1) 778 765 INTEGER :: iloc 779 # endif780 766 !!---------------------------------------------------------------------- 781 767 … … 907 893 ! Right hand side of the streamfunction equation 908 894 909 # if defined key_mpp 910 911 ! north fold treatment 912 IF( npolj == 3 .OR. npolj == 5) iloc=jpiglo-(nimpp-1+nimppt(nono+1)-1) 913 IF( npolj == 4 .OR. npolj == 6) iloc=jpiglo-2*(nimpp-1) 914 t2p1(:,1,1) = 0.e0 915 ! north and south grid-points 916 DO jii = 1, 2 917 DO jnp = 1, mnisl(jii,jni) 918 ii = miisl(jnp,jii,jni) 919 ij = mjisl(jnp,jii,jni) 920 IF( ( npolj == 3 .OR. npolj == 4 ) .AND. & 921 ( ij == nlcj-1 .AND. jii == 1) ) THEN 922 iju=iloc-ii+1 923 t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 924 ELSEIF( ( npolj == 5 .OR. npolj == 6 ) .AND. & 925 ( ij == nlcj-1 .AND. jii == 1) ) THEN 926 iju=iloc-ii 927 gcb(ii,ij) = gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 928 t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 929 ELSE 930 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 931 ENDIF 932 END DO 933 END DO 934 935 ! east and west grid-points 936 937 DO jii = 3, 4 938 DO jnp = 1, mnisl(jii,jni) 939 ii = miisl(jnp,jii,jni) 940 ij = mjisl(jnp,jii,jni) 941 gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 942 END DO 943 END DO 944 CALL mpplnks( gcb ) 945 946 # else 947 948 ! north and south grid-points 949 DO jii = 1, 2 950 DO jnp = 1, mnisl(jii,jni) 951 ii = miisl(jnp,jii,jni) 952 ij = mjisl(jnp,jii,jni) 953 IF( ( nperio == 3 .OR. nperio == 4 ) .AND. & 954 ( ij == jpj-1 .AND. jii == 1) ) THEN 955 gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 956 ELSEIF( ( nperio == 5 .OR. nperio == 6 ) .AND. & 957 ( ij == jpj-1 .AND. jii == 1) ) THEN 958 gcb(ii,ij) = gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 959 gcb(jpi-ii,ij) = gcb(jpi-ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 960 ELSE 961 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 962 ENDIF 963 END DO 964 END DO 965 966 ! east and west grid-points 967 DO jii = 3, 4 968 DO jnp = 1, mnisl(jii,jni) 969 ii = miisl(jnp,jii,jni) 970 ij = mjisl(jnp,jii,jni) 971 IF( bmask(ii-jii+3,ij) /= 0. ) THEN 895 IF( lk_mpp ) THEN 896 897 ! north fold treatment 898 IF( npolj == 3 .OR. npolj == 5) iloc=jpiglo-(nimpp-1+nimppt(nono+1)-1) 899 IF( npolj == 4 .OR. npolj == 6) iloc=jpiglo-2*(nimpp-1) 900 t2p1(:,1,1) = 0.e0 901 ! north and south grid-points 902 DO jii = 1, 2 903 DO jnp = 1, mnisl(jii,jni) 904 ii = miisl(jnp,jii,jni) 905 ij = mjisl(jnp,jii,jni) 906 IF( ( npolj == 3 .OR. npolj == 4 ) .AND. & 907 ( ij == nlcj-1 .AND. jii == 1) ) THEN 908 iju=iloc-ii+1 909 t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 910 ELSEIF( ( npolj == 5 .OR. npolj == 6 ) .AND. & 911 ( ij == nlcj-1 .AND. jii == 1) ) THEN 912 iju=iloc-ii 913 gcb(ii,ij) = gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 914 t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 915 ELSE 916 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 917 ENDIF 918 END DO 919 END DO 920 921 ! east and west grid-points 922 923 DO jii = 3, 4 924 DO jnp = 1, mnisl(jii,jni) 925 ii = miisl(jnp,jii,jni) 926 ij = mjisl(jnp,jii,jni) 972 927 gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 973 ELSE 974 975 ! east-west cyclic boundary conditions 976 IF( ii-jii+3 == 1 ) THEN 977 gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 928 END DO 929 END DO 930 931 IF( lk_mpp ) CALL mpplnks( gcb ) !!bug ? should use an lbclnk ? is it possible? 932 933 ELSE 934 935 ! north and south grid-points 936 DO jii = 1, 2 937 DO jnp = 1, mnisl(jii,jni) 938 ii = miisl(jnp,jii,jni) 939 ij = mjisl(jnp,jii,jni) 940 IF( ( nperio == 3 .OR. nperio == 4 ) .AND. & 941 ( ij == jpj-1 .AND. jii == 1) ) THEN 942 gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 943 ELSEIF( ( nperio == 5 .OR. nperio == 6 ) .AND. & 944 ( ij == jpj-1 .AND. jii == 1) ) THEN 945 gcb(ii,ij) = gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 946 gcb(jpi-ii,ij) = gcb(jpi-ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 947 ELSE 948 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 978 949 ENDIF 979 ENDIF 980 END DO 981 END DO 982 983 # endif 950 END DO 951 END DO 952 953 ! east and west grid-points 954 DO jii = 3, 4 955 DO jnp = 1, mnisl(jii,jni) 956 ii = miisl(jnp,jii,jni) 957 ij = mjisl(jnp,jii,jni) 958 IF( bmask(ii-jii+3,ij) /= 0. ) THEN 959 gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 960 ELSE 961 ! east-west cyclic boundary conditions 962 IF( ii-jii+3 == 1 ) THEN 963 gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 964 ENDIF 965 ENDIF 966 END DO 967 END DO 968 ENDIF 984 969 985 970 ! Preconditioned right hand side and absolute precision … … 1011 996 END DO 1012 997 END DO 1013 # if defined key_mpp 1014 CALL mpp_sum( rnorme ) 1015 # endif 998 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 999 1016 1000 IF(lwp) WRITE(numout,*) 'rnorme ', rnorme 1017 1001 epsr = epsisl * epsisl * rnorme … … 1070 1054 END DO 1071 1055 ENDIF 1072 # if defined key_mpp 1073 CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. ) 1074 # endif 1056 IF( lk_mpp ) CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. ) ! link at G-point 1075 1057 1076 1058 … … 1212 1194 END DO 1213 1195 END DO 1214 # if defined key_mpp 1215 ! Mpp : global sum to obtain global dot from local ones 1216 CALL mpp_sum( bisl, jpisl ) 1217 # endif 1196 IF( lk_mpp ) CALL mpp_sum( bisl, jpisl ) ! sum over the global domain 1197 1218 1198 DO jni = 1, jpisl ! Island stream function trend 1219 1199 visl(jni) = 0.e0 … … 1270 1250 zfact = 1.e-6 * bsfn(miisl(1,0,jni),mjisl(1,0,jni)) 1271 1251 ENDIF 1272 # if defined key_mpp 1273 CALL mpp_isl( zfact ) 1274 # endif 1252 IF( lk_mpp ) CALL mpp_isl( zfact ) 1253 1275 1254 IF(lwp) WRITE(numisp,9300) kt, jni, zfact, visl(jni) 1276 1255 IF( MOD( kt, nwrite ) == 0 .OR. kindic < 0 & … … 1293 1272 !! Default option Empty module 1294 1273 !!---------------------------------------------------------------------- 1295 LOGICAL, PUBLIC :: l_isl = .FALSE. !'key_islands' flag1274 LOGICAL, PUBLIC, PARAMETER :: l_isl = .FALSE. !: 'key_islands' flag 1296 1275 CONTAINS 1297 1276 SUBROUTINE isl_dom ! Empty routine … … 1304 1283 END SUBROUTINE isl_dyn_spg 1305 1284 SUBROUTINE isl_stp_ctl( kt, kindic ) ! Empty routine 1306 WRITE(*,*) kt, kindic ! no compilation warning1285 WRITE(*,*) 'isl_stp_ctl: You should not have seen this print! error?', kt, kindic 1307 1286 END SUBROUTINE isl_stp_ctl 1308 1287 #endif -
trunk/NEMO/OPA_SRC/SOL/solisl_fdir.h90
r3 r16 109 109 110 110 END DO 111 # if defined key_mpp 112 CALL mpp_sum( aisl, jpisl*jpisl ) 113 # endif 111 IF( lk_mpp ) CALL mpp_sum( aisl, jpisl*jpisl ) ! sum over the global domain 114 112 115 113 ! 1.3 Control print 116 117 114 IF(lwp) THEN 118 115 WRITE(numout,*) … … 296 293 ! 1.2 Right hand side of the stream FUNCTION equation 297 294 298 # if defined key_mpp 299 300 ! north fold treatment 301 IF( npolj == 3 ) iloc = jpiglo -(nimpp-1+nimppt(nono+1)-1) 302 IF( npolj == 4 ) iloc = jpiglo - 2*(nimpp-1) 303 t2p1(:,1,1) = 0. 304 ! north and south grid-points 305 DO jii = 1, 2 306 DO jnp = 1, mnisl(jii,jni) 307 ii = miisl(jnp,jii,jni) 308 ij = mjisl(jnp,jii,jni) 309 IF( ( npolj == 3 .OR. npolj == 4 ) .AND. ( ij == nlcj-1 .AND. jii == 1) ) THEN 310 iju=iloc-ii+1 311 t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 312 ELSE 313 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 314 ENDIF 315 END DO 316 END DO 317 318 ! east and west grid-points 319 320 DO jii = 3, 4 321 DO jnp = 1, mnisl(jii,jni) 322 ii = miisl(jnp,jii,jni) 323 ij = mjisl(jnp,jii,jni) 324 gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 325 END DO 326 END DO 327 CALL mpplnks( gcb ) 328 329 # else 330 331 ! north and south grid-points 332 DO jii = 1, 2 333 DO jnp = 1, mnisl(jii,jni) 334 ii = miisl(jnp,jii,jni) 335 ij = mjisl(jnp,jii,jni) 336 IF( ( nperio == 3 .OR. nperio == 4 ) .AND. ( ij == jpj-1 .AND. jii == 1) ) THEN 337 gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 338 ELSE 339 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 340 ENDIF 341 END DO 342 END DO 343 344 ! east and west grid-points 345 DO jii = 3, 4 346 DO jnp = 1, mnisl(jii,jni) 347 ii = miisl(jnp,jii,jni) 348 ij = mjisl(jnp,jii,jni) 349 IF( bmask(ii-jii+3,ij) /= 0. ) THEN 295 IF( lk_mpp ) THEN 296 ! north fold treatment 297 IF( npolj == 3 ) iloc = jpiglo -(nimpp-1+nimppt(nono+1)-1) 298 IF( npolj == 4 ) iloc = jpiglo - 2*(nimpp-1) 299 t2p1(:,1,1) = 0. 300 ! north and south grid-points 301 DO jii = 1, 2 302 DO jnp = 1, mnisl(jii,jni) 303 ii = miisl(jnp,jii,jni) 304 ij = mjisl(jnp,jii,jni) 305 IF( ( npolj == 3 .OR. npolj == 4 ) .AND. ( ij == nlcj-1 .AND. jii == 1) ) THEN 306 iju=iloc-ii+1 307 t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 308 ELSE 309 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 310 ENDIF 311 END DO 312 END DO 313 314 ! east and west grid-points 315 316 DO jii = 3, 4 317 DO jnp = 1, mnisl(jii,jni) 318 ii = miisl(jnp,jii,jni) 319 ij = mjisl(jnp,jii,jni) 350 320 gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 351 ELSE 352 ! east-west cyclic boundary conditions 353 IF( ii-jii+3 == 1 ) THEN 354 gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 321 END DO 322 END DO 323 324 IF( lk_mpp ) CALL mpplnks( gcb ) !!bug ? should use an lbclnk ? is it possible??? 325 326 ELSE 327 ! north and south grid-points 328 DO jii = 1, 2 329 DO jnp = 1, mnisl(jii,jni) 330 ii = miisl(jnp,jii,jni) 331 ij = mjisl(jnp,jii,jni) 332 IF( ( nperio == 3 .OR. nperio == 4 ) .AND. ( ij == jpj-1 .AND. jii == 1) ) THEN 333 gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 334 ELSE 335 gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 355 336 ENDIF 356 ENDIF 357 END DO 358 END DO 359 360 # endif 337 END DO 338 END DO 339 340 ! east and west grid-points 341 DO jii = 3, 4 342 DO jnp = 1, mnisl(jii,jni) 343 ii = miisl(jnp,jii,jni) 344 ij = mjisl(jnp,jii,jni) 345 IF( bmask(ii-jii+3,ij) /= 0. ) THEN 346 gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 347 ELSE 348 ! east-west cyclic boundary conditions 349 IF( ii-jii+3 == 1 ) THEN 350 gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) 351 ENDIF 352 ENDIF 353 END DO 354 END DO 355 ENDIF 361 356 362 357 ! 1.4 Preconditioned right hand side and absolute precision … … 388 383 END DO 389 384 END DO 390 #if defined key_mpp 391 CALL mpp_sum( rnorme ) 392 #endif 385 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 386 393 387 IF(lwp) WRITE(numout,*) 'rnorme ', rnorme 394 388 epsr = epsisl * epsisl * rnorme … … 451 445 END DO 452 446 ENDIF 453 #if defined key_mpp 454 CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. ) 455 #endif 447 IF( lk_mpp ) CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. ) ! link at G-point 456 448 457 449 -
trunk/NEMO/OPA_SRC/SOL/solmat.F90
r3 r16 18 18 USE obc_oce ! ocean open boundary conditions 19 19 USE lib_mpp ! distributed memory computing 20 USE dynspg_rl 21 USE dynspg_fsc 20 22 21 23 IMPLICIT NONE … … 38 40 !! 39 41 !! ** Method : The matrix depends on the type of free surface: 40 !! * default option: rigid lid and bsf42 !! * lk_dynspg_rl=T: rigid lid formulation 41 43 !! The matrix is built for the barotropic stream function system. 42 44 !! a diagonal preconditioning matrix is also defined. 43 !! * 'key_dynspg_fsc' defined: free surface45 !! * lk_dynspg_fsc=T: free surface formulation 44 46 !! The matrix is built for the divergence of the transport system 45 47 !! a diagonal preconditioning matrix is also defined. … … 67 69 INTEGER :: ii, ij, iiend, ijend ! temporary integers 68 70 REAL(wp) :: zcoefs, zcoefw, zcoefe, zcoefn ! temporary scalars 69 REAL(wp) :: z2dt 70 #if defined key_dynspg_fsc 71 REAL(wp) :: zcoef 72 #endif 71 REAL(wp) :: z2dt, zcoef 73 72 !!---------------------------------------------------------------------- 74 73 … … 83 82 84 83 ! initialize to zero 84 zcoef = 0.e0 85 85 gcp(:,:,1) = 0.e0 86 86 gcp(:,:,2) = 0.e0 … … 94 94 95 95 #if defined key_dynspg_fsc && ! defined key_obc 96 !!cr IF( lk_dynspg_fsc .AND. .NOT.lk_obc ) THEN 96 97 97 98 ! defined the coefficients for free surface elliptic system … … 99 100 DO jj = 2, jpjm1 100 101 DO ji = 2, jpim1 101 zcoef = z2dt * z2dt * g * rnu * bmask(ji,jj)102 zcoef = z2dt * z2dt * grav * rnu * bmask(ji,jj) 102 103 zcoefs = -zcoef * hv(ji ,jj-1) * e1v(ji ,jj-1) / e2v(ji ,jj-1) ! south coefficient 103 104 zcoefw = -zcoef * hu(ji-1,jj ) * e2u(ji-1,jj ) / e1u(ji-1,jj ) ! west coefficient … … 109 110 gcp(ji,jj,4) = zcoefn 110 111 gcdmat(ji,jj) = e1t(ji,jj) * e2t(ji,jj) * bmask(ji,jj) & ! diagonal coefficient 111 112 & - zcoefs -zcoefw -zcoefe -zcoefn 112 113 END DO 113 114 END DO 114 115 115 116 # elif defined key_dynspg_fsc && defined key_obc 117 !!cr ELSEIF( lk_dynspg_fsc .AND. lk_obc ) THEN 116 118 117 119 ! defined gcdmat in the case of open boundaries … … 119 121 DO jj = 2, jpjm1 120 122 DO ji = 2, jpim1 121 zcoef = z2dt * z2dt * g * rnu * bmask(ji,jj)123 zcoef = z2dt * z2dt * grav * rnu * bmask(ji,jj) 122 124 ! south coefficient 123 125 IF( ( lpsouthobc ) .AND. ( jj == njs0p1 ) ) THEN … … 159 161 160 162 # else 163 !!cr ELSE 161 164 162 165 ! defined the coefficients for bsf elliptic system … … 173 176 gcp(ji,jj,4) = zcoefn 174 177 gcdmat(ji,jj) = -zcoefs -zcoefw -zcoefe -zcoefn ! diagonal coefficient 175 176 178 END DO 177 179 END DO 178 180 181 !!cr ENDIF 179 182 #endif 180 183 … … 194 197 ! account for the existence of the south symmetric bassin. 195 198 199 !!cr IF( .NOT.lk_dynspg_fsc ) THEN 196 200 #if ! defined key_dynspg_fsc 197 201 IF( nperio == 2 ) THEN … … 203 207 END DO 204 208 ENDIF 209 !!cr ENDIF 205 210 #endif 206 211 … … 225 230 gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:) 226 231 gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:) 232 IF( nsolv == 2 ) gccd(:,:) = sor * gcp(:,:,2) 227 233 228 234 ELSE … … 467 473 nnitot = nni 468 474 469 CALL mpp_sum( nnitot,1,numit0ete)475 CALL mpp_sum( nnitot, 1, numit0ete ) 470 476 CALL feti_creadr(malxm,malxmax,nxm,npe*npe,maae,'ae') 471 477 -
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 -
trunk/NEMO/OPA_SRC/SOL/solsor.F90
r3 r16 22 22 !! * Routine accessibility 23 23 PUBLIC sol_sor ! ??? 24 24 25 !!---------------------------------------------------------------------- 25 26 !! OPA 9.0 , LODYC-IPSL (2003) … … 28 29 CONTAINS 29 30 30 SUBROUTINE sol_sor( k t, kindic )31 SUBROUTINE sol_sor( kindic ) 31 32 !!---------------------------------------------------------------------- 32 33 !! *** ROUTINE sol_sor *** 33 34 !! 34 35 !! ** Purpose : Solve the ellipic equation for the barotropic stream 35 !! function system ( default option) or the transport divergence36 !! system ( key_dynspg_fsc =T) using a successive-over-relaxation36 !! function system (lk_dynspg_rl=T) or the transport divergence 37 !! system (lk_dynspg_fsc=T) using a successive-over-relaxation 37 38 !! method. 38 39 !! In the former case, the barotropic stream function trend has a … … 59 60 !!---------------------------------------------------------------------- 60 61 !! * Arguments 61 INTEGER, INTENT( in ) :: kt ! ocean time-step62 62 INTEGER, INTENT( inout ) :: kindic ! solver indicator, < 0 if the conver- 63 63 ! ! gence is not reached: the model is … … 67 67 !! * Local declarations 68 68 INTEGER :: ji, jj, jn ! dummy loop indices 69 REAL(wp) :: zgwgt ! temporary scalar70 69 !!---------------------------------------------------------------------- 71 70 72 73 ! Iterative loop 74 ! ============== 75 76 IF( kt == nit000 ) gccd(:,:) = sor * gcp(:,:,2) 71 ! ! ============== 72 DO jn = 1, nmax ! Iterative loop 73 ! ! ============== 77 74 78 79 DO jn = 1, nmax 80 81 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 75 CALL lbc_lnk( gcx, c_solver_pt, 1. ) ! applied the lateral boubary conditions 82 76 83 ! boundary conditions (at each sor iteration) only cyclic b. c. are required 84 #if defined key_dynspg_fsc 85 # if defined key_mpp 86 ! Mpp: export boundary values to neighbouring processors 87 CALL lbc_lnk( gcx, 'S', 1. ) ! S=T with special staff ??? 88 # else 89 CALL lbc_lnk( gcx, 'T', 1. ) 90 # endif 91 #else 92 # if defined key_mpp 93 ! Mpp: export boundary values to neighbouring processors 94 CALL lbc_lnk( gcx, 'G', 1. ) ! G= F with special staff ??? 95 # else 96 CALL lbc_lnk( gcx, 'F', 1. ) 97 # endif 98 #endif 99 100 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 101 102 ! 1. Residus 103 ! ---------- 104 77 ! Residus 78 ! ------- 105 79 DO jj = 2, jpjm1 106 80 DO ji = 2, jpim1 107 gcr(ji,jj) = gcb(ji,jj ) - gcx(ji ,jj ) &108 - gcp(ji,jj,1)*gcx(ji ,jj-1) &109 - gcp(ji,jj,2)*gcx(ji-1,jj ) &110 - gcp(ji,jj,3)*gcx(ji+1,jj ) &111 - gcp(ji,jj,4)*gcx(ji ,jj+1)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) 112 86 END DO 113 87 END DO 88 CALL lbc_lnk( gcr, c_solver_pt, 1. ) ! applied the lateral boubary conditions 114 89 115 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 116 117 ! 1.1 Boundary conditions (at each sor iteration) only cyclic b. c. are required 118 #if defined key_dynspg_fsc 119 # if defined key_mpp 120 ! Mpp: export boundary values to neighbouring processors 121 CALL lbc_lnk( gcr, 'S', 1. ) 122 # else 123 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 124 CALL lbc_lnk( gcr, 'T', 1. ) 125 # endif 126 #else 127 # if defined key_mpp 128 ! Mpp: export boundary values to neighbouring processors 129 CALL lbc_lnk( gcr, 'G', 1. ) 130 # else 131 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 132 CALL lbc_lnk( gcr, 'F', 1. ) 133 # endif 134 #endif 135 136 ! 1.2 Successive over relaxation 137 90 ! Successive over relaxation 138 91 DO jj = 2, jpj 139 92 DO ji = 1, jpi 140 gcr(ji,jj) = gcr(ji,jj) - sor *gcp(ji,jj,1)*gcr(ji,jj-1)93 gcr(ji,jj) = gcr(ji,jj) - sor * gcp(ji,jj,1) * gcr(ji,jj-1) 141 94 END DO 142 95 DO ji = 2, jpi 143 gcr(ji,jj) = gcr(ji,jj) - sor *gcp(ji,jj,2)*gcr(ji-1,jj)96 gcr(ji,jj) = gcr(ji,jj) - sor * gcp(ji,jj,2) * gcr(ji-1,jj) 144 97 END DO 145 98 END DO 146 99 147 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,148 149 100 ! gcx guess 150 151 101 DO jj = 2, jpjm1 152 102 DO ji = 1, jpi 153 gcx(ji,jj) = ( gcx(ji,jj)+sor*gcr(ji,jj))*bmask(ji,jj)103 gcx(ji,jj) = ( gcx(ji,jj) + sor * gcr(ji,jj) ) * bmask(ji,jj) 154 104 END DO 155 105 END DO 156 157 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 158 159 ! boundary conditions (at each sor iteration) only cyclic b. c. are required 160 #if defined key_dynspg_fsc 161 # if defined key_mpp 162 ! Mpp: export boundary values to neighbouring processors 163 CALL lbc_lnk( gcx, 'S', 1. ) 164 # else 165 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 166 CALL lbc_lnk( gcx, 'T', 1. ) 167 # endif 168 #else 169 # if defined key_mpp 170 ! Mpp: export boundary values to neighbouring processors 171 CALL lbc_lnk( gcx, 'G', 1. ) 172 # else 173 ! mono- or macro-tasking: W-point, >0, 2D array, no slab 174 CALL lbc_lnk( gcx, 'F', 1. ) 175 # endif 176 #endif 177 178 ! maximal residu (old exit test on the maximum value of residus) 179 ! 180 ! imax = isamax( jpi*jpj, gcr, 1 ) 181 182 ! avoid an out of bound in no bounds compilation 183 184 ! iimax1 = mod( imax, jpi ) 185 ! ijmax1 = int( float(imax) / float(jpi)) + 1 186 ! resmax = abs( gcr(iimax1,ijmax1) ) 106 CALL lbc_lnk( gcx, c_solver_pt, 1. ) 187 107 188 108 ! relative precision 189 190 rnorme = 0. 191 DO jj = 1, jpj 192 DO ji = 1, jpi 193 zgwgt = gcdmat(ji,jj) * gcr(ji,jj) 194 rnorme= rnorme + gcr(ji,jj)*zgwgt 195 END DO 196 END DO 197 198 #if defined key_mpp 199 ! mpp sum over all the global domain 200 CALL mpp_sum( rnorme ) 201 #endif 109 rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) ) 110 IF( lk_mpp ) CALL mpp_sum( rnorme ) ! sum over the global domain 202 111 203 112 ! test of convergence … … 216 125 !**** 217 126 218 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,219 220 127 ! indicator of non-convergence or explosion 221 222 128 IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2 223 129 IF( ncut == 999 ) GOTO 999 224 130 225 226 ! END of iterative loop 227 ! ===================== 228 229 END DO 230 131 ! ! ===================== 132 END DO ! END of iterative loop 133 ! ! ===================== 231 134 232 135 999 CONTINUE 233 136 234 137 235 ! 2.Output in gcx236 ! ------------- ----138 ! Output in gcx 139 ! ------------- 237 140 238 ! boundary conditions (est-ce necessaire? je ne crois pas!!!!) 239 240 #if defined key_dynspg_fsc 241 # if defined key_mpp 242 ! Mpp: export boundary values to neighbouring processors 243 CALL lbc_lnk( gcx, 'S', 1. ) 244 # else 245 IF( nperio /= 0 ) THEN 246 CALL lbc_lnk( gcx, 'T', 1. ) 247 ENDIF 248 # endif 249 #else 250 # if defined key_mpp 251 ! Mpp: export boundary values to neighbouring processors 252 CALL lbc_lnk( gcx, 'G', 1. ) 253 # else 254 IF( nperio /= 0 ) THEN 255 CALL lbc_lnk( gcx, 'F', 1. ) 256 ENDIF 257 # endif 258 #endif 141 CALL lbc_lnk( gcx, c_solver_pt, 1. ) ! boundary conditions (est-ce necessaire? je ne crois pas!!!!) 259 142 260 143 END SUBROUTINE sol_sor -
trunk/NEMO/OPA_SRC/SOL/solver.F90
r3 r16 18 18 USE in_out_manager ! I/O manager 19 19 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 20 USE lib_mpp 21 USE dynspg_rl 22 USE dynspg_fsc 20 23 21 24 IMPLICIT NONE … … 34 37 !! * default option: barotropic stream function system 35 38 !! and islands initialization (if l_isl=T) 36 !! * key_dynspg_fsc = T : transport divergence system. No specific39 !! * lk_dynspg_fsc = T : transport divergence system. No specific 37 40 !! treatment of islands. 38 41 !! 39 42 !! ** Method : 40 43 !! - Compute the local depth of the water column at u- and v-point 41 !! ( key_dynspg_fsc = T) or its inverse (key_dynspg_rl = T).44 !! (lk_dynspg_fsc = T) or its inverse (lk_dynspg_rl = T). 42 45 !! The local depth of the water column is computed by summing 43 46 !! the vertical scale factors. For its inverse, the thickness of … … 56 59 !! 57 60 !! ** Action : - hur, hvr : masked inverse of the local depth at 58 !! u- and v-point. ( key_dynspg_rl = T)61 !! u- and v-point. (lk_dynspg_rl = T) 59 62 !! - hu, hv : masked local depth at u- and v- points 60 !! (key_dynspg_fsc = T) 63 !! (lk_dynspg_fsc = T) 64 !! - c_solver_pt : nature of the gridpoint at which the 65 !! solver is applied 61 66 !! References : 62 67 !! Jensen, 1986: adv. phys. oceanogr. num. mod.,ed. o brien,87-110. … … 115 120 ENDIF 116 121 117 #if defined key_dynspg_fsc 122 IF( lk_dynspg_fsc ) THEN 118 123 IF(lwp) WRITE(numout,*) 119 124 IF(lwp) WRITE(numout,*) ' *** free surface formulation' … … 123 128 nstop = nstop + 1 124 129 ENDIF 125 #endif 126 #if defined key_dynspg_rl 130 ELSEIF( lk_dynspg_rl ) THEN 127 131 IF(lwp) WRITE(numout,*) 128 132 IF(lwp) WRITE(numout,*) ' *** Rigid lid formulation' 129 #endif 130 #if defined key_dynspg_fsc && defined key_dynspg_rl 133 ELSE 134 IF(lwp) WRITE(numout,cform_err) 135 IF(lwp) WRITE(numout,*) ' Chose at least one surface pressure gradient calculation: free surface or rigid-lid' 136 nstop = nstop + 1 137 ENDIF 138 IF( lk_dynspg_fsc .AND. lk_dynspg_rl ) THEN 131 139 IF(lwp) WRITE(numout,cform_err) 132 140 IF(lwp) WRITE(numout,*) ' Chose between free surface or rigid-lid, not both' 133 141 nstop = nstop + 1 134 #endif 142 ENDIF 135 143 136 144 SELECT CASE ( nsolv ) … … 144 152 CASE ( 3 ) ! FETI solver 145 153 IF(lwp) WRITE(numout,*) ' use the FETI solver' 146 #if ! defined key_mpp 147 IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp option' 148 nstop = nstop + 1 149 #else 150 IF( jpnij == 1 ) THEN 151 IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor' 154 IF( .NOT.lk_mpp ) THEN 155 IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option' 152 156 nstop = nstop + 1 153 ENDIF 154 #endif 157 ELSE 158 IF( jpnij == 1 ) THEN 159 IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor' 160 nstop = nstop + 1 161 ENDIF 162 ENDIF 155 163 156 164 CASE DEFAULT … … 161 169 END SELECT 162 170 171 ! Grid-point at which the solver is applied 172 ! ----------------------------------------- 173 174 IF( lk_dynspg_rl ) THEN ! rigid-lid 175 IF( lk_mpp ) THEN 176 c_solver_pt = 'G' ! G= F with special staff ??? which one? 177 ELSE 178 c_solver_pt = 'F' 179 ENDIF 180 ELSE ! free surface T-point 181 IF( lk_mpp ) THEN 182 c_solver_pt = 'S' ! S=T with special staff ??? which one? 183 ELSE 184 c_solver_pt = 'T' 185 ENDIF 186 ENDIF 187 163 188 164 189 ! Construction of the elliptic system matrix
Note: See TracChangeset
for help on using the changeset viewer.
