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