C $Header: /u/gcmpack/MITgcm/pkg/exch2/w2_set_f2f_index.F,v 1.2 2010/04/23 20:21:06 jmc Exp $ C $Name: $ #include "CPP_EEOPTIONS.h" #include "W2_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP 0 C !ROUTINE: W2_SET_F2F_INDEX C !INTERFACE: SUBROUTINE W2_SET_F2F_INDEX( myThid ) C !DESCRIPTION: C Set-up index correspondence matrix for connected Facet-Edges C !USES: IMPLICIT NONE C Tile topology settings data structures #include "SIZE.h" #include "EEPARAMS.h" #include "W2_EXCH2_SIZE.h" #include "W2_EXCH2_PARAMS.h" #include "W2_EXCH2_TOPOLOGY.h" C !INPUT PARAMETERS: C myThid :: my Thread Id number C (Note: not relevant since threading has not yet started) INTEGER myThid C !LOCAL VARIABLES: C === Local variables === C msgBuf :: Informational/error message buffer CHARACTER*(MAX_LEN_MBUF) msgBuf CHARACTER*1 edge(4) INTEGER i, j, ii, jj, i1, j1, k, lo, ll INTEGER errCnt INTEGER chk1, chk2, chk3, chk4, chk5, chk6 LOGICAL prtFlag CEOP DATA edge / 'N' , 'S' , 'E' , 'W' / WRITE(msgBuf,'(2A)') 'W2_SET_F2F_INDEX:', & ' index matrix for connected Facet-Edges:' CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid ) prtFlag = ABS(W2_printMsg).GE.2 & .OR. ( W2_printMsg .NE.0 .AND. myProcId.EQ.0 ) C-- Check Topology errCnt = 0 DO j=1,nFacets DO i=1,4 C- connected to: jj = INT(facet_link(i,j)) ii = MOD( NINT(facet_link(i,j)*10.), 10 ) IF ( facet_link(i,j).EQ. 0. ) THEN WRITE(msgBuf,'(3A,I3,A,F6.2,A)') & '** WARNING ** ', edge(i), '.Edge of facet #', & j, ' disconnected (facet_link=',facet_link(i,j),')' CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid ) CALL PRINT_MESSAGE( msgBuf,errorMessageUnit,SQUEEZE_RIGHT,1 ) ELSEIF ( jj.LT.1 .OR. jj.GT.nFacets & .OR. ii.LT.1 .OR. ii.GT.4 ) THEN WRITE(msgBuf,'(2A,I3,A,F6.2,A)') edge(i), '.Edge of facet #', & j, ' : bad connection (facet_link=',facet_link(i,j),')' CALL PRINT_ERROR( msgBuf, myThid ) errCnt = errCnt + 1 ELSE C- check connection is set-up both ways: j1 = INT(facet_link(ii,jj)) i1 = MOD( NINT(facet_link(ii,jj)*10.), 10 ) IF ( j1.NE.j .OR. i1.NE.i ) THEN WRITE(msgBuf,'(2(2A,I3,A),F5.2,A)') & edge(i), '.Edge facet #', j,' connect to: ', & edge(ii),'.Edge facet #',jj,' (',facet_link(i,j),' )' CALL PRINT_ERROR( msgBuf, myThid ) IF ( i1.GE.1 .AND. i1.LE.4 ) THEN WRITE(msgBuf,'(A,2(2A,I3,A),F5.2,A)') ' but ', & edge(ii),'.Edge facet #',jj,' connect to: ', & edge(i1),'.Edge facet #',j1,' (',facet_link(ii,jj),' )' ELSE WRITE(msgBuf,'(A,2(2A,I3,A),F5.2,A)') ' but ', & edge(ii),'.Edge facet #',jj,' connect to: ', & '?' ,'.Edge facet #',j1,' (',facet_link(ii,jj),' )' ENDIF CALL PRINT_ERROR( msgBuf, myThid ) errCnt = errCnt + 1 ENDIF ENDIF ENDDO ENDDO IF ( errCnt.GT.0 ) THEN WRITE(msgBuf,'(A,I3,A)') & ' W2_SET_F2F_INDEX: found', errCnt, ' Topology errors' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R W2_SET_F2F_INDEX' ENDIF C-- Check length of connection (facet size) between connected facet-edges errCnt = 0 DO j=1,nFacets DO i=1,4 C- Length of N or S Edge = x-size, E or W Edge = y-size lo = 2*(j-1) + (i+1)/2 lo = facet_dims( lo ) C- connected to: jj = INT(facet_link(i,j)) ii = MOD( NINT(facet_link(i,j)*10.), 10 ) IF ( jj.GE.1 ) THEN ll = 2*(jj-1)+(ii+1)/2 ll = facet_dims( ll ) IF ( lo.NE.ll ) THEN WRITE(msgBuf,'(3A,I3,A,I8,A)') 'Length of connection: ', & edge(i), '.Edge facet #', j, ' (=',lo,')' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(3A,I3,A,I8,A)') ' to: ', & edge(ii),'.Edge facet #', jj, ' (=',ll,') are different' CALL PRINT_ERROR( msgBuf, myThid ) errCnt = errCnt + 1 ENDIF C-- Set indices correspondence matrix (facet_pij): C suffix "so" for indices of source facet j ; C suffix "tg" for indices of target facet jj= INT(facet_link(i,j)) C pij(:,i,j) : matrix which gives so indices when applied to tg indices C iso = pij(1)*itg + pij(2)*jtg + oi C jso = pij(3)*itg + pij(4)*jtg + oj C- Default: Same orientation (works for N <-> S or E <-> W) facet_pij(1,i,j) = 1 facet_pij(2,i,j) = 0 facet_pij(3,i,j) = 0 facet_pij(4,i,j) = 1 C-- 1rst: cases with same orientation IF ( i.EQ.1 .AND. ii.EQ.2 ) THEN C- N <-- S ; [:,jHi]_so <-- [:,0]_tg facet_oi(i,j) = 0 facet_oj(i,j) = +facet_dims(2*j) ELSEIF ( i.EQ.2 .AND. ii.EQ.1 ) THEN C- S <-- N ; [:, 1 ]_so <-- [:,1+jHi]_tg facet_oi(i,j) = 0 facet_oj(i,j) = -facet_dims(2*jj) ELSEIF ( i.EQ.3 .AND. ii.EQ.4 ) THEN C- E <-- W ; [iHi,:]_so <-- [0,:]_tg facet_oi(i,j) = +facet_dims(2*j-1) facet_oj(i,j) = 0 ELSEIF ( i.EQ.4 .AND. ii.EQ.3 ) THEN C- W <-- E (i=4 & ii=3); [ 1 ,:]_so <-- [iHi,:]_tg facet_oi(i,j) = -facet_dims(2*jj-1) facet_oj(i,j) = 0 C-- 2nd : cases where orientation changes ELSEIF ( i.EQ.1 .AND. ii.EQ.4 ) THEN C- N <-- W ; mapping W.face target indices to N.face source indices C- [:,jHi]_so <-- [0,:]_tg facet_pij(1,i,j) = 0 facet_pij(2,i,j) =-1 facet_pij(3,i,j) = 1 facet_pij(4,i,j) = 0 facet_oi(i,j) = lo+1 facet_oj(i,j) = +facet_dims(2*j) ELSEIF ( i.EQ.2 .AND. ii.EQ.3 ) THEN C- S <-- E ; mapping E.face target indices to S.face source indices C- [:,1]_so <-- [1+iHi,:]_tg facet_pij(1,i,j) = 0 facet_pij(2,i,j) =-1 facet_pij(3,i,j) = 1 facet_pij(4,i,j) = 0 facet_oi(i,j) = lo+1 facet_oj(i,j) = -facet_dims(2*jj-1) ELSEIF ( i.EQ.3 .AND. ii.EQ.2 ) THEN C- E <-- S ; mapping S.face target indices to E.face source indices C- [iHi,:]_so <-- [:,0]_tg facet_pij(1,i,j) = 0 facet_pij(2,i,j) = 1 facet_pij(3,i,j) =-1 facet_pij(4,i,j) = 0 facet_oi(i,j) = +facet_dims(2*j-1) facet_oj(i,j) = lo+1 ELSEIF ( i.EQ.4 .AND. ii.EQ.1 ) THEN C- W <-- N ; mapping N.face target indices to W.face source indices C- [ 1 ,:]_so <-- [:,1+jHi]_tg facet_pij(1,i,j) = 0 facet_pij(2,i,j) = 1 facet_pij(3,i,j) =-1 facet_pij(4,i,j) = 0 facet_oi(i,j) = -facet_dims(2*jj) facet_oj(i,j) = lo+1 ELSE WRITE(msgBuf,'(2(3A,I3),A)') ' connect ', & edge(i), '.Edge (facet#', j, ' ) to: ', & edge(ii),'.Edge (facet#', jj,' )' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A)') & ' W2_SET_F2F_INDEX: Connection not supported' CALL PRINT_ERROR( msgBuf, myThid ) errCnt = errCnt + 1 ENDIF C-- Print resulting index matrix: IF ( prtFlag ) WRITE(W2_oUnit,'(2(3A,I3),A,4I3,A,2I6)') & ' ', edge(i), '.Edge Facet', j, ' <-- ', & edge(ii),'.Edge Facet', jj, & ' : pij=', (facet_pij(k,i,j),k=1,4), & ' ; oi,oj=', facet_oi(i,j), facet_oj(i,j) ENDIF ENDDO ENDDO IF ( errCnt.GT.0 ) THEN WRITE(msgBuf,'(A,I3,A)') & ' W2_SET_F2F_INDEX: found', errCnt, ' Connection errors' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R W2_SET_F2F_INDEX' ENDIF C-- Check indices correspondence matrix reciprocity : C This is not necessary (if code above is right) ; However: C a) reciprocity is used later-on => good to check; C b) easy to add an error in setting offset => worth to check. errCnt = 0 DO j=1,nFacets DO i=1,4 C- connected to: jj = INT(facet_link(i,j)) ii = MOD( NINT(facet_link(i,j)*10.), 10 ) IF ( jj.GE.1 ) THEN C iso = pij(1)*( Pij(1)*iso + Pij(2)*jso + Oi ) C + pij(2)*( Pij(3)*iso + Pij(4)*jso + Oj ) C + oi C jso = pij(3)*( Pij(1)*iso + Pij(2)*jso + Oi ) C + pij(4)*( Pij(3)*iso + Pij(4)*jso + Oj ) C + oj C- Matrix product: chk1 = facet_pij(1,i,j)*facet_pij(1,ii,jj) & + facet_pij(2,i,j)*facet_pij(3,ii,jj) chk2 = facet_pij(1,i,j)*facet_pij(2,ii,jj) & + facet_pij(2,i,j)*facet_pij(4,ii,jj) chk3 = facet_pij(3,i,j)*facet_pij(1,ii,jj) & + facet_pij(4,i,j)*facet_pij(3,ii,jj) chk4 = facet_pij(3,i,j)*facet_pij(2,ii,jj) & + facet_pij(4,i,j)*facet_pij(4,ii,jj) C- Offsets: chk5 = facet_pij(1,i,j)*facet_oi(ii,jj) & + facet_pij(2,i,j)*facet_oj(ii,jj) & + facet_oi(i,j) chk6 = facet_pij(3,i,j)*facet_oi(ii,jj) & + facet_pij(4,i,j)*facet_oj(ii,jj) & + facet_oj(i,j) IF ( chk1.NE.1 .OR. chk2.NE.0 .OR. chk5.NE.0 .OR. & chk3.NE.0 .OR. chk4.NE.1 .OR. chk6.NE.0 ) THEN WRITE(msgBuf,'(2(3A,I3),A)') ' connect ', & edge(i), '.Edge (facet#', j, ' ) to: ', & edge(ii),'.Edge (facet#', jj,' )' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(A,4I4,2I8)') & ' Bug in Matrix Product:',chk1,chk2,chk3,chk4,chk5,chk6 CALL PRINT_ERROR( msgBuf, myThid ) errCnt = errCnt + 1 ENDIF ENDIF ENDDO ENDDO IF ( errCnt.GT.0 ) THEN WRITE(msgBuf,'(A,I3,A)') & ' W2_SET_F2F_INDEX: found', errCnt, ' bugs in Matrix product' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R W2_SET_F2F_INDEX' ENDIF RETURN END