!!---------------------------------------------------------------------- !! *** solisl_fdir.h90 *** !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! islmat : Compute the matrix associated with islands !! and save it in a direct access file !! islbsf : Compute the barotropic streamfunction of each island !! and save it in a direct access file !!---------------------------------------------------------------------- SUBROUTINE islmat !!---------------------------------------------------------------------- !! *** ROUTINE islmat *** !! !! ** Purpose : Compute and invert the island matrix !! !! ** Method : aisl(jni,jnj) matrix constituted by bsfisl circulation !! !! ** Action : - aisl : island matrix !! - aislm1 : invert of the island matrix !! !! ** file : islands: contain bsfisl, aisl and aislm1 !! !! History : !! 1.0 ! 88-10 (G. Madec) Original code !! 7.0 ! 96-01 (G. Madec) suppression of common work arrays !! 8.5 ! 02-08 (G. Madec) Free form, F90 !!---------------------------------------------------------------------- !! * Local declarations INTEGER :: ji, jj, jni, jnj, ios, jl REAL(wp) :: zwx(jpi,jpj,2) !!---------------------------------------------------------------------- ! I. Island matrix lecture in numisl (if it exists) ! ================================== ! file already open in islbsf routine ! Lecture READ( numisl, REC=jpisl+2, IOSTAT=ios ) aisl, aislm1 IF( ios /= 0 ) GOTO 110 ! Control print IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' islmat: lecture aisl/aislm1 in numisl done' WRITE(numout,*) ' ~~~~~~' WRITE(numout,*) WRITE(numout,*) ' island matrix' WRITE(numout,*) ' -------------' WRITE(numout,*) DO jnj = 1, jpisl WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl ) END DO WRITE(numout,*) WRITE(numout,*) ' inverse of the island matrix' WRITE(numout,*) ' ----------------------------' WRITE(numout,*) DO jnj = 1, jpisl WRITE(numout,'(12e11.3)') ( aislm1(jni,jnj), jni=1, jpisl ) END DO ENDIF CLOSE(numisl) RETURN 110 CONTINUE ! II. Island matrix computation ! ============================= DO jnj = 1, jpisl ! 1.1 Circulation of bsf(jnj) around island jni DO jj = 2, jpj zwx(:,jj,1) = -( bsfisl(:,jj,jnj) - bsfisl(:,jj-1,jnj) ) END DO DO ji = 1, jpi zwx(ji,1,1) = 0.e0 END DO DO jj = 1, jpj DO ji = 2, jpi zwx(ji,jj,2) = ( bsfisl(ji,jj,jnj) - bsfisl(ji-1,jj,jnj) ) END DO END DO DO jj = 1, jpj zwx(1,jj,2) = 0.e0 END DO ! 1.2 Island matrix DO jni = 1, jpisl aisl(jni,jnj) = 0. DO jl = 1, 2 DO jj = 1, jpj DO ji = 1, jpi aisl(jni,jnj) = aisl(jni,jnj) + acisl2(ji,jj,jl,jni) * zwx(ji,jj,jl) END DO END DO END DO END DO END DO # if defined key_mpp CALL mpp_sum( aisl, jpisl*jpisl ) # endif ! 1.3 Control print IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' islmat: island matrix' WRITE(numout,*) ' ~~~~~~' WRITE(numout,*) DO jnj = 1, jpisl WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl ) END DO ENDIF ! 2. Invertion of the island matrix ! --------------------------------- ! 2.1 Call of an imsl routine for the matrix invertion CALL linrg( jpisl, aisl, jpisl, aislm1, jpisl ) ! 2.2 Control print IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' islmat: inverse of the island matrix' WRITE(numout,*) ' ~~~~~~' WRITE(numout,*) DO jnj = 1, jpisl WRITE(numout, '(12e11.3)') ( aislm1(jni,jnj), jni=1, jpisl ) END DO ENDIF ! 3. Output of aisl and aislm1 in numisl ! -------------------------------------- IF(lwp) WRITE(numisl,REC=jpisl+2) aisl, aislm1 CLOSE( numisl ) END SUBROUTINE islmat SUBROUTINE islbsf !!---------------------------------------------------------------------- !! *** ROUTINE islbsf *** !! !! ** Purpose : Compute the barotropic stream function associated with !! each island using FETI or preconditioned conjugate gradient !! method or read them in numisl. !! !! ** Method : Direct acces file !! !! ** file : numisl: barotropic stream function associated !! with each island of the domain !! !! ** Action : - bsfisl, the streamfunction which takes the value 1 !! over island ni, and 0 over the others islands !! - update 'numisl' file !! !! History : !! ! 87-10 (G. Madec) Original code !! ! 91-11 (G. Madec) !! ! 93-03 (G. Madec) release 7.1 !! ! 93-04 (M. Guyon) loops and suppress pointers !! ! 96-11 (A. Weaver) correction to preconditioning !! ! 98-02 (M. Guyon) FETI method !! ! 99-11 (M. Imbard) NetCDF FORMAT with IOIPSL !! 8.5 ! 02-08 (G. Madec) Free form, F90 !!---------------------------------------------------------------------- !! * Local declarations INTEGER :: ji, jj, jni, jii, jnp INTEGER :: iimlu, ijmlu, inmlu, ios, iju, iloc INTEGER :: ii, ij, icile, icut, inmax, indic REAL(wp) :: zepsr, zeplu, zgwgt REAL(wp) :: zep(jpisl) !!---------------------------------------------------------------------- ! 0. Initializations ! ================== ! set to zero the convergence indicator icile = 0 ! set to zero of bsfisl bsfisl(:,:,1:jpisl) = 0.e0 ! I. Lecture of bsfisl in numisl (if it exists) ! ============================================= ! open islands file ibloc = 4096 ilglo = ibloc*((jpiglo*jpjglo*jpbyt-1 )/ibloc+1) clname = 'islands' CALL ctlopn( numisl, 'islands', 'UNKNOWN', 'UNFORMATTED', 'DIRECT', & ilglo,numout,lwp,0) icut = 0 iimlu = 0 ijmlu = 0 inmlu = 0 zeplu = 0. READ (numisl, REC=1, IOSTAT=ios) iimlu, ijmlu, inmlu, zeplu IF( ios /= 0 ) GOTO 110 DO jni = 1, jpisl CALL read2( numisl, bsfisl(1,1,jni), 1, jni+1 ) END DO 110 CONTINUE ! the read precision is not the required one : icut=888 IF( zeplu > epsisl ) icut = 888 ! the read domain does not correspond to the model one : icut=999 IF( iimlu /= jpiglo .OR. ijmlu /= jpjglo .OR. inmlu /= jpisl ) icut = 999 ! Control print IF( icut == 999 ) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' islbsf : lecture bsfisl in numisl failed' WRITE(numout,*) ' ~~~~~~' WRITE(numout,*) ' icut= ', icut WRITE(numout,*) ' imlu= ', iimlu, ' jmlu= ', ijmlu, & ' nilu= ', inmlu, ' epsisl lu= ', zeplu WRITE(numout,*) ' the bsfisl are computed from zero' ENDIF ELSEIF( icut == 888 ) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' islbsf : lecture bsfisl in numisl done' WRITE(numout,*) ' ~~~~~~' WRITE(numout,*) ' the required accuracy is not reached' WRITE(numout,*) ' epsisl lu= ', zeplu,' epsisl required ', epsisl WRITE(numout,*) ' the bsfisl are computed from the read values' ENDIF ELSE IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' islbsf : lecture bsfisl in numisl done' WRITE(numout,*) ' ~~~~~~' ENDIF RETURN ENDIF ! II. Compute the bsfisl (if icut=888 or 999) ! ====================== ! save nmax inmax = nmax ! set the number of iteration of island computation nmax = nmisl ! Loop over islands ! ----------------- DO jni = 1, jpisl ! 1. Initalizations of island computation ! --------------------------------------- ! 1.0 Set the pcg solution gcb to zero gcb(:,:) = 0.e0 ! 1.1 Set first guess gcx either to zero or to the read bsfisl IF( icut == 999 ) THEN gcx(:,:) = 0.e0 ELSEIF( icut == 888 ) THEN IF(lwp) WRITE(numout,*) ' islbsf: bsfisl read are used as first guess' ! C a u t i o n : bsfisl masked because it contains 1 along island seaside gcx(:,:) = bsfisl(:,:,jni) * bmask(:,:) ENDIF ! 1.2 Right hand side of the stream FUNCTION equation # if defined key_mpp ! north fold treatment IF( npolj == 3 ) iloc = jpiglo -(nimpp-1+nimppt(nono+1)-1) IF( npolj == 4 ) iloc = jpiglo - 2*(nimpp-1) t2p1(:,1,1) = 0. ! north and south grid-points DO jii = 1, 2 DO jnp = 1, mnisl(jii,jni) ii = miisl(jnp,jii,jni) ij = mjisl(jnp,jii,jni) IF( ( npolj == 3 .OR. npolj == 4 ) .AND. ( ij == nlcj-1 .AND. jii == 1) ) THEN iju=iloc-ii+1 t2p1(iju,1,1) = t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) ELSE gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) ENDIF END DO END DO ! east and west grid-points DO jii = 3, 4 DO jnp = 1, mnisl(jii,jni) ii = miisl(jnp,jii,jni) ij = mjisl(jnp,jii,jni) gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) END DO END DO CALL mpplnks( gcb ) # else ! north and south grid-points DO jii = 1, 2 DO jnp = 1, mnisl(jii,jni) ii = miisl(jnp,jii,jni) ij = mjisl(jnp,jii,jni) IF( ( nperio == 3 .OR. nperio == 4 ) .AND. ( ij == jpj-1 .AND. jii == 1) ) THEN gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) ELSE gcb(ii,ij-jii+1) = gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) ENDIF END DO END DO ! east and west grid-points DO jii = 3, 4 DO jnp = 1, mnisl(jii,jni) ii = miisl(jnp,jii,jni) ij = mjisl(jnp,jii,jni) IF( bmask(ii-jii+3,ij) /= 0. ) THEN gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) ELSE ! east-west cyclic boundary conditions IF( ii-jii+3 == 1 ) THEN gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij) ENDIF ENDIF END DO END DO # endif ! 1.4 Preconditioned right hand side and absolute precision IF( nsolv == 3 ) THEN ! FETI method ncut = 0 rnorme=0. DO jj = 1, jpj DO ji = 1, jpi gcb (ji,jj) = bmask(ji,jj) * gcb(ji,jj) rnorme = rnorme + gcb(ji,jj) * gcb(ji,jj) END DO END DO CALL mpp_sum( rnorme ) IF(lwp) WRITE(numout,*) 'rnorme ', rnorme epsr = epsisl * epsisl * rnorme indic = 0 ELSE ncut = 0 rnorme = 0. DO jj = 1, jpj DO ji = 1, jpi gcb (ji,jj) = gcdprc(ji,jj) * gcb(ji,jj) zgwgt = gcdmat(ji,jj) * gcb(ji,jj) rnorme = rnorme + gcb(ji,jj) * zgwgt END DO END DO #if defined key_mpp CALL mpp_sum( rnorme ) #endif IF(lwp) WRITE(numout,*) 'rnorme ', rnorme epsr = epsisl * epsisl * rnorme indic = 0 ENDIF ! 3. PCG solver for gcp.gcx=gcb in monotask ! ------------------------------------------- IF(nsolv == 3) THEN ! FETI method ! precision to compute Islands matrix A epsilo = epsisl CALL solfet( indic ) ! ---> ! precision to compute grad PS epsilo = eps ELSE CALL solpcg( indic ) ENDIF ! 4. Save the solution in bsfisl ! ------------------------------ DO jj = 1, jpj DO ji = 1, jpi bsfisl(ji,jj,jni) = gcx(ji,jj) END DO END DO ! 5. Boundary conditions ! ---------------------- ! set to 1. coastal gridpoints of the island DO jnp = 1, mnisl(0,jni) ii = miisl(jnp,0,jni) ij = mjisl(jnp,0,jni) bsfisl(ii,ij,jni) = 1. END DO ! cyclic boundary conditions IF( nperio == 1 .OR. nperio == 4 ) THEN DO jj = 1, jpj bsfisl( 1 ,jj,jni) = bsfisl(jpim1,jj,jni) bsfisl(jpi,jj,jni) = bsfisl( 2 ,jj,jni) END DO ENDIF IF( nperio == 3 .OR. nperio == 4 ) THEN DO ji = 1, jpim1 iju = jpi-ji+1 bsfisl(ji,jpj-1,jni) = bsfisl(iju,jpj-2,jni) bsfisl(ji, jpj ,jni) = bsfisl(iju,jpj-3,jni) END DO ENDIF #if defined key_mpp CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. ) #endif ! 6. Control print ! ---------------- IF(lwp)WRITE(numout,*) IF(lwp)WRITE(numout,*) ' islbsf: island number: ', jni IF(lwp)WRITE (numout,9290) niter, res, SQRT(epsr)/epsisl zep(jni) = MAX(epsisl, res/(SQRT(epsr)/epsisl)) IF(indic < 0) THEN icile = icile-1 IF(lwp)WRITE(numout,*) ' pcg do not converge for island: ', jni IF(lwp)WRITE(numout,*) ' Precision reached: ',zep(jni) ENDIF 9290 FORMAT(' niter :',i4,' , res :',e20.10,' , gcb :',e20.10) ! End loop for islands ! -------------------- END DO ! 7. Reset PCG ! ------------ ! reset the number of iteration for pcg nmax = inmax ! reset to zero pcg arrays gcx (:,:) = 0.e0 gcxb (:,:) = 0.e0 gcb (:,:) = 0.e0 gcr (:,:) = 0.e0 gcdes(:,:) = 0.e0 gccd (:,:) = 0.e0 ! III. Output of bsfisl in numisl ! =============================== IF( icile == 0 .AND. icut /= 0 ) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*)' islbsf : write bsfisl in numisl ', numisl WRITE(numout,*)' ~~~~~~' ENDIF IF(lwp) WRITE(numisl,REC=1) jpiglo, jpjglo, jpisl, epsisl DO jni = 1, jpisl CALL write2( numisl, bsfisl(1,1,jni), 1, jni+1 ) END DO ENDIF IF( icile < 0 ) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' islbsf : number of island without convergence: ',ABS(icile) WRITE(numout,*) ' ~~~~~~' ENDIF zepsr = epsisl DO jni = 1, jpisl IF(lwp)WRITE(numout,*) ' isl ',jni,' precision reached ', zep(jni) zepsr = MAX( zep(jni), zepsr ) END DO IF( zepsr == 0. ) zepsr = epsisl IF(lwp) THEN WRITE(numout,*) ' save value of precision reached: ',zepsr WRITE(numout,*) WRITE(numout,*)' islbsf : save bsfisl in numisl ',numisl WRITE(numout,*)' ~~~~~~' ENDIF IF(lwp) WRITE(numisl,REC=1) jpiglo, jpjglo, jpisl, zepsr DO jni = 1, jpisl CALL write2( numisl, bsfisl(1,1,jni), 1, jni+1 ) END DO nstop = nstop + 1 ENDIF END SUBROUTINE islbsf