text: ##########################################################################
  #
  # MPI+Vectorized seq code generator
  #
  # This routine is called by op2_fortran which parses the input files
  #
  # It produces a file xxx_veckernel.F90 for each kernel,
  # plus a master kernel file
  #
  ##########################################################################
  
  import re
  import datetime
  import os
  import glob
  import util
  
  def comm(line):
    global file_text, FORTRAN
    global depth
    if len(line) == 0:
      prefix = ''
    else:
      prefix = ' '*depth
    if len(line) == 0:
      file_text +='\n'
    elif FORTRAN:
      file_text +=prefix+'! '+line+'\n'
    elif CPP:
      file_text +=prefix+'//'+line+'\n'
  
  def rep(line,m):
    global dims, idxs, typs, indtyps, inddims
  
    if FORTRAN:
      if m < len(inddims):
        line = re.sub('INDDIM',str(inddims[m]),line)
        line = re.sub('INDTYP',str(indtyps[m]),line)
  
      line = re.sub('INDARG','ind_arg'+str(m+1),line)
      line = re.sub('DIMS',str(dims[m]),line)
      line = re.sub('ARG','arg'+str(m+1),line)
      line = re.sub('TYP',typs[m],line)
      line = re.sub('IDX',str(int(idxs[m])),line)
    elif CPP:
      line = re.sub('INDDIM',str(inddims[m]),line)
      line = re.sub('INDTYP',str(indtyps[m]),line)
  
      line = re.sub('INDARG','ind_arg'+str(m),line)
      line = re.sub('DIM',str(dims[m]),line)
      line = re.sub('ARG','arg'+str(m),line)
      line = re.sub('TYP',typs[m],line)
      line = re.sub('IDX',str(int(idxs[m])),line)
    return line
  
  def code(text):
    global file_text, FORTRAN, CPP, g_m
    global depth
    if len(text) == 0:
      file_text += '\n'
      return
    if len(text) == 0:
      prefix = ''
    else:
      prefix = ' '*depth
    if FORTRAN:
      file_text += prefix+rep(text,g_m)+'\n'
    elif CPP:
      file_text += prefix+rep(text,g_m)+'\n'
  
  def code_pre(text):
    global file_text, FORTRAN, CPP, g_m
    if FORTRAN:
      file_text += rep(text,g_m)+'\n'
    elif CPP:
      file_text += rep(text,g_m)+'\n'
  
  def DO(i,start,finish):
    global file_text, FORTRAN, CPP, g_m
    global depth
    if FORTRAN:
      code('DO '+i+' = '+start+', '+finish+'-1, 1')
    elif CPP:
      code('for ( int '+i+'='+start+'; '+i+'<'+finish+'; '+i+'++ ){')
    depth += 2
  
  def DO2(i,start,finish,step):
    global file_text, FORTRAN, CPP, g_m
    global depth
    if FORTRAN:
      code('DO '+i+' = '+start+', '+finish+'-1, '+step)
    elif CPP:
      code('for ( int '+i+'='+start+'; '+i+'<'+finish+'; '+i+'++ ){')
    depth += 2
  
  def DO3(i,start,finish):
    global file_text, FORTRAN, CPP, g_m
    global depth
    if FORTRAN:
      code('DO '+i+' = '+start+', '+finish+', 1')
    elif CPP:
      code('for ( int '+i+'='+start+'; '+i+'<'+finish+'; '+i+'++ ){')
    depth += 2
  
  def FOR(i,start,finish):
    global file_text, FORTRAN, CPP, g_m
    global depth
    if FORTRAN:
      code('DO '+i+' = '+start+', '+finish+'-1')
    elif CPP:
      code('for ( int '+i+'='+start+'; '+i+'<'+finish+'; '+i+'++ ){')
    depth += 2
  
  def ENDDO():
    global file_text, FORTRAN, CPP, g_m
    global depth
    depth -= 2
    if FORTRAN:
      code('END DO')
    elif CPP:
      code('}')
  
  def ENDFOR():
    global file_text, FORTRAN, CPP, g_m
    global depth
    depth -= 2
    if FORTRAN:
      code('END DO')
    elif CPP:
      code('}')
  
  def IF(line):
    global file_text, FORTRAN, CPP, g_m
    global depth
    if FORTRAN:
      code('IF ('+line+') THEN')
    elif CPP:
      code('if ('+ line + ') {')
    depth += 2
  
  def ELSE():
    global file_text, FORTRAN, CPP, g_m
    global depth
    depth -= 2
    if FORTRAN:
      code('ELSE')
    elif CPP:
      code('else {')
    depth += 2
  
  def ENDIF():
    global file_text, FORTRAN, CPP, g_m
    global depth
    depth -= 2
    if FORTRAN:
      code('END IF')
    elif CPP:
      code('}')
  
  def para_parse(text, j, op_b, cl_b):
      """Parsing code block, i.e. text to find the correct closing brace"""
  
      depth = 0
      loc2 = j
  
      while 1:
        if text[loc2] == op_b:
              depth = depth + 1
  
        elif text[loc2] == cl_b:
              depth = depth - 1
              if depth == 0:
                  return loc2
        loc2 = loc2 + 1
  
  arg_parse=util.arg_parse
  
  def op2_gen_mpivec(master, date, consts, kernels, hydra, bookleaf):
  
    global dims, idxs, typs, indtyps, inddims
    global FORTRAN, CPP, g_m, file_text, depth
  
    OP_ID   = 1;  OP_GBL   = 2;  OP_MAP = 3;
  
    OP_READ = 1;  OP_WRITE = 2;  OP_RW  = 3;
    OP_INC  = 4;  OP_MAX   = 5;  OP_MIN = 6;
  
    accsstring = ['OP_READ','OP_WRITE','OP_RW','OP_INC','OP_MAX','OP_MIN' ]
    typestrigs = ['INTEGER','INT','REAL','DOUBLE','CHAR','FLOAT' ]
  
    any_soa = 0
    for nk in range (0,len(kernels)):
      any_soa = any_soa or sum(kernels[nk]['soaflags'])
  
  ##########################################################################
  #  create new kernel file
  ##########################################################################
  
    for nk in range (0,len(kernels)):
      name  = kernels[nk]['name']
      nargs = kernels[nk]['nargs']
      dims  = kernels[nk]['dims']
      maps  = kernels[nk]['maps']
      var   = kernels[nk]['var']
      typs  = kernels[nk]['typs']
      accs  = kernels[nk]['accs']
      idxs  = kernels[nk]['idxs']
      inds  = kernels[nk]['inds']
      soaflags = kernels[nk]['soaflags']
      optflags = kernels[nk]['optflags']
      ninds   = kernels[nk]['ninds']
      inddims = kernels[nk]['inddims']
      indaccs = kernels[nk]['indaccs']
      indtyps = kernels[nk]['indtyps']
      invinds = kernels[nk]['invinds']
      mapnames = kernels[nk]['mapnames']
      invmapinds = kernels[nk]['invmapinds']
      mapinds = kernels[nk]['mapinds']
      nmaps = 0
      if ninds > 0:
        nmaps = max(mapinds)+1
  
      optidxs = [0]*nargs
      indopts = [-1]*nargs
      nopts = 0
      for i in range(0,nargs):
        if optflags[i] == 1 and maps[i] == OP_ID:
          optidxs[i] = nopts
          nopts = nopts+1
        elif optflags[i] == 1 and maps[i] == OP_MAP:
          if i == invinds[inds[i]-1]: #i.e. I am the first occurence of this dat+map combination
            optidxs[i] = nopts
            indopts[inds[i]-1] = i
            nopts = nopts+1
          else:
            optidxs[i] = optidxs[invinds[inds[i]-1]]
  
      for g_m in range(0,nargs):
        dims[g_m] = dims[g_m].replace('NPDE','6')
        dims[g_m] = dims[g_m].replace('DNTQMU','3')
        dims[g_m] = dims[g_m].replace('DNFCROW','3')
        dims[g_m] = dims[g_m].replace('DMAXZONE','500')
  
      needDimList = []
      for g_m in range(0,nargs):
        try:
          dims[g_m] = str(eval(dims[g_m]))
        except:
          needDimList = needDimList + [g_m]
          if maps[g_m] == OP_MAP:
            dims[g_m] = 'opDat'+str(inds[g_m])+'Dim'
          else:
            dims[g_m] = 'opDat'+str(g_m+1)+'Dim'
  
  #      if (not dims[g_m].isdigit()) and not (dims[g_m] in ['NPDE','DNTQMU','DNFCROW']):
  #        needDimList = needDimList + [g_m]
  #
  # set three logicals
  #
      j = -1
      for i in range(0,nargs):
        if maps[i] == OP_MAP and accs[i] == OP_INC:
          j = i
      ind_inc = j >= 0
  
      j = -1
      for i in range(0,nargs):
        if maps[i] == OP_GBL and accs[i] <> OP_READ:
          j = i
      reduct = j >= 0
  
      j = -1
      for i in range(0,nargs):
        if maps[i] == OP_MAP :
          j = i
      indirect_kernel = j > -1
  
      if bookleaf:
        for i in range(0,nargs):
          if 'LI' in dims[i]:
            dims[i] = dims[i].replace('LI','100')
      FORTRAN = 1;
      CPP     = 0;
      g_m = 0;
      file_text = ''
      depth = 0
  
  ##########################################################################
  #  Generate Header
  ##########################################################################
  #    if hydra :
  #      if indirect_kernel:
  #        #if name <> 'VFLUX_EDGEF': #'ACCUMEDGES':
  #          #print "skipping indirect kernel :", name
  #          continue
  #      elif name <> 'GRAD_VOLAPF': #UPDATEK - problems with op_wirtes, SRCSA_NODE
  #        print "skipping unspecified kernel :", name
  #        continue
  
      if hydra:
        code('MODULE '+kernels[nk]['mod_file'][4:]+'_MODULE')
      else:
        code('MODULE '+name.upper()+'_MODULE')
      code('USE OP2_FORTRAN_DECLARATIONS')
      code('USE OP2_FORTRAN_RT_SUPPORT')
      code('USE ISO_C_BINDING')
      if hydra == 0:
        code('USE OP2_CONSTANTS')
      if bookleaf:
        code('USE kinds_mod,    ONLY: ink,rlk')
        code('USE parameters_mod,ONLY: LI')
      code('')
  
  ####################################################################################
  #  generate the user kernel function - creating versions for vectorisation as needed
  ####################################################################################
      code('')
      code('CONTAINS')
      code('')
      if 1:#hydra == 0:
  #
  # First original version
  #
        funcall_in_kernel = 0
        comm('user function')
        if bookleaf:
          modfile = kernels[nk]['mod_file']
          prefixes=['./','ale/','utils/','io/','eos/','hydro/','mods/']
          prefix_i=0
          while (prefix_i<7 and (not os.path.exists(prefixes[prefix_i]+modfile))):
            prefix_i=prefix_i+1
          fid = open(prefixes[prefix_i]+modfile, 'r')
          text = fid.read()
          i = re.search('SUBROUTINE '+name+'\\b',text).start() #text.find('SUBROUTINE '+name)
          j = i + 10 + text[i+10:].find('SUBROUTINE '+name) + 12 + len(name)
          kernel_text = text[i:j]
        elif hydra:
          file_text += '!DEC$ ATTRIBUTES FORCEINLINE :: ' + name + '\n'
          modfile = kernels[nk]['mod_file'][4:]
          filename = modfile.split('_')[1].lower() + '/' + modfile.split('_')[0].lower() + '/' + name + '.F95'
          if not os.path.isfile(filename):
            filename = modfile.split('_')[1].lower() + '/' + modfile.split('_')[0].lower() + '/' + name[:-1] + '.F95'
          fid = open(filename, 'r')
          text = fid.read()
          fid.close()
          p = re.compile('SUBROUTINE\\s+\\b'+name+'\\b',re.IGNORECASE)
          i = p.search(text).start()
          p = re.compile('END\\s+SUBROUTINE\\s+.*\\n',re.IGNORECASE)
          j = p.search(text).end()
          text= text[i:j]
          #get rid of #ifdef COMPEX
          i = text.find('#ifdef COMPLEX')
          while i>=0:
            j = i+text[i:].find('#else')
            k = i+text[i:].find('#endif')
            #if there is a nested macro, bail
            l = i+text[i:].find('#if')
            if l>i and l < k:
              print 'Nested macro in '+name+'prevents parsing for vectorization'
              exit(-1)
            #otherwise remove the body
            if j>i and j<k:
              j = j+text[j:].find('\n')+1
              k2 = k+text[k:].find('\n')+1
              text=text[:i]+text[j:k]+text[k2:]
            else:
              k2 = k+text[k:].find('\n')+1
              text=text[:i]+text[k2:]
            i = text.find('#ifdef COMPLEX')
            
          res=re.search('\\bcall\\b',text)
          if res <> None:
            funcall_in_kernel = 1
  
          kernel_text = text
        else:
          file_name = name+'.inc'
          f = open(file_name, 'r')
          kernel_text = f.read()
          f.close()
        file_text += kernel_text
        code('')
  
        #code('#define SIMD_VEC 8')
        code('#define SIMD_VEC 4')
  
        if not funcall_in_kernel:
          if 'VISCWALL' in name or 'IFLUX_EDGE' in name or 'BCSSETK' in name or 'BCSRESK' in name or 'JACOB_WALL_CORRECTIONS' in name:
            code('#define VDECTORIZE') #disable these
          else:
            code('#define VECTORIZE')
        else:
          print 'Function call in '+name+' prevents vectorisation - disabled for this kernel'
  #
  # Modified vectorisable version if its an indirect kernel
  # - direct kernels can be vectorised without modification
  #
        if indirect_kernel:
          code('#ifdef VECTORIZE')
          comm('user function -- modified for vectorisation')
  
          p = re.compile('SUBROUTINE\\s+\\b'+name+'\\b',re.IGNORECASE)
          i = p.search(kernel_text).start()
  
          if(i < 0):
            print "\n********"
            print "Error: cannot locate user kernel function name: "+name+" - Aborting code generation"
            exit(2)
          i2 = i
          j = kernel_text[i:].find('(')
          k = para_parse(kernel_text, i+j, '(', ')')
          l = kernel_text[k:].find('END'+'\\s+\\b'+'SUBROUTINE')
          paramlist = kernel_text[j+1:k]
          paramlist = paramlist.replace('&','')
          paramlist = paramlist.replace('\n','')
          paramlist = paramlist.replace(' ','')
          para = paramlist.split(',')
          #remove direct vars from para
          para_ind = []
          for i in range(0,nargs):
            if maps[i] == OP_ID or (maps[i]==OP_GBL and accs[i]==OP_READ) :
              para[i] = 'DIRECT'
            else:
              para[i] = para[i].strip()
  
          code('SUBROUTINE '+name+'_vec('+kernel_text[j+1:k]+',idx)')
          depth = depth + 2
          code('!dir$ attributes vector :: '+name+'_vec')
  
          #locate and remove non-indirection parameters and global parameters
          body_lines = kernel_text[k+1:l].split('\n')
          typestr_list = ['r\*8','r\*4','i\*4','real','integer','logical','integer\*4','real\*8','real\*4','logical\*1']
          current_type = ''
          typelist = []
          varlist = []
          dimlist = []
          typelist.append('INTEGER(KIND=4) ::')
          varlist.append('idx')
          dimlist.append('')
          last_bl_idx=0
          bl_idx = 0
          while bl_idx < len(body_lines):
            bl = body_lines[bl_idx].strip()
            if any((re.search(r'\b'+typestr+r'\b',bl.lower()) <> None and \
                    re.search(r'\b'+typestr+r'\b',bl.lower()).start()==0) for typestr in typestr_list):
              #concatenate lines
              while bl.find('&') <> -1: #TODO: commented out &?
                bl_idx = bl_idx + 1
                if body_lines[bl_idx].strip().find('#') == 0:
                  bl_idx = bl_idx + 1
                bl = bl[:bl.find('&')] + body_lines[bl_idx].strip()[1:]
  
              colons = bl.find('::') 
              if colons > -1:
                current_type = bl[:colons+2]
              else:
                idx=0
                while not bl[idx].isspace() and not bl[idx]=='(':
                  idx = idx+1
                if bl[idx] == '(':
                  idx = arg_parse(bl,idx)
                current_type = bl[:idx+1] 
            else:
              current_type = ''
  
            if current_type <> '':
              last_bl_idx=bl_idx
              idx = len(current_type)
              prev_idx = len(current_type)
              curr_dim = ''
              while idx < len(bl):
                if bl[idx].isspace():
                  idx = idx+1
                  if len(bl[prev_idx:idx].strip())==0:
                    prev_idx = idx
                elif bl[idx] == '(':
                  idx2 = arg_parse(bl,idx)
                  dimlist.append(bl[idx+1:idx2])
                  idx = idx2+1
                  typelist.append(current_type)
                  varlist.append(bl[prev_idx:idx])
                  prev_idx = idx
                elif bl[idx] == ',':
                  if len(bl[prev_idx:idx].strip()):
                    typelist.append(current_type)
                    varlist.append(bl[prev_idx:idx].strip())
                    dimlist.append(curr_dim)
                  idx = idx + 1
                  prev_idx = idx
                elif idx == len(bl)-1:
                  idx = idx + 1
                  typelist.append(current_type)
                  varlist.append(bl[prev_idx:idx].strip())
                  dimlist.append(curr_dim)
                else:
                  idx = idx+1
            bl_idx = bl_idx + 1
          
          depth = depth - 2
          for i in range(0,len(varlist)):
            for j in range(0,len(para)):
              if (not para[j] == 'DIRECT') and (not re.search(r'\b'+para[j]+r'\b', varlist[i]) == None):
                if dimlist[i] <> '':
                  dimsstr=dimlist[i]
                elif j in needDimList:
                  dimsstr = '*'
                else:
                  dimsstr = dims[j]
                if maps[j] == OP_MAP:
                  if (accs[j] == OP_INC or accs[j] == OP_RW or accs[j] == OP_WRITE):
                    typelist[i] = typs[j]+', DIMENSION(SIMD_VEC,'+dimsstr+') ::'
                    varlist[i] = para[j]
                  if (accs[j] == OP_READ):
                    typelist[i] = typs[j]+', DIMENSION(SIMD_VEC,'+dimsstr+'), INTENT(IN) ::'
                    varlist[i] = para[j]
                elif maps[j] == OP_GBL:
                  if ninds>0:
                    typelist[i] = typs[j]+', DIMENSION(SIMD_VEC,'+dimsstr+') ::'
                  else:
                    typelist[i] = typs[j]+', DIMENSION('+dimsstr+') ::'
                  varlist[i] = para[j]
  
          types_inserted = 0
          bl_idx = 0
          while bl_idx < len(body_lines):
            bl = body_lines[bl_idx].strip()
            if any((re.search(r'\b'+typestr+r'\b',bl.lower()) <> None and \
                    re.search(r'\b'+typestr+r'\b',bl.lower()).start()==0) for typestr in typestr_list):
              while bl.find('&') <> -1: #TODO: commented out &?
                bl_idx = bl_idx + 1
                if body_lines[bl_idx].strip().find('#') == 0:
                  bl_idx = bl_idx + 1
                bl = bl[:bl.find('&')] + body_lines[bl_idx].strip()[1:]
              if types_inserted == 0:
                for i in range(0,len(varlist)):
                  code('  '+typelist[i]+' '+varlist[i])
                types_inserted = 1
              bl_idx = bl_idx + 1
              continue
            temp = body_lines[bl_idx]
            for p in range(0,len(para)):
              temp = re.sub(r'(\b'+para[p]+r'\b\s*'+'\('+')', r'\1'+'idx,', temp)
              #TODO: only if dim is 1, otherwise vector op - fail?
              temp = re.sub(r'\b'+para[p]+r'\b(\s*\n|[^\(])', para[p]+'(idx,1)'+r'\1', temp+'\n')[:-1]
            kernel_line = temp
            if bl_idx == len(body_lines)-1:
              kernel_line = kernel_line.lower().replace(name.lower(),name+'_vec')
            if hydra and types_inserted==1:
              kernel_line = re.sub('\\bnpdes\\b','NPDE',kernel_line)
            code(kernel_line)
            bl_idx = bl_idx + 1
            
  
          code('#endif')
  
      ###  Hydra specific user kernel #### should not be in code generator .. to be fixed
      else:
  
        #
        # First original version
        #
  
        file_text += '!DEC$ ATTRIBUTES FORCEINLINE :: ' + name + '\n'
        modfile = kernels[nk]['mod_file'][4:]
        #print modfile
        modfile = modfile.replace('INIT_INIT','INIT')
        name2 = name.replace('INIT_INIT','INIT')
        #print modfile
        file_name = modfile.split('_')[1].lower() + '/' + modfile.split('_')[0].lower() + '/' + name2 + '.F95'
        if not os.path.isfile(file_name):
          file_name = modfile.split('_')[1].lower() + '/' + modfile.split('_')[0].lower() + '/' + name + '.F95'
        if not os.path.isfile(file_name):
          file_name = modfile.split('_')[1].lower() + '/' + modfile.split('_')[0].lower() + '/' + name2[:-1] + '.F95'
        fid = open(file_name, 'r')
        text = fid.read()
        fid.close()
        text = text.replace('recursive subroutine','subroutine')
        text = text.replace('module','!module')
        text = text.replace('contains','!contains')
        text = text.replace('end !module','!end module')
        file_text += text
        #code('#define SIMD_VEC 8')
        code('#define SIMD_VEC 4')
  
        #
        # Modified vectorisable version if its an indirect kernel
        # - direct kernels can be vectorised without modification
        #
  
        if indirect_kernel:
          code('#ifdef VECTORIZE')
          comm('user function -- modified for vectorisation')
          kernel_text = text
          p = re.compile('SUBROUTINE\\s+\\b'+name+'\\b',re.IGNORECASE)
          i = p.search(kernel_text).start()
  
          if(i < 0):
            print "\n********"
            print "Error: cannot locate user kernel function name: "+name+" - Aborting code generation"
            exit(2)
          i2 = i
          j = kernel_text[i:].find('(')
          k = para_parse(kernel_text, i+j, '(', ')')
          l = kernel_text[k:].find('END'+'\\s+\\b'+'SUBROUTINE')
          para = kernel_text[j+1:k].split(',')
          #remove direct vars from para
          para_ind = []
          for i in range(0,nargs):
            if maps[i] == OP_ID:
              para[i] = 'DIRECT'
  
          code('SUBROUTINE '+name+'_vec('+kernel_text[j+1:k]+',idx)')
          depth = depth + 2
          code('!dir$ attributes vector :: '+name+'_vec')
          code('IMPLICIT NONE')
          print needDimList
          for g_m in range(0,nargs):
            if g_m in needDimList:
              dimsstr = '*'
            else:
              dimsstr = '(DIMS)'
            if maps[g_m] == OP_MAP:
              if (accs[g_m] == OP_INC or accs[g_m] == OP_RW or accs[g_m] == OP_WRITE):
                code('TYP, DIMENSION(SIMD_VEC,'+dimsstr+') :: '+para[g_m])
              if (accs[g_m] == OP_READ):
                code('TYP, DIMENSION(SIMD_VEC,'+dimsstr+'), INTENT(IN) :: '+para[g_m])
            elif maps[g_m] == OP_GBL:
                code('TYP DIMENSION('+dimsstr+') :: '+para[g_m])
          code('INTEGER(4) :: idx')
  
  
      code('')
  
  ##########################################################################
  #  Generate wrapper to iterate over set
  ##########################################################################
  
      code('SUBROUTINE op_wrap_'+name+'( &')
      depth = depth + 2
      if nopts >0:
        code('&  optflags,        &')
      for g_m in range(0,ninds):
        if invinds[g_m] in needDimList:
          code('& opDat'+str(invinds[g_m]+1)+'Dim, &')
        code('& opDat'+str(invinds[g_m]+1)+'Local, &')
      for g_m in range(0,nargs):
        if maps[g_m] <> OP_MAP:
          if g_m in needDimList:
            code('& opDat'+str(g_m+1)+'Dim, &')
          code('& opDat'+str(g_m+1)+'Local, &')
      if nmaps > 0:
        k = []
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP and (not mapnames[g_m] in k):
            k = k + [mapnames[g_m]]
            code('& opDat'+str(invinds[inds[g_m]-1]+1)+'Map, &')
            code('& opDat'+str(invinds[inds[g_m]-1]+1)+'MapDim, &')
      code('& bottom,top)')
      code('implicit none')
      if nopts>0:
        code('INTEGER(kind=4), VALUE :: optflags')
      for g_m in range(0,ninds):
        if invinds[g_m] in needDimList:
          code('INTEGER(kind=4) opDat'+str(invinds[g_m]+1)+'Dim')
          code(typs[invinds[g_m]]+' opDat'+str(invinds[g_m]+1)+'Local(opDat'+str(invinds[g_m]+1)+'Dim,*)')
        else:
          code(typs[invinds[g_m]]+' opDat'+str(invinds[g_m]+1)+'Local('+str(dims[invinds[g_m]])+',*)')
      for g_m in range(0,nargs):
        if maps[g_m] <> OP_MAP:
          if g_m in needDimList:
            code('INTEGER(kind=4) opDat'+str(g_m+1)+'Dim')
        if maps[g_m] == OP_ID:
          if g_m in needDimList:
            code(typs[g_m]+' opDat'+str(g_m+1)+'Local(opDat'+str(g_m+1)+'Dim,*)')
          else:
            code(typs[g_m]+' opDat'+str(g_m+1)+'Local('+str(dims[g_m])+',*)')
        elif maps[g_m] == OP_GBL:
          if g_m in needDimList:
            code(typs[g_m]+' opDat'+str(g_m+1)+'Local(opDat'+str(g_m+1)+'Dim)')
          else:
            code(typs[g_m]+' opDat'+str(g_m+1)+'Local('+str(dims[g_m])+')')
      if nmaps > 0:
        k = []
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP and (not mapnames[g_m] in k):
            k = k + [mapnames[g_m]]
            code('INTEGER(kind=4) opDat'+str(invinds[inds[g_m]-1]+1)+'Map(*)')
            code('INTEGER(kind=4) opDat'+str(invinds[inds[g_m]-1]+1)+'MapDim')
  
      code('INTEGER(kind=4) bottom,top,i1, i2')
      if len(needDimList)>0:
        code('INTEGER(KIND=4) i3')
      if nmaps > 0:
        k = []
        line = 'INTEGER(kind=4) '
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP and (not mapinds[g_m] in k):
            k = k + [mapinds[g_m]]
            line += 'map'+str(mapinds[g_m]+1)+'idx, '
        code(line[:-2])
  
      #vars for globals - used when called with vectorisation
      for g_m in range(0,nargs):
        if maps[g_m] == OP_GBL and (accs[g_m] == OP_INC or accs[g_m] == OP_WRITE\
          or accs[g_m] == OP_MAX or accs[g_m] == OP_MIN):
          code('TYP dat'+str(g_m+1)+'(SIMD_VEC*'+dims[g_m]+')')
          code('!dir$ attributes align: 64:: dat'+str(g_m+1))
      code('')
  
  #
  # kernel call for indirect version
  #
      #If indirect kernel then add vector gather/scatter variables
      if indirect_kernel:
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP and (accs[g_m] == OP_READ \
            or accs[g_m] == OP_RW or accs[g_m] == OP_WRITE \
            or accs[g_m] == OP_INC):
            code('TYP dat'+str(g_m+1)+'(SIMD_VEC,DIMS)')
  
        code('')
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP and (accs[g_m] == OP_READ \
            or accs[g_m] == OP_RW or accs[g_m] == OP_WRITE \
            or accs[g_m] == OP_INC):
            code('!dir$ attributes align: 64:: dat'+str(g_m+1))
        code('')
  
      for g_m in range(0,ninds):
        code('!DIR$ ASSUME_ALIGNED opDat'+str(invinds[g_m]+1)+'Local : 64')
      for g_m in range(0,nargs):
        if maps[g_m] == OP_ID:
          code('!DIR$ ASSUME_ALIGNED opDat'+str(g_m+1)+'Local : 64')
      if nmaps > 0:
        k = []
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP and (not mapnames[g_m] in k):
            k = k + [mapnames[g_m]]
            code('!DIR$ ASSUME_ALIGNED opDat'+str(invinds[inds[g_m]-1]+1)+'Map : 64')
  
  
        code_pre('#ifdef VECTORIZE')
        DO2('i1','bottom','((top-1)/SIMD_VEC)*SIMD_VEC','SIMD_VEC')
        code('!DIR$ SIMD')
        DO3('i2','1','SIMD_VEC')
        k = []
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP :
            if (accs[g_m] == OP_READ or accs[g_m] == OP_RW or accs[g_m] == OP_WRITE) and (not mapinds[g_m] in k):
              k = k + [mapinds[g_m]]
              code('map'+str(mapinds[g_m]+1)+'idx = opDat'+str(invmapinds[inds[g_m]-1]+1)+'Map(1 + (i1+i2-1) * opDat'+str(invmapinds[inds[g_m]-1]+1)+'MapDim + '+str(int(idxs[g_m])-1)+') + 1')
  
        code('')
  
        #setup gathers
        code('')
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP :
            if (accs[g_m] == OP_READ or accs[g_m] == OP_RW):#and (not mapinds[g_m] in k):
              if optflags[g_m]==1:
                IF('BTEST(optflags,'+str(optidxs[g_m])+')')
              if g_m in needDimList:
                DO3('i3','1',dims[g_m])
                code('dat'+str(g_m+1)+'(i2,i3) = opDat'+\
                    str(invinds[inds[g_m]-1]+1)+'Local(i3,map'+str(mapinds[g_m]+1)+'idx)')
                ENDDO()
              else:
                for d in range(0,int(eval(dims[g_m]))):
                  code('dat'+str(g_m+1)+'(i2,'+str(d+1)+') = opDat'+\
                    str(invinds[inds[g_m]-1]+1)+'Local('+str(d+1)+',map'+str(mapinds[g_m]+1)+'idx)')
              if optflags[g_m]==1:
                ENDIF()
              code('')
            elif (accs[g_m] == OP_INC):
              code('dat'+str(g_m+1)+'(i2,:) = 0.0')
        ENDDO()
  
        #initialize globals
        for g_m in range(0,nargs):
          if maps[g_m] == OP_GBL:
            if accs[g_m] == OP_INC:
              code('dat'+str(g_m+1)+' = 0.0_8')
            elif accs[g_m] == OP_MAX:
              code('dat'+str(g_m+1)+' = -HUGE(dat'+str(g_m+1)+')')
            elif accs[g_m] == OP_MIN:
              code('dat'+str(g_m+1)+' = HUGE(dat'+str(g_m+1)+')')
  
        #vectorized kernel call
        code('!DIR$ SIMD')
        code('!DIR$ FORCEINLINE')
        DO3('i2','1','SIMD_VEC')
        comm('vectorized kernel call')
        line = 'CALL '+name+'_vec( &'
        indent = '\n'+' '*depth
        for g_m in range(0,nargs):
          if maps[g_m] == OP_ID:
            line = line + indent + '& opDat'+str(g_m+1)+'Local(1,(i1+i2-1)+1), &'
          if maps[g_m] == OP_MAP:
            line = line +indent + '& dat'+str(g_m+1)+', &'
          if maps[g_m] == OP_GBL:
            if accs[g_m] == OP_READ:
              line = line + indent +'& opDat'+str(g_m+1)+'Local(1), &'
            elif dims[g_m].isdigit() and eval(dims[g_m]) == 1:
              line = line + indent +'& dat'+str(g_m+1)+'(i2), &'
            else:
              line = line + indent +'& dat'+str(g_m+1)+'((DIMS)*(i2-1)+1:(DIMS)*(i2-1)+(DIMS)), &'
        line = line + indent +'& i2)'
        code(line)
        ENDDO()
  
        #do reductions
        for g_m in range(0,nargs):
          if maps[g_m] == OP_GBL:
            leftarg = 'opDat'+str(g_m+1)+'Local(1:(DIMS))'
            rightarg = 'dat'+str(g_m+1)+'((DIMS)*(i2-1)+1:(DIMS)*(i2-1)+(DIMS))'
            if accs[g_m] == OP_INC:
              DO3('i2','1','SIMD_VEC')
              code(leftarg+' = '+leftarg+' + '+rightarg)
              ENDDO()
            elif accs[g_m] == OP_MAX:
              DO3('i2','1','SIMD_VEC')
              code(leftarg+' = MAX('+leftarg+' , '+rightarg+')')
              ENDDO()
            elif accs[g_m] == OP_MIN:
              DO3('i2','1','SIMD_VEC')
              code(leftarg+' = MIN('+leftarg+' , '+rightarg+')')
              ENDDO()
        #do the scatters
        DO3('i2','1','SIMD_VEC')
        if nmaps > 0:
          for g_m in range(0,nargs):
            if maps[g_m] == OP_MAP :
              if (accs[g_m] == OP_INC or accs[g_m] == OP_RW or accs[g_m] == OP_WRITE):
                code('map'+str(mapinds[g_m]+1)+'idx = opDat'+str(invmapinds[inds[g_m]-1]+1)+'Map(1 + (i1+i2-1) * opDat'+str(invmapinds[inds[g_m]-1]+1)+'MapDim + '+str(int(idxs[g_m])-1)+') + 1')
          code('')
          for g_m in range(0,nargs):
            if maps[g_m] == OP_MAP :
              if optflags[g_m]==1:
                IF('BTEST(optflags,'+str(optidxs[g_m])+')')
              if (accs[g_m] == OP_INC ):
                if g_m in needDimList:
                  DO3('i3','1',dims[g_m])
                  code('opDat'+str(invinds[inds[g_m]-1]+1)+\
                      'Local(i3,map'+str(mapinds[g_m]+1)+'idx) = opDat'+str(invinds[inds[g_m]-1]+1)+\
                      'Local(i3,map'+str(mapinds[g_m]+1)+'idx) + dat'+str(g_m+1)+'(i2,i3)')
                  ENDDO()
                else:
                  for d in range(0,int(eval(dims[g_m]))):
                    code('opDat'+str(invinds[inds[g_m]-1]+1)+\
                      'Local('+str(d+1)+',map'+str(mapinds[g_m]+1)+'idx) = opDat'+str(invinds[inds[g_m]-1]+1)+\
                      'Local('+str(d+1)+',map'+str(mapinds[g_m]+1)+'idx) + dat'+str(g_m+1)+'(i2,'+str(d+1)+')')
                  code('')
              if (accs[g_m] == OP_WRITE or accs[g_m] == OP_RW):
                if g_m in needDimList:
                  DO3('i3','1',dims[g_m])
                  code('opDat'+str(invinds[inds[g_m]-1]+1)+\
                      'Local(i3,map'+str(mapinds[g_m]+1)+'idx) = dat'+str(g_m+1)+'(i2,i3)')
                  ENDDO()
                else:
                  for d in range(0,int(eval(dims[g_m]))):
                    code('opDat'+str(invinds[inds[g_m]-1]+1)+\
                      'Local('+str(d+1)+',map'+str(mapinds[g_m]+1)+'idx) = dat'+str(g_m+1)+'(i2,'+str(d+1)+')')
                  code('')
              if optflags[g_m]==1:
                ENDIF()
  
  
  
        ENDDO()
  
  
      #do reductions
      #TODO -- need exmple code
  
  #
  # kernel call for direct version
  #
      else:
        code_pre('#ifdef VECTORIZE')
        DO2('i1','bottom','((top-1)/SIMD_VEC)*SIMD_VEC','SIMD_VEC')
        #initialize globals
        for g_m in range(0,nargs):
          if maps[g_m] == OP_GBL:
            if accs[g_m] == OP_INC:
              code('dat'+str(g_m+1)+' = 0.0_8')
            elif accs[g_m] == OP_MAX:
              code('dat'+str(g_m+1)+' = -HUGE(dat'+str(g_m+1)+')')
            elif accs[g_m] == OP_MIN:
              code('dat'+str(g_m+1)+' = HUGE(dat'+str(g_m+1)+')')
        #vectorized kernel call
        code('!DIR$ SIMD')
        code('!DIR$ FORCEINLINE')
        DO3('i2','1','SIMD_VEC')
        comm('vectorized kernel call')
        line = '  CALL '+name+'( &'
        indent = '\n'+' '*depth
        for g_m in range(0,nargs):
          if maps[g_m] == OP_ID:
            line = line + indent + '& opDat'+str(g_m+1)+'Local(1,(i1+i2-1)+1)'
          if maps[g_m] == OP_GBL:
            if accs[g_m] == OP_READ:
              line = line + indent +'& opDat'+str(g_m+1)+'Local(1)'
            elif dims[g_m].isdigit() and eval(dims[g_m]) == 1:
              line = line + indent +'& dat'+str(g_m+1)+'(i2)'
            elif accs[g_m] <> OP_READ:
              line = line + indent +'& dat'+str(g_m+1)+'('+dims[g_m]+\
                '*(i2-1)+1:'+dims[g_m]+'*(i2-1)+'+dims[g_m]+')'
          if g_m < nargs-1:
            line = line +', &'
          else:
            line = line +' &'
        depth = depth - 2
        code(line + indent + '& )')
        depth = depth + 2
        ENDDO()
        code('')
  
        #do reductions
        for g_m in range(0,nargs):
          if maps[g_m] == OP_GBL:
            leftarg = 'opDat'+str(g_m+1)+'Local(1:(DIMS))'
            rightarg = 'dat'+str(g_m+1)+'((DIMS)*(i2-1)+1:(DIMS)*(i2-1)+(DIMS))'
            if accs[g_m] == OP_INC:
              DO3('i2','1','SIMD_VEC')
              code(leftarg+' = '+leftarg+' + '+rightarg)
              ENDDO()
            elif accs[g_m] == OP_MAX:
              DO3('i2','1','SIMD_VEC')
              code(leftarg+' = MAX('+leftarg+' , '+rightarg+')')
              ENDDO()
            elif accs[g_m] == OP_MIN:
              DO3('i2','1','SIMD_VEC')
              code(leftarg+' = MIN('+leftarg+' , '+rightarg+')')
              ENDDO()
  
  
      ENDDO()# end of SIMD_VEC length strided loop
  #
  # remainder of loop
  #
  
      comm('remainder')
      DO('i1','((top-1)/SIMD_VEC)*SIMD_VEC','top')
      depth = depth - 2
      code_pre('#else')
      code('!DIR$ FORCEINLINE')
      DO('i1','bottom','top')
      depth = depth - 2
      code_pre('#endif')
      depth = depth + 2
      k = []
      for g_m in range(0,nargs):
        if maps[g_m] == OP_MAP and (not mapinds[g_m] in k):
          k = k + [mapinds[g_m]]
          code('map'+str(mapinds[g_m]+1)+'idx = opDat'+str(invmapinds[inds[g_m]-1]+1)+'Map(1 + i1 * opDat'+str(invmapinds[inds[g_m]-1]+1)+'MapDim + '+str(int(idxs[g_m])-1)+')+1')
      comm('kernel call')
      line = '  CALL '+name+'( &'
      indent = '\n'+' '*depth
      for g_m in range(0,nargs):
        if maps[g_m] == OP_ID:
          line = line + indent + '& opDat'+str(g_m+1)+'Local(1,i1+1)'
        if maps[g_m] == OP_MAP:
          line = line +indent + '& opDat'+str(invinds[inds[g_m]-1]+1)+'Local(1,map'+str(mapinds[g_m]+1)+'idx)'
        if maps[g_m] == OP_GBL:
          line = line + indent +'& opDat'+str(g_m+1)+'Local(1)'
        if g_m < nargs-1:
          line = line +', &'
        else:
           line = line +' &'
      depth = depth - 2
      code(line + indent + '& )')
      depth = depth + 2
  
  
      ENDDO()
      depth = depth - 2
      code('END SUBROUTINE')
  
  ##########################################################################
  #  Generate SEQ host stub
  ##########################################################################
      code('SUBROUTINE '+name+'_host( userSubroutine, set, &'); depth = depth + 2
      for g_m in range(0,nargs):
        if g_m == nargs-1:
          code('& opArg'+str(g_m+1)+' )')
        else:
          code('& opArg'+str(g_m+1)+', &')
  
      code('')
      code('IMPLICIT NONE')
      code('character(kind=c_char,len=*), INTENT(IN) :: userSubroutine')
      code('type ( op_set ) , INTENT(IN) :: set')
      code('')
  
      for g_m in range(0,nargs):
        code('type ( op_arg ) , INTENT(IN) :: opArg'+str(g_m+1))
      code('')
  
      code('type ( op_arg ) , DIMENSION('+str(nargs)+') :: opArgArray')
      code('INTEGER(kind=4) :: numberOfOpDats')
      code('REAL(kind=4) :: dataTransfer')
      code('INTEGER(kind=4), DIMENSION(1:8) :: timeArrayStart')
      code('INTEGER(kind=4), DIMENSION(1:8) :: timeArrayEnd')
      code('REAL(kind=8) :: startTime')
      code('REAL(kind=8) :: endTime')
      code('INTEGER(kind=4) :: returnSetKernelTiming')
      code('INTEGER(kind=4) :: n_upper')
      code('type ( op_set_core ) , POINTER :: opSetCore')
      code('')
      for g_m in range(0,ninds):
        code('INTEGER(kind=4), POINTER, DIMENSION(:) :: opDat'+str(invinds[g_m]+1)+'Map')
        code('INTEGER(kind=4) :: opDat'+str(invinds[g_m]+1)+'MapDim')
        code(typs[invinds[g_m]]+', POINTER, DIMENSION(:) :: opDat'+str(invinds[g_m]+1)+'Local')
        code('INTEGER(kind=4) :: opDat'+str(invinds[g_m]+1)+'Cardinality')
        code('')
      for g_m in range(0,nargs):
        if maps[g_m] == OP_ID:
          code(typs[g_m]+', POINTER, DIMENSION(:) :: opDat'+str(g_m+1)+'Local')
          code('INTEGER(kind=4) :: opDat'+str(g_m+1)+'Cardinality')
          code('')
        if maps[g_m] == OP_GBL:
          code(typs[g_m]+', POINTER, DIMENSION(:) :: opDat'+str(g_m+1)+'Local')
  
      for g_m in range(0,nargs):
        if maps[g_m] == OP_MAP and optflags[g_m]==1:
          code(typs[g_m]+', POINTER, DIMENSION(:) :: opDat'+str(g_m+1)+'OptPtr')
  
  
      code('')
      code('INTEGER(kind=4) :: i1')
     
      if nopts>0:
        code('INTEGER(kind=4) :: optflags')
        code('optflags = 0')
        for i in range(0,nargs):
          if optflags[i] == 1:
            IF('opArg'+str(i+1)+'%opt == 1')
            code('optflags = IBSET(optflags,'+str(optidxs[i])+')')
            ENDIF()
      if nopts > 30:
        print 'ERROR: too many optional arguments to store flags in an integer'
  
  
      code('')
      code('numberOfOpDats = '+str(nargs))
      code('')
  
      for g_m in range(0,nargs):
        code('opArgArray('+str(g_m+1)+') = opArg'+str(g_m+1))
      code('')
  
      code('returnSetKernelTiming = setKernelTime('+str(nk)+' , userSubroutine//C_NULL_CHAR, &')
      code('& 0.d0, 0.00000,0.00000, 0)')
  
      code('call op_timers_core(startTime)')
      code('')
      #mpi halo exchange call
      code('n_upper = op_mpi_halo_exchanges(set%setCPtr,numberOfOpDats,opArgArray)')
  
      code('')
      code('opSetCore => set%setPtr')
      code('')
      for g_m in range(0,ninds):
        code('opDat'+str(invinds[g_m]+1)+'Cardinality = opArg'+str(invinds[g_m]+1)+'%dim * getSetSizeFromOpArg(opArg'+str(invinds[g_m]+1)+')')
        code('opDat'+str(invinds[g_m]+1)+'MapDim = getMapDimFromOpArg(opArg'+str(invinds[g_m]+1)+')')
      for g_m in range(0,nargs):
        if maps[g_m] == OP_ID:
          code('opDat'+str(g_m+1)+'Cardinality = opArg'+str(g_m+1)+'%dim * getSetSizeFromOpArg(opArg'+str(g_m+1)+')')
  
      for g_m in range(0,ninds):
        code('CALL c_f_pointer(opArg'+str(invinds[g_m]+1)+'%data,opDat'+str(invinds[g_m]+1)+'Local,(/opDat'+str(invinds[g_m]+1)+'Cardinality/))')
        code('CALL c_f_pointer(opArg'+str(invinds[g_m]+1)+'%map_data,opDat'+str(invinds[g_m]+1)+'Map,(/opSetCore%size*opDat'+str(invinds[g_m]+1)+'MapDim/))')
      for g_m in range(0,nargs):
        if maps[g_m] == OP_ID:
          code('CALL c_f_pointer(opArg'+str(g_m+1)+'%data,opDat'+str(g_m+1)+'Local,(/opDat'+str(g_m+1)+'Cardinality/))')
        elif maps[g_m] == OP_GBL:
          code('CALL c_f_pointer(opArg'+str(g_m+1)+'%data,opDat'+str(g_m+1)+'Local, (/opArg'+str(g_m+1)+'%dim/))')
      code('')
  
      code('')
      if 0:
        code('CALL op_wrap_'+name+'( &')
        if nopts>0:
          code('& optflags, &')
        for g_m in range(0,ninds):
          if invinds[g_m] in needDimList:
            code('& opArg'+str(invinds[g_m]+1)+'%dim, &')
          code('& opDat'+str(invinds[g_m]+1)+'Local, &')
        for g_m in range(0,nargs):
          if maps[g_m] <> OP_MAP:
            if g_m in needDimList:
              code('& opArg'+str(g_m+1)+'%dim, &')
            code('& opDat'+str(g_m+1)+'Local, &')
        if nmaps > 0:
          k = []
          for g_m in range(0,nargs):
            if maps[g_m] == OP_MAP and (not mapnames[g_m] in k):
              k = k + [mapnames[g_m]]
              code('& opDat'+str(invinds[inds[g_m]-1]+1)+'Map, &')
              code('& opDat'+str(invinds[inds[g_m]-1]+1)+'MapDim, &')
        code('& 0, opSetCore%core_size)')
      code('CALL op_mpi_wait_all(numberOfOpDats,opArgArray)')
      code('CALL op_wrap_'+name+'( &')
      if nopts>0:
        code('& optflags, &')
      for g_m in range(0,ninds):
        if invinds[g_m] in needDimList:
          code('& opArg'+str(invinds[g_m]+1)+'%dim, &')
        code('& opDat'+str(invinds[g_m]+1)+'Local, &')
      for g_m in range(0,nargs):
        if maps[g_m] <> OP_MAP:
          if g_m in needDimList:
            code('& opArg'+str(g_m+1)+'%dim, &')
          code('& opDat'+str(g_m+1)+'Local, &')
      if nmaps > 0:
        k = []
        for g_m in range(0,nargs):
          if maps[g_m] == OP_MAP and (not mapnames[g_m] in k):
            k = k + [mapnames[g_m]]
            code('& opDat'+str(invinds[inds[g_m]-1]+1)+'Map, &')
            code('& opDat'+str(invinds[inds[g_m]-1]+1)+'MapDim, &')
      code('& 0, n_upper)')
  #    code('& opSetCore%core_size, n_upper)')
  
  
  #    IF('(n_upper .EQ. 0) .OR. (n_upper .EQ. opSetCore%core_size)')
  #    code('CALL op_mpi_wait_all(numberOfOpDats,opArgArray)')
  #    ENDIF()
  #    code('')
  
  
  
      code('')
      code('CALL op_mpi_set_dirtybit(numberOfOpDats,opArgArray)')
      code('')
  
      #reductions
      for g_m in range(0,nargs):
        if maps[g_m] == OP_GBL and (accs[g_m] == OP_INC or accs[g_m] == OP_MIN or accs[g_m] == OP_MAX or accs[g_m] == OP_WRITE):
          if typs[g_m] == 'real(8)' or typs[g_m] == 'REAL(kind=8)':
            code('CALL op_mpi_reduce_double(opArg'+str(g_m+1)+',opArg'+str(g_m+1)+'%data)')
          elif typs[g_m] == 'real(4)' or typs[g_m] == 'REAL(kind=4)':
            code('CALL op_mpi_reduce_float(opArg'+str(g_m+1)+',opArg'+str(g_m+1)+'%data)')
          elif typs[g_m] == 'integer(4)' or typs[g_m] == 'INTEGER(kind=4)':
            code('CALL op_mpi_reduce_int(opArg'+str(g_m+1)+',opArg'+str(g_m+1)+'%data)')
          elif typs[g_m] == 'logical' or typs[g_m] == 'logical*1':
            code('CALL op_mpi_reduce_bool(opArg'+str(g_m+1)+',opArg'+str(g_m+1)+'%data)')
          code('')
  
      code('call op_timers_core(endTime)')
      code('')
      code('dataTransfer = 0.0')
      if ninds == 0:
        for g_m in range(0,nargs):
          if optflags[g_m] == 1:
                IF('opArg'+str(g_m+1)+'%opt == 1')
          if accs[g_m] == OP_READ or accs[g_m] == OP_WRITE:
            if maps[g_m] == OP_GBL:
              code('dataTransfer = dataTransfer + opArg'+str(g_m+1)+'%size')
            else:
              code('dataTransfer = dataTransfer + opArg'+str(g_m+1)+'%size * opSetCore%size')
          else:
            if maps[g_m] == OP_GBL:
              code('dataTransfer = dataTransfer + opArg'+str(g_m+1)+'%size * 2.d0')
            else:
              code('dataTransfer = dataTransfer + opArg'+str(g_m+1)+'%size * opSetCore%size * 2.d0')
          if optflags[g_m] == 1:
            ENDIF()
      else:
        names = []
        for g_m in range(0,ninds):
          mult=''
          if indaccs[g_m] <> OP_WRITE and indaccs[g_m] <> OP_READ:
            mult = ' * 2.d0'
          if not var[invinds[g_m]] in names:
            if optflags[invinds[g_m]] == 1:
              IF('opArg'+str(g_m+1)+'%opt == 1')
            code('dataTransfer = dataTransfer + opArg'+str(invinds[g_m]+1)+'%size * MIN(n_upper,getSetSizeFromOpArg(opArg'+str(invinds[g_m]+1)+'))'+mult)
            names = names + [var[invinds[g_m]]]
            if optflags[invinds[g_m]] == 1:
              ENDIF()
        for g_m in range(0,nargs):
          mult=''
          if accs[g_m] <> OP_WRITE and accs[g_m] <> OP_READ:
            mult = ' * 2.d0'
          if not var[g_m] in names:
            if optflags[g_m] == 1:
              IF('opArg'+str(g_m+1)+'%opt == 1')
            names = names + [var[invinds[g_m]]]
            if maps[g_m] == OP_ID:
              code('dataTransfer = dataTransfer + opArg'+str(g_m+1)+'%size * MIN(n_upper,getSetSizeFromOpArg(opArg'+str(g_m+1)+'))'+mult)
            elif maps[g_m] == OP_GBL:
              code('dataTransfer = dataTransfer + opArg'+str(g_m+1)+'%size'+mult)
            if optflags[g_m] == 1:
              ENDIF()
        if nmaps > 0:
          k = []
          for g_m in range(0,nargs):
            if maps[g_m] == OP_MAP and (not mapnames[g_m] in k):
              k = k + [mapnames[g_m]]
              code('dataTransfer = dataTransfer + n_upper * opDat'+str(invinds[inds[g_m]-1]+1)+'MapDim * 4.d0')
  
      code('returnSetKernelTiming = setKernelTime('+str(nk)+' , userSubroutine//C_NULL_CHAR, &')
      code('& endTime-startTime, dataTransfer, 0.00000_4, 1)')
      code('')
      #code('returnSetKernelTiming = setKernelTime('+str(nk)+' , userSubroutine//C_NULL_CHAR, &')
      #code('& endTime-startTime,0.00000,0.00000, 1)')
      depth = depth - 2
      code('END SUBROUTINE')
      code('END MODULE')
      code('')
  
  ##########################################################################
  #  output individual kernel file
  ##########################################################################
      if hydra:
        name = 'kernels/'+kernels[nk]['master_file']+'/'+name
        #fid = open(name+'_seqkernel.F95','w')
        fid = open(name+'_veckernel.F95','w')
      elif bookleaf:
        fid = open(prefixes[prefix_i]+name+'_veckernel.F90','w')
      else:
        #fid = open(name+'_seqkernel.F90','w')
        fid = open(name+'_veckernel.F90','w')
      date = datetime.datetime.now()
      fid.write('!\n! auto-generated by op2_fortran.py\n!\n\n')
      fid.write(file_text.strip())
      fid.close()
  