1:   
    2:         MODULE types
    3:         
    4:   !    symbolic names for 4, 2, and 1 byte integers
    5:   	
    6:   	INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
    7:   	INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
    8:   	INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
    9:   
   10:   ! and for single and double precision reals
   11:   !	INTEGER, PARAMETER :: SP = KIND(1.0)
   12:   !	INTEGER, PARAMETER :: DP = KIND(1.0D0)
   13:   	INTEGER, PARAMETER :: SP = KIND(1.0D0)
   14:   
   15:         END MODULE types
   16:   
   17:   
   18:   	
   19:   module hpfprocessors
   20:   !HPF$ PROCESSORS PROC(4)
   21:   end module hpfprocessors
   22:   
   23:   
   24:   module diff_interface
   25:   
   26:   INTERFACE
   27:   subroutine diff_x(nx,ny,nz,h,f,df)
   28:   use types
   29:   use hpfprocessors
   30:   integer(I4B)                :: nx,ny,nz
   31:   real(SP) h
   32:   real(SP), dimension(:,:,:), intent(IN)  :: f
   33:   real(SP), dimension(:,:,:), intent(OUT) :: df
   34:   
   35:   !HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
   36:   !HPF$ ALIGN WITH *f                      :: df
   37:   end subroutine diff_x
   38:   
   39:   subroutine diff_y(nx,ny,nz,h,f,df)
   40:   use types
   41:   use hpfprocessors
   42:   integer(I4B)                :: nx,ny,nz
   43:   real(SP) h
   44:   real(SP), dimension(:,:,:), intent(IN)  :: f
   45:   real(SP), dimension(:,:,:), intent(OUT) :: df
   46:   
   47:   !HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
   48:   !HPF$ ALIGN WITH *f                      :: df
   49:   end subroutine diff_y
   50:   
   51:   subroutine diff_z(nx,ny,nz,h,f,df)
   52:   use types
   53:   use hpfprocessors
   54:   integer(I4B)                :: nx,ny,nz
   55:   real(SP) h
   56:   real(SP), dimension(:,:,:), intent(IN)  :: f
   57:   real(SP), dimension(:,:,:), intent(OUT) :: df
   58:   
   59:   !HPF$ DISTRIBUTE *(*,*,BLOCK) ONTO *PROC :: f
   60:   !HPF$ ALIGN WITH *f                      :: df
   61:   end subroutine diff_z
   62:   
   63:   END INTERFACE
   64:   end module diff_interface
   65:   
   66:   
   67:   
   68:   
   69:   
   70:   
   71:   
   72:   subroutine connect(g11,g12,g13,g22,g23,g33, &
   73:                      gup11,gup12,gup13,gup22,gup23,gup33, &
   74:                      c111,c112,c113,c122,c123,c133, & 
   75:                      c211,c212,c213,c222,c223,c233, & 
   76:                      c311,c312,c313,c322,c323,c333, &
   77:                      nx,ny,nz,h)
   78:   
   79:   !-----------------------------------------------------------------------
   80:   ! This routine returns covariant form of ricci tensor. Raise index before
   81:   ! using in evolution code.
   82:   !-----------------------------------------------------------------------
   83:   
   84:   use types
   85:   use hpfprocessors
   86:   
   87:   real(SP), dimension(:,:,:), intent(IN) :: g11,g12,g13,g22,g23,g33
   88:   real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
   89:   real(SP), dimension(:,:,:), intent(OUT):: c111,c112,c113,c122,c123,c133
   90:   real(SP), dimension(:,:,:), intent(OUT):: c211,c212,c213,c222,c223,c233
   91:   real(SP), dimension(:,:,:), intent(OUT):: c311,c312,c313,c322,c323,c333
   92:   
   93:   !HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
   94:   !HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
   95:   !HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
   96:   !HPF$ ALIGN WITH *g11 :: c111,c112,c113,c122,c123,c133
   97:   !HPF$ ALIGN WITH *g11 :: c211,c212,c213,c222,c223,c233
   98:   !HPF$ ALIGN WITH *g11 :: c311,c312,c313,c322,c323,c333
   99:   
  100:   integer(I4B)  nx, ny, nz
  101:   real(SP) h
  102:   
  103:   real(SP), dimension(size(g11,1),size(g11,2),size(g11,3)) :: & 
  104:        t11,t12,t13,t21,t22,t23,t31,t32,t33
  105:   
  106:   !HPF$ ALIGN WITH  g11 ::  t11,t12,t13,t21,t22,t23,t31,t32,t33
  107:   
  108:   integer(I4B) i,j,k
  109:   
  110:   !............................................................................
  111:   INTERFACE
  112:   subroutine connect1(g11,g12,g13,g22,g23,g33, &
  113:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  114:                      c111,c112,c113,c122,c123,c133, & 
  115:                      c211,c212,c213,c222,c223,c233, & 
  116:                      c311,c312,c313,c322,c323,c333, &
  117:                      nx,ny,nz,h)
  118:   
  119:   use types
  120:   use hpfprocessors
  121:   use diff_interface
  122:   
  123:   real(SP), dimension(:,:,:), intent(IN) :: g11,g12,g13,g22,g23,g33
  124:   real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
  125:   real(SP), dimension(:,:,:), intent(INOUT):: c111,c112,c113,c122,c123,c133
  126:   real(SP), dimension(:,:,:), intent(INOUT):: c211,c212,c213,c222,c223,c233
  127:   real(SP), dimension(:,:,:), intent(INOUT):: c311,c312,c313,c322,c323,c333
  128:   
  129:   !HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
  130:   !HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
  131:   !HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
  132:   !HPF$ ALIGN WITH *g11 :: c111,c112,c113,c122,c123,c133
  133:   !HPF$ ALIGN WITH *g11 :: c211,c212,c213,c222,c223,c233
  134:   !HPF$ ALIGN WITH *g11 :: c311,c312,c313,c322,c323,c333
  135:   
  136:   integer(I4B)  nx, ny, nz
  137:   real(SP) h
  138:   end subroutine connect1
  139:   
  140:   subroutine connect2(g11,g12,g13,g22,g23,g33, &
  141:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  142:                      c111,c112,c113,c122,c123,c133, & 
  143:                      c211,c212,c213,c222,c223,c233, & 
  144:                      c311,c312,c313,c322,c323,c333, &
  145:                      nx,ny,nz,h)
  146:   
  147:   use types
  148:   use hpfprocessors
  149:   use diff_interface
  150:   
  151:   real(SP), dimension(:,:,:), intent(IN) :: g11,g12,g13,g22,g23,g33
  152:   real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
  153:   real(SP), dimension(:,:,:), intent(INOUT):: c111,c112,c113,c122,c123,c133
  154:   real(SP), dimension(:,:,:), intent(INOUT):: c211,c212,c213,c222,c223,c233
  155:   real(SP), dimension(:,:,:), intent(INOUT):: c311,c312,c313,c322,c323,c333
  156:   
  157:   !HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
  158:   !HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
  159:   !HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
  160:   !HPF$ ALIGN WITH *g11 :: c111,c112,c113,c122,c123,c133
  161:   !HPF$ ALIGN WITH *g11 :: c211,c212,c213,c222,c223,c233
  162:   !HPF$ ALIGN WITH *g11 :: c311,c312,c313,c322,c323,c333
  163:   
  164:   integer(I4B)  nx, ny, nz
  165:   real(SP) h
  166:   end subroutine connect2
  167:   
  168:   subroutine connect3(g11,g12,g13,g22,g23,g33, &
  169:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  170:                      c111,c112,c113,c122,c123,c133, & 
  171:                      c211,c212,c213,c222,c223,c233, & 
  172:                      c311,c312,c313,c322,c323,c333, &
  173:                      nx,ny,nz,h)
  174:   
  175:   use types
  176:   use hpfprocessors
  177:   use diff_interface
  178:   
  179:   real(SP), dimension(:,:,:), intent(IN) :: g11,g12,g13,g22,g23,g33
  180:   real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
  181:   real(SP), dimension(:,:,:), intent(INOUT):: c111,c112,c113,c122,c123,c133
  182:   real(SP), dimension(:,:,:), intent(INOUT):: c211,c212,c213,c222,c223,c233
  183:   real(SP), dimension(:,:,:), intent(INOUT):: c311,c312,c313,c322,c323,c333
  184:   
  185:   !HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
  186:   !HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
  187:   !HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
  188:   !HPF$ ALIGN WITH *g11 :: c111,c112,c113,c122,c123,c133
  189:   !HPF$ ALIGN WITH *g11 :: c211,c212,c213,c222,c223,c233
  190:   !HPF$ ALIGN WITH *g11 :: c311,c312,c313,c322,c323,c333
  191:   
  192:   integer(I4B)  nx, ny, nz
  193:   real(SP) h
  194:   end subroutine connect3
  195:   
  196:   subroutine connect4(g11,g12,g13,g22,g23,g33, &
  197:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  198:                      c111,c112,c113,c122,c123,c133, & 
  199:                      c211,c212,c213,c222,c223,c233, & 
  200:                      c311,c312,c313,c322,c323,c333, &
  201:                      nx,ny,nz,h)
  202:   
  203:   use types
  204:   use hpfprocessors
  205:   use diff_interface
  206:   
  207:   real(SP), dimension(:,:,:), intent(IN) :: g11,g12,g13,g22,g23,g33
  208:   real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
  209:   real(SP), dimension(:,:,:), intent(INOUT):: c111,c112,c113,c122,c123,c133
  210:   real(SP), dimension(:,:,:), intent(INOUT):: c211,c212,c213,c222,c223,c233
  211:   real(SP), dimension(:,:,:), intent(INOUT):: c311,c312,c313,c322,c323,c333
  212:   
  213:   !HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
  214:   !HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
  215:   !HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
  216:   !HPF$ ALIGN WITH *g11 :: c111,c112,c113,c122,c123,c133
  217:   !HPF$ ALIGN WITH *g11 :: c211,c212,c213,c222,c223,c233
  218:   !HPF$ ALIGN WITH *g11 :: c311,c312,c313,c322,c323,c333
  219:   
  220:   integer(I4B)  nx, ny, nz
  221:   real(SP) h
  222:   end subroutine connect4
  223:   
  224:   subroutine connect5(g11,g12,g13,g22,g23,g33, &
  225:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  226:                      c111,c112,c113,c122,c123,c133, & 
  227:                      c211,c212,c213,c222,c223,c233, & 
  228:                      c311,c312,c313,c322,c323,c333, &
  229:                      nx,ny,nz,h)
  230:   
  231:   use types
  232:   use hpfprocessors
  233:   use diff_interface
  234:   
  235:   real(SP), dimension(:,:,:), intent(IN) :: g11,g12,g13,g22,g23,g33
  236:   real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
  237:   real(SP), dimension(:,:,:), intent(INOUT):: c111,c112,c113,c122,c123,c133
  238:   real(SP), dimension(:,:,:), intent(INOUT):: c211,c212,c213,c222,c223,c233
  239:   real(SP), dimension(:,:,:), intent(INOUT):: c311,c312,c313,c322,c323,c333
  240:   
  241:   !HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
  242:   !HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
  243:   !HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
  244:   !HPF$ ALIGN WITH *g11 :: c111,c112,c113,c122,c123,c133
  245:   !HPF$ ALIGN WITH *g11 :: c211,c212,c213,c222,c223,c233
  246:   !HPF$ ALIGN WITH *g11 :: c311,c312,c313,c322,c323,c333
  247:   
  248:   integer(I4B)  nx, ny, nz
  249:   real(SP) h
  250:   end subroutine connect5
  251:   
  252:   subroutine connect6(g11,g12,g13,g22,g23,g33, &
  253:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  254:                      c111,c112,c113,c122,c123,c133, & 
  255:                      c211,c212,c213,c222,c223,c233, & 
  256:                      c311,c312,c313,c322,c323,c333, &
  257:                      nx,ny,nz,h)
  258:   
  259:   use types
  260:   use hpfprocessors
  261:   use diff_interface
  262:   
  263:   real(SP), dimension(:,:,:), intent(IN) :: g11,g12,g13,g22,g23,g33
  264:   real(SP), dimension(:,:,:), intent(IN) :: gup11,gup12,gup13,gup22,gup23,gup33
  265:   real(SP), dimension(:,:,:), intent(INOUT):: c111,c112,c113,c122,c123,c133
  266:   real(SP), dimension(:,:,:), intent(INOUT):: c211,c212,c213,c222,c223,c233
  267:   real(SP), dimension(:,:,:), intent(INOUT):: c311,c312,c313,c322,c323,c333
  268:   
  269:   !HPF$ DISTRIBUTE g11 *(*,*,BLOCK) ONTO *PROC
  270:   !HPF$ ALIGN WITH *g11 :: g12,g13,g22,g23,g33
  271:   !HPF$ ALIGN WITH *g11 :: gup11,gup12,gup13,gup22,gup23,gup33
  272:   !HPF$ ALIGN WITH *g11 :: c111,c112,c113,c122,c123,c133
  273:   !HPF$ ALIGN WITH *g11 :: c211,c212,c213,c222,c223,c233
  274:   !HPF$ ALIGN WITH *g11 :: c311,c312,c313,c322,c323,c333
  275:   
  276:   integer(I4B)  nx, ny, nz
  277:   real(SP) h
  278:   end subroutine connect6
  279:   
  280:   END INTERFACE
  281:   !............................................................................
  282:   
  283:   ! calculate  Christoffell tensor piecewise.
  284:   
  285:         c111=0
  286:         c112=0
  287:         c113=0
  288:         c122=0
  289:         c123=0
  290:         c133=0
  291:   
  292:         c211=0
  293:         c212=0
  294:         c213=0
  295:         c222=0
  296:         c223=0
  297:         c233=0
  298:   
  299:         c311=0
  300:         c312=0
  301:         c313=0
  302:         c322=0
  303:         c323=0
  304:         c333=0
  305:   
  306:         call connect1(g11,g12,g13,g22,g23,g33, &
  307:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  308:                      c111,c112,c113,c122,c123,c133, & 
  309:                      c211,c212,c213,c222,c223,c233, & 
  310:                      c311,c312,c313,c322,c323,c333, &
  311:                      nx,ny,nz,h)
  312:   
  313:         call connect2(g11,g12,g13,g22,g23,g33, &
  314:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  315:                      c111,c112,c113,c122,c123,c133, & 
  316:                      c211,c212,c213,c222,c223,c233, & 
  317:                      c311,c312,c313,c322,c323,c333, &
  318:                      nx,ny,nz,h)
  319:   
  320:         call connect3(g11,g12,g13,g22,g23,g33, &
  321:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  322:                      c111,c112,c113,c122,c123,c133, & 
  323:                      c211,c212,c213,c222,c223,c233, & 
  324:                      c311,c312,c313,c322,c323,c333, &
  325:                      nx,ny,nz,h)
  326:   
  327:         call connect4(g11,g12,g13,g22,g23,g33, &
  328:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  329:                      c111,c112,c113,c122,c123,c133, & 
  330:                      c211,c212,c213,c222,c223,c233, & 
  331:                      c311,c312,c313,c322,c323,c333, &
  332:                      nx,ny,nz,h)
  333:   
  334:         call connect5(g11,g12,g13,g22,g23,g33, &
  335:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  336:                      c111,c112,c113,c122,c123,c133, & 
  337:                      c211,c212,c213,c222,c223,c233, & 
  338:                      c311,c312,c313,c322,c323,c333, &
  339:                      nx,ny,nz,h)
  340:   
  341:         call connect6(g11,g12,g13,g22,g23,g33, &
  342:                      gup11,gup12,gup13,gup22,gup23,gup33, &
  343:                      c111,c112,c113,c122,c123,c133, & 
  344:                      c211,c212,c213,c222,c223,c233, & 
  345:                      c311,c312,c313,c322,c323,c333, &
  346:                      nx,ny,nz,h)
  347:   
  348:   return 
  349:   end subroutine connect
  350:   
  351:   
  352:   
  353:   
  354:   
  355:   
  356:   
  357: