[mpas-developers] additions to MPAS registry and I/O

Michael Duda duda at ucar.edu
Tue May 4 16:35:50 MDT 2010


Hi, Developers.

Currently, it's not possible to define a time-invariant, real
scalar field in the registry; however, I've got a set of changes
to the registry and framework code that enable this, and I've
tested them in the context of the mpas_cam_coupling branch.
Essentially, the changes would make it possible to define
something like

      var real    ptop ( ) iro ptop - -

that would be accessible in the code as 

   block % mesh % ptop % scalar

Attached are the updated versions of three files as I'd like to
commit them to the trunk. If anyone has any objections or
suggestions for improvements, please feel free to comment.

Cheers,
Michael
-------------- next part --------------
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dictionary.h"
#include "registry_types.h"
#include "gen_inc.h"
#include "fortprintf.h"

int is_derived_dim(char * d)
{
   if (strchr(d, (int)'+')) return 1;
   if (strchr(d, (int)'-')) return 1;

   return 0;
}

void split_derived_dim_string(char * dim, char ** p1, char ** p2)
{
   char * cp, * cm, * c;
   int n;

   cp = strchr(dim, (int)'+');
   cm = strchr(dim, (int)'-');
   if (!cp) 
      c = cm;
   else if (!cm) 
      c = cp;
   else if (cm < cp) 
      c = cm;
   else 
      c = cp;

   n = c - dim;
   *p1 = (char *)malloc(n*sizeof(char));
   snprintf(*p1, n, "%s", dim+1);

   *p2 = (char *)malloc((strlen(dim)-n+1)*sizeof(char));
   sprintf(*p2, "%s", dim+n);
}

void gen_namelists(struct namelist * nls)
{
   struct namelist * nls_ptr;
   struct dtable * dictionary;
   int done;
   char nlrecord[1024];
   FILE * fd;

   /*
    *  Generate config_type.inc
    */
   fd = fopen("config_defs.inc", "w");

   nls_ptr = nls;
   while (nls_ptr) {
      if (nls_ptr->vtype == INTEGER)   fortprintf(fd, "   integer :: %s\n",nls_ptr->name);
      if (nls_ptr->vtype == REAL)      fortprintf(fd, "   real (KIND=RKIND) :: %s\n",nls_ptr->name);
      if (nls_ptr->vtype == LOGICAL)   fortprintf(fd, "   logical :: %s\n",nls_ptr->name);
      if (nls_ptr->vtype == CHARACTER) fortprintf(fd, "   character (len=32) :: %s\n",nls_ptr->name);

      nls_ptr = nls_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate namelist_defs.inc
    */
   fd = fopen("config_namelist_defs.inc", "w");
   dict_alloc(&dictionary);

   done = 0;
  
   while (!done) {
      nls_ptr = nls;
      while (nls_ptr && dict_search(dictionary, nls_ptr->record))
         nls_ptr = nls_ptr->next;

      if (nls_ptr) {
         dict_insert(dictionary, nls_ptr->record);
         strncpy(nlrecord, nls_ptr->record, 1024);
         fortprintf(fd, "      namelist /%s/ %s", nls_ptr->record, nls_ptr->name);
         nls_ptr = nls_ptr->next;
         while(nls_ptr) {
            if (strncmp(nls_ptr->record, nlrecord, 1024) == 0)
               fortprintf(fd, ", &\n                    %s", nls_ptr->name);
            nls_ptr = nls_ptr->next;
         }
         fortprintf(fd, "\n");
      }
      else
         done = 1;
   }
   

   dict_free(&dictionary);
   fclose(fd);


   /*
    *  Generate namelist_reads.inc
    */
   fd = fopen("config_set_defaults.inc", "w");
   nls_ptr = nls;
   while (nls_ptr) {
      if (nls_ptr->vtype == INTEGER) fortprintf(fd, "      %s = %i\n", nls_ptr->name, nls_ptr->defval.ival);
      if (nls_ptr->vtype == REAL)    fortprintf(fd, "      %s = %f\n", nls_ptr->name, nls_ptr->defval.rval);
      if (nls_ptr->vtype == LOGICAL) {
         if (nls_ptr->defval.lval == 0) 
            fortprintf(fd, "      %s = .false.\n", nls_ptr->name);
         else
            fortprintf(fd, "      %s = .true.\n", nls_ptr->name);
      }
      if (nls_ptr->vtype == CHARACTER)
         fortprintf(fd, "      %s = \"%s\"\n", nls_ptr->name, nls_ptr->defval.cval);
      nls_ptr = nls_ptr->next;
   }
   fortprintf(fd, "\n");
   fclose(fd);


   fd = fopen("config_namelist_reads.inc", "w");
   dict_alloc(&dictionary);
   nls_ptr = nls;
   while (nls_ptr) {
      if (!dict_search(dictionary, nls_ptr->record)) {
         fortprintf(fd, "         read(funit,%s)\n", nls_ptr->record);
         dict_insert(dictionary, nls_ptr->record);
      }
      nls_ptr = nls_ptr->next;
   }
   fortprintf(fd, "\n");

   dict_free(&dictionary);
   fclose(fd);


   fd = fopen("config_bcast_namelist.inc", "w");
   nls_ptr = nls;
   while (nls_ptr) {
      if (nls_ptr->vtype == INTEGER)   fortprintf(fd, "      call dmpar_bcast_int(dminfo, %s)\n", nls_ptr->name);
      if (nls_ptr->vtype == REAL)      fortprintf(fd, "      call dmpar_bcast_real(dminfo, %s)\n", nls_ptr->name);
      if (nls_ptr->vtype == LOGICAL)   fortprintf(fd, "      call dmpar_bcast_logical(dminfo, %s)\n", nls_ptr->name);
      if (nls_ptr->vtype == CHARACTER) fortprintf(fd, "      call dmpar_bcast_char(dminfo, %s)\n", nls_ptr->name);
      nls_ptr = nls_ptr->next;
   }
   fortprintf(fd, "\n");
   fclose(fd);
}


void gen_field_defs(struct variable * vars, struct dimension * dims)
{
   struct variable * var_ptr;
   struct variable * var_ptr2;
   struct dimension * dim_ptr;
   struct dimension_list * dimlist_ptr;
   FILE * fd;
   char super_array[1024];
   char array_class[1024];
   int i;
   int class_start, class_end;
   int vtype;

   /*
    * Generate indices for super arrays
    */
   fd = fopen("super_array_indices.inc", "w");
   var_ptr = vars;
   memcpy(super_array, var_ptr->super_array, 1024);
   i = 1;
   while (var_ptr) {
      if (strncmp(super_array, var_ptr->super_array, 1024) != 0) {
         memcpy(super_array, var_ptr->super_array, 1024);
         i = 1;
       }
      if (strncmp(var_ptr->array_class, "-", 1024) != 0) fortprintf(fd, "      integer :: index_%s = %i\n", var_ptr->name_in_code, i++);
      var_ptr = var_ptr->next;
   }
   var_ptr = vars;
   memcpy(super_array, var_ptr->super_array, 1024);
   memcpy(array_class, var_ptr->array_class, 1024);
   class_start = 1;
   class_end = 1;
   i = 1;
   while (var_ptr) {
      if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
         if (strncmp(super_array, var_ptr->super_array, 1024) != 0) {
            if (strncmp(super_array, "-", 1024) != 0) fortprintf(fd, "      integer :: %s_end = %i\n", array_class, class_end);
            if (strncmp(super_array, "-", 1024) != 0) fortprintf(fd, "      integer :: num_%s = %i\n", super_array, i);
            class_start = 1;
            class_end = 1;
            i = 1;
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            fortprintf(fd, "      integer :: %s_start = %i\n", array_class, class_start);
         }
         else if (strncmp(array_class, var_ptr->array_class, 1024) != 0) {
            fortprintf(fd, "      integer :: %s_end = %i\n", array_class, class_end);
            class_start = class_end+1;
            class_end = class_start;
            memcpy(array_class, var_ptr->array_class, 1024);
            fortprintf(fd, "      integer :: %s_start = %i\n", array_class, class_start);
            i++;
         }
         else {
            class_end++;
            i++;
         }
      }
      var_ptr = var_ptr->next;
   }
   if (strncmp(super_array, "-", 1024) != 0) fortprintf(fd, "      integer :: %s_end = %i\n", array_class, class_end);
   if (strncmp(super_array, "-", 1024) != 0) fortprintf(fd, "      integer :: num_%s = %i\n", super_array, i);
   fclose(fd);


   /*
    *  Generate declarations of dimensions
    */
   fd = fopen("field_dimensions.inc", "w");
   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      integer :: %s\n", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      integer :: %sSolve\n", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate dummy dimension argument list
    */
   fd = fopen("dim_dummy_args.inc", "w");
   dim_ptr = dims;
   if (dim_ptr && dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
      fortprintf(fd, "                            %s", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, " &\n");

   fclose(fd);


   /*
    *  Generate dummy dimension argument declaration list
    */
   fd = fopen("dim_dummy_decls.inc", "w");
   dim_ptr = dims;
   if (dim_ptr && dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
      fortprintf(fd, "      integer, intent(in) :: %s", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %s", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   fclose(fd);


   /*
    *  Generate declarations of dimensions
    */
   fd = fopen("dim_decls.inc", "w");
   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      integer :: %s\n", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate calls to read dimensions from input file
    */
   fd = fopen("read_dims.inc", "w");
   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      call io_input_get_dimension(input_obj, \'%s\', %s)\n", dim_ptr->name_in_file, dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate declarations of time-invariant fields
    */
   fd = fopen("time_invariant_fields.inc", "w");
   var_ptr = vars;
   while (var_ptr) {
      if (var_ptr->timedim == 0) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            if (vtype == INTEGER)  fortprintf(fd, "      type (field%idInteger), pointer :: %s\n", (var_ptr2->ndims)+1, var_ptr2->super_array);
            if (vtype == REAL)     fortprintf(fd, "      type (field%idReal), pointer :: %s\n", (var_ptr2->ndims)+1, var_ptr2->super_array);
         }
         else {
            if (var_ptr->vtype == INTEGER)  fortprintf(fd, "      type (field%idInteger), pointer :: %s\n", var_ptr->ndims, var_ptr->name_in_code);
            if (var_ptr->vtype == REAL)     fortprintf(fd, "      type (field%idReal), pointer :: %s\n", var_ptr->ndims, var_ptr->name_in_code);
            var_ptr = var_ptr->next;
         }
      }
      else
         var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate declarations of time-invariant fields
    */
   fd = fopen("time_varying_fields.inc", "w");
   var_ptr = vars;
   while (var_ptr) {
      if (var_ptr->timedim == 1) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            if (vtype == INTEGER)  fortprintf(fd, "      type (field%idInteger), pointer :: %s\n", (var_ptr2->ndims)+1, var_ptr2->super_array);
            if (vtype == REAL)     fortprintf(fd, "      type (field%idReal), pointer :: %s\n", (var_ptr2->ndims)+1, var_ptr2->super_array);
         }
         else {
            if (var_ptr->vtype == INTEGER)  fortprintf(fd, "      type (field%idInteger), pointer :: %s\n", var_ptr->ndims, var_ptr->name_in_code);
            if (var_ptr->vtype == REAL)     fortprintf(fd, "      type (field%idReal), pointer :: %s\n", var_ptr->ndims, var_ptr->name_in_code);
            var_ptr = var_ptr->next;
         }
      }
      else
         var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate grid metadata allocations
    */
   fd = fopen("grid_meta_allocs.inc", "w");

   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      g %% %s = %s\n", dim_ptr->name_in_code, dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   var_ptr = vars;
   while (var_ptr) {
      if (var_ptr->timedim == 0) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            i = 0;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               i++;
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            fortprintf(fd, "      allocate(g %% %s)\n", var_ptr2->super_array);
            fortprintf(fd, "      allocate(g %% %s %% ioinfo)\n", var_ptr2->super_array);
            fortprintf(fd, "      allocate(g %% %s %% array(%i, ", var_ptr2->super_array, i);
            dimlist_ptr = var_ptr2->dimlist;
            fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
            dimlist_ptr = dimlist_ptr->next;
            while (dimlist_ptr) {
               fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
               dimlist_ptr = dimlist_ptr->next;
            }
            fortprintf(fd, "))\n");

            if (var_ptr2->iostreams & INPUT0) 
               fortprintf(fd, "      g %% %s %% ioinfo %% input = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      g %% %s %% ioinfo %% input = .false.\n", var_ptr2->super_array);

            if (var_ptr2->iostreams & RESTART0) 
               fortprintf(fd, "      g %% %s %% ioinfo %% restart = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      g %% %s %% ioinfo %% restart = .false.\n", var_ptr2->super_array);

            if (var_ptr2->iostreams & OUTPUT0) 
               fortprintf(fd, "      g %% %s %% ioinfo %% output = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      g %% %s %% ioinfo %% output = .false.\n", var_ptr2->super_array);
            fortprintf(fd, "\n");
         }
         else {
            fortprintf(fd, "      allocate(g %% %s)\n", var_ptr->name_in_code);
            fortprintf(fd, "      allocate(g %% %s %% ioinfo)\n", var_ptr->name_in_code);
            if (var_ptr->ndims > 0) {
               fortprintf(fd, "      allocate(g %% %s %% array(", var_ptr->name_in_code);
               dimlist_ptr = var_ptr->dimlist;
               fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
               dimlist_ptr = dimlist_ptr->next;
               while (dimlist_ptr) {
                  fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
                  dimlist_ptr = dimlist_ptr->next;
               }
               fortprintf(fd, "))\n");
   
               if (var_ptr->iostreams & INPUT0) 
                  fortprintf(fd, "      g %% %s %% ioinfo %% input = .true.\n", var_ptr->name_in_code);
               else
                  fortprintf(fd, "      g %% %s %% ioinfo %% input = .false.\n", var_ptr->name_in_code);
   
               if (var_ptr->iostreams & RESTART0) 
                  fortprintf(fd, "      g %% %s %% ioinfo %% restart = .true.\n", var_ptr->name_in_code);
               else
                  fortprintf(fd, "      g %% %s %% ioinfo %% restart = .false.\n", var_ptr->name_in_code);
   
               if (var_ptr->iostreams & OUTPUT0) 
                  fortprintf(fd, "      g %% %s %% ioinfo %% output = .true.\n", var_ptr->name_in_code);
               else
                  fortprintf(fd, "      g %% %s %% ioinfo %% output = .false.\n", var_ptr->name_in_code);
               fortprintf(fd, "\n");
            }
            else { 
               if (var_ptr->iostreams & INPUT0) 
                  fortprintf(fd, "      g %% %s %% ioinfo %% input = .true.\n", var_ptr->name_in_code);
               else
                  fortprintf(fd, "      g %% %s %% ioinfo %% input = .false.\n", var_ptr->name_in_code);
   
               if (var_ptr->iostreams & RESTART0) 
                  fortprintf(fd, "      g %% %s %% ioinfo %% restart = .true.\n", var_ptr->name_in_code);
               else
                  fortprintf(fd, "      g %% %s %% ioinfo %% restart = .false.\n", var_ptr->name_in_code);
   
               if (var_ptr->iostreams & OUTPUT0) 
                  fortprintf(fd, "      g %% %s %% ioinfo %% output = .true.\n", var_ptr->name_in_code);
               else
                  fortprintf(fd, "      g %% %s %% ioinfo %% output = .false.\n", var_ptr->name_in_code);
               fortprintf(fd, "\n");
            }
            var_ptr = var_ptr->next;
         }
      }
      else
         var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate grid metadata deallocations
    */
   fd = fopen("grid_meta_deallocs.inc", "w");

   var_ptr = vars;
   while (var_ptr) {
      if (var_ptr->timedim == 0) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            i = 0;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               i++;
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            fortprintf(fd, "      deallocate(g %% %s %% array)\n", var_ptr2->super_array);
            fortprintf(fd, "      deallocate(g %% %s %% ioinfo)\n", var_ptr2->super_array);
            fortprintf(fd, "      deallocate(g %% %s)\n\n", var_ptr2->super_array);
         }
         else {
            if (var_ptr->ndims > 0) {
               fortprintf(fd, "      deallocate(g %% %s %% array)\n", var_ptr->name_in_code);
               fortprintf(fd, "      deallocate(g %% %s %% ioinfo)\n", var_ptr->name_in_code);
               fortprintf(fd, "      deallocate(g %% %s)\n\n", var_ptr->name_in_code);
            }
            else {
               fortprintf(fd, "      deallocate(g %% %s %% ioinfo)\n", var_ptr->name_in_code);
               fortprintf(fd, "      deallocate(g %% %s)\n\n", var_ptr->name_in_code);
            }
            var_ptr = var_ptr->next;
         }
      }
      else
         var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate grid state allocations
    */
   fd = fopen("grid_state_allocs.inc", "w");

   var_ptr = vars;
   while (var_ptr) {
      if (var_ptr->timedim == 1 && var_ptr->ndims > 0) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            i = 0;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               i++;
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            fortprintf(fd, "      allocate(s %% %s)\n", var_ptr2->super_array);
            fortprintf(fd, "      allocate(s %% %s %% ioinfo)\n", var_ptr2->super_array);
            fortprintf(fd, "      allocate(s %% %s %% array(%i, ", var_ptr2->super_array, i);
            dimlist_ptr = var_ptr2->dimlist;
            fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
            dimlist_ptr = dimlist_ptr->next;
            while (dimlist_ptr) {
               fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
               dimlist_ptr = dimlist_ptr->next;
            }
            fortprintf(fd, "))\n");
            fortprintf(fd, "      s %% %s %% block => b\n", var_ptr2->super_array);
   
            if (var_ptr2->iostreams & INPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .false.\n", var_ptr2->super_array);
   
            if (var_ptr2->iostreams & RESTART0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .false.\n", var_ptr2->super_array);
   
            if (var_ptr2->iostreams & OUTPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .false.\n", var_ptr2->super_array);
            fortprintf(fd, "\n");
         }
         else {
            fortprintf(fd, "      allocate(s %% %s)\n", var_ptr->name_in_code);
            fortprintf(fd, "      allocate(s %% %s %% ioinfo)\n", var_ptr->name_in_code);
            fortprintf(fd, "      allocate(s %% %s %% array(", var_ptr->name_in_code);
            dimlist_ptr = var_ptr->dimlist;
            fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
            dimlist_ptr = dimlist_ptr->next;
            while (dimlist_ptr) {
               fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
               dimlist_ptr = dimlist_ptr->next;
            }
            fortprintf(fd, "))\n");
            fortprintf(fd, "      s %% %s %% block => b\n", var_ptr->name_in_code);
   
            if (var_ptr->iostreams & INPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .true.\n", var_ptr->name_in_code);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .false.\n", var_ptr->name_in_code);
   
            if (var_ptr->iostreams & RESTART0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .true.\n", var_ptr->name_in_code);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .false.\n", var_ptr->name_in_code);
   
            if (var_ptr->iostreams & OUTPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .true.\n", var_ptr->name_in_code);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .false.\n", var_ptr->name_in_code);
            fortprintf(fd, "\n");
            var_ptr = var_ptr->next;
         }
      }
      else if (var_ptr->timedim == 1) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            i = 0;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               i++;
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            fortprintf(fd, "      allocate(s %% %s)\n", var_ptr2->super_array);
            fortprintf(fd, "      allocate(s %% %s %% ioinfo)\n", var_ptr2->super_array);
            fortprintf(fd, "      allocate(s %% %s %% array(%i)", var_ptr->name_in_code, i);
            fortprintf(fd, "      s %% %s %% block => b\n", var_ptr2->super_array);
   
            if (var_ptr2->iostreams & INPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .false.\n", var_ptr2->super_array);
   
            if (var_ptr2->iostreams & RESTART0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .false.\n", var_ptr2->super_array);
   
            if (var_ptr2->iostreams & OUTPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .true.\n", var_ptr2->super_array);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .false.\n", var_ptr2->super_array);
            fortprintf(fd, "\n");
         }
         else {
            fortprintf(fd, "      allocate(s %% %s)\n", var_ptr->name_in_code);
            fortprintf(fd, "      allocate(s %% %s %% ioinfo)\n", var_ptr->name_in_code);
            fortprintf(fd, "      s %% %s %% block => b\n", var_ptr->name_in_code);
   
            if (var_ptr->iostreams & INPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .true.\n", var_ptr->name_in_code);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% input = .false.\n", var_ptr->name_in_code);
   
            if (var_ptr->iostreams & RESTART0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .true.\n", var_ptr->name_in_code);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% restart = .false.\n", var_ptr->name_in_code);
   
            if (var_ptr->iostreams & OUTPUT0) 
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .true.\n", var_ptr->name_in_code);
            else
               fortprintf(fd, "      s %% %s %% ioinfo %% output = .false.\n", var_ptr->name_in_code);
            fortprintf(fd, "\n");
            var_ptr = var_ptr->next;
         }
      } 
      else
         var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate grid state deallocations
    */
   fd = fopen("grid_state_deallocs.inc", "w");

   var_ptr = vars;
   while (var_ptr) {
      if (var_ptr->timedim == 1) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            i = 0;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               i++;
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            fortprintf(fd, "      deallocate(s %% %s %% array)\n", var_ptr2->super_array);
            fortprintf(fd, "      deallocate(s %% %s %% ioinfo)\n", var_ptr2->super_array);
            fortprintf(fd, "      deallocate(s %% %s)\n\n", var_ptr2->super_array);
         }
         else {
            if (var_ptr->ndims > 0) fortprintf(fd, "      deallocate(s %% %s %% array)\n", var_ptr->name_in_code);
            fortprintf(fd, "      deallocate(s %% %s %% ioinfo)\n", var_ptr->name_in_code);
            fortprintf(fd, "      deallocate(s %% %s)\n\n", var_ptr->name_in_code);
            var_ptr = var_ptr->next;
         }
      }
      else
         var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate copies of state arrays
    */
   fd = fopen("copy_state.inc", "w");

   var_ptr = vars;
   while (var_ptr) {
      if (var_ptr->timedim == 1) {
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            memcpy(super_array, var_ptr->super_array, 1024);
            memcpy(array_class, var_ptr->array_class, 1024);
            vtype = var_ptr->vtype;
            i = 0;
            while (var_ptr && strncmp(super_array, var_ptr->super_array, 1024) == 0) {
               i++;
               var_ptr2 = var_ptr;
               var_ptr = var_ptr->next;
            }
            if (var_ptr2->ndims > 0) 
               fortprintf(fd, "      dest %% %s %% array = src %% %s %% array\n", var_ptr2->super_array, var_ptr2->super_array);
            else
               fortprintf(fd, "      dest %% %s %% scalar = src %% %s %% scalar\n", var_ptr2->super_array, var_ptr2->super_array);
         }
         else {
            if (var_ptr->ndims > 0) 
               fortprintf(fd, "      dest %% %s %% array = src %% %s %% array\n", var_ptr->name_in_code, var_ptr->name_in_code);
            else
               fortprintf(fd, "      dest %% %s %% scalar = src %% %s %% scalar\n", var_ptr->name_in_code, var_ptr->name_in_code);
            var_ptr = var_ptr->next;
         }
      }
      else
         var_ptr = var_ptr->next;
   }

   fclose(fd);
}


void gen_reads(struct variable * vars, struct dimension * dims)
{
   struct variable * var_ptr;
   struct dimension * dim_ptr;
   struct dimension_list * dimlist_ptr, * lastdim;
   struct dtable * dictionary;
   FILE * fd;
   char vtype[5];
   char fname[32];
   char * cp1, * cp2;
   int i, j;
   int ivtype;
   int has_vert_dim, vert_dim;


   /*
    *  Generate declarations of IDs belonging in io_input_object
    */
   fd = fopen("io_input_obj_decls.inc", "w");

   fortprintf(fd, "      integer :: rdDimIDTime\n");
   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      integer :: rdDimID%s\n", dim_ptr->name_in_file);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   fortprintf(fd, "      integer :: rdLocalTime\n");
   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      integer :: rdLocal%s\n", dim_ptr->name_in_file);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   var_ptr = vars;
   while (var_ptr) {
      fortprintf(fd, "      integer :: rdVarID%s\n", var_ptr->name_in_file);
      var_ptr = var_ptr->next;
   }

   fclose(fd);
   

   /*
    *  Generate read and distribute code
    */
   fd = fopen("io_input_fields.inc", "w");

   var_ptr = vars;
   while (var_ptr) {
      i = 1;
      dimlist_ptr = var_ptr->dimlist;
      if (var_ptr->vtype == INTEGER) sprintf(vtype, "int"); 
      else if (var_ptr->vtype == REAL) sprintf(vtype, "real"); 

      if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
         if (var_ptr->timedim) {
            fortprintf(fd, "      if ((block %% time_levs(1) %% state %% %s %% ioinfo %% input .and. .not. config_do_restart) .or. &\n", var_ptr->super_array);
            fortprintf(fd, "          (block %% time_levs(1) %% state %% %s %% ioinfo %% restart .and. config_do_restart)) then\n", var_ptr->super_array);
         }
         else {
            fortprintf(fd, "      if ((block %% mesh %% %s %% ioinfo %% input .and. .not. config_do_restart) .or. &\n", var_ptr->super_array);
            fortprintf(fd, "          (block %% mesh %% %s %% ioinfo %% restart .and. config_do_restart)) then\n", var_ptr->super_array);
         }
      }
      else {
         if (var_ptr->timedim) {
            fortprintf(fd, "      if ((block %% time_levs(1) %% state %% %s %% ioinfo %% input .and. .not. config_do_restart) .or. &\n", var_ptr->name_in_code);
            fortprintf(fd, "          (block %% time_levs(1) %% state %% %s %% ioinfo %% restart .and. config_do_restart)) then\n", var_ptr->name_in_code);
         }
         else {
            fortprintf(fd, "      if ((block %% mesh %% %s %% ioinfo %% input .and. .not. config_do_restart) .or. &\n", var_ptr->name_in_code);
            fortprintf(fd, "          (block %% mesh %% %s %% ioinfo %% restart .and. config_do_restart)) then\n", var_ptr->name_in_code);
         }
      }
      vert_dim = 0;
      while (dimlist_ptr) {
            if (i < var_ptr->ndims) {
               has_vert_dim = !strcmp( "nVertLevels", dimlist_ptr->dim->name_in_code);
               fortprintf(fd, "      %s%id %% ioinfo %% start(%i) = 1\n", vtype, var_ptr->ndims, i);
               if (has_vert_dim) {
                  vert_dim = i;
                  fortprintf(fd, "#ifdef EXPAND_LEVELS\n");
                  fortprintf(fd, "      if (.not. config_do_restart) then\n");
                  fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = 1\n", vtype, var_ptr->ndims, i);
                  fortprintf(fd, "      else\n");
                  fortprintf(fd, "#endif\n");
               }
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = block %% mesh %% %s\n", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = %s\n", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
               if (has_vert_dim) {
                  fortprintf(fd, "#ifdef EXPAND_LEVELS\n");
                  fortprintf(fd, "      end if\n");
                  fortprintf(fd, "#endif\n");
               }
            }
            else {
               if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
                  split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
                  fortprintf(fd, "      %s%id %% ioinfo %% start(%i) = read%sStart\n", vtype, var_ptr->ndims, i, cp1);
                  fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = read%sCount%s\n", vtype, var_ptr->ndims, i, cp1, cp2);
                  free(cp1);
                  free(cp2);
               }
               else {
                  fortprintf(fd, "      %s%id %% ioinfo %% start(%i) = read%sStart\n", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code+1);
                  fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = read%sCount\n", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code+1);
               }
            }
         dimlist_ptr = dimlist_ptr->next;
         i++;
      }

      if (var_ptr->ndims > 0) {
         fortprintf(fd, "      allocate(%s%id %% array(", vtype, var_ptr->ndims);
         i = 1;
         dimlist_ptr = var_ptr->dimlist;
   
         if (i < var_ptr->ndims) {
            if (dimlist_ptr->dim->constant_value < 0)
               fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
            else
               fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
         }
         else {
            if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
               split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
               fortprintf(fd, "read%sCount%s", cp1, cp2);
               free(cp1);
               free(cp2);
            }
            else
               fortprintf(fd, "read%sCount", dimlist_ptr->dim->name_in_code+1);
         }
    
         dimlist_ptr = dimlist_ptr->next;
         i++;
         while (dimlist_ptr) {
            if (i < var_ptr->ndims) {
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
            }
            else {
               if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
                  split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
                  fortprintf(fd, ", read%sCount%s", cp1, cp2);
                  free(cp1);
                  free(cp2);
               }
               else
                  fortprintf(fd, ", read%sCount", dimlist_ptr->dim->name_in_code+1);
            }
            dimlist_ptr = dimlist_ptr->next;
            i++;
         }
         fortprintf(fd, "))\n\n");

         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            fortprintf(fd, "      allocate(super_%s%id(", vtype, var_ptr->ndims);
            i = 1;
            dimlist_ptr = var_ptr->dimlist;
      
            if (i < var_ptr->ndims) {
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, "block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
            }
            dimlist_ptr = dimlist_ptr->next;
            i++;
            while (dimlist_ptr) {
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
               dimlist_ptr = dimlist_ptr->next;
               i++;
            }
            fortprintf(fd, "))\n\n");
         }
      }

      fortprintf(fd, "      %s%id %% ioinfo %% fieldName = \'%s\'\n", vtype, var_ptr->ndims, var_ptr->name_in_file);
      if (var_ptr->timedim)
         fortprintf(fd, "      call io_input_field_time(input_obj, %s%id)\n", vtype, var_ptr->ndims);
      else
         fortprintf(fd, "      call io_input_field(input_obj, %s%id)\n", vtype, var_ptr->ndims);

      if (vert_dim > 0) {
         fortprintf(fd, "#ifdef EXPAND_LEVELS\n");
         fortprintf(fd, "      if (.not. config_do_restart) then\n");
         fortprintf(fd, "         do k=2,EXPAND_LEVELS\n");
         fortprintf(fd, "            %s%id %% array(", vtype, var_ptr->ndims);
         for (i=1; i<=var_ptr->ndims; i++) {
            if (i > 1) fortprintf(fd, ",");
            fortprintf(fd, "%s", i == vert_dim ? "k" : ":");
         }
         fortprintf(fd, ") = %s%id %% array(", vtype, var_ptr->ndims);
         for (i=1; i<=var_ptr->ndims; i++) {
            if (i > 1) fortprintf(fd, ",");
            fortprintf(fd, "%s", i == vert_dim ? "1" : ":");
         }
         fortprintf(fd, ")\n");
         fortprintf(fd, "         end do\n");
         fortprintf(fd, "      end if\n");
         fortprintf(fd, "#endif\n");
      }

      if (var_ptr->ndims > 0) {
         fortprintf(fd, "      call dmpar_alltoall_field(dminfo, &\n");
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            if (var_ptr->timedim) 
               fortprintf(fd, "                                %s%id %% array, super_%s%id, &\n", vtype, var_ptr->ndims, vtype, var_ptr->ndims);
            else
               fortprintf(fd, "                                %s%id %% array, super_%s%id, &\n", vtype, var_ptr->ndims, vtype, var_ptr->ndims);
         }
         else {
            if (var_ptr->timedim) 
               fortprintf(fd, "                                %s%id %% array, block %% time_levs(1) %% state %% %s %% array, &\n", vtype, var_ptr->ndims, var_ptr->name_in_code);
            else
               fortprintf(fd, "                                %s%id %% array, block %% mesh %% %s %% array, &\n", vtype, var_ptr->ndims, var_ptr->name_in_code);
         }
   
         i = 1;
         dimlist_ptr = var_ptr->dimlist;
         
         if (i < var_ptr->ndims)
            if (dimlist_ptr->dim->constant_value < 0)
               fortprintf(fd, "                                block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
            else
               fortprintf(fd, "                                %s", dimlist_ptr->dim->name_in_code);
         else {
            lastdim = dimlist_ptr;
            if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
               split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
               fortprintf(fd, "                                read%sCount%s", cp1, cp2);
               free(cp1);
               free(cp2);
            }
            else
               fortprintf(fd, "                                read%sCount", dimlist_ptr->dim->name_in_code+1);
         }
    
         dimlist_ptr = dimlist_ptr->next;
         i++;
         while (dimlist_ptr) {
            if (i < var_ptr->ndims)
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, ", block %% mesh %% %s", dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
            else {
               lastdim = dimlist_ptr;
               if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
                  split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
                  fortprintf(fd, ", read%sCount%s", cp1, cp2);
                  free(cp1);
                  free(cp2);
               }
               else
                  fortprintf(fd, ", read%sCount", dimlist_ptr->dim->name_in_code+1);
            }
            dimlist_ptr = dimlist_ptr->next;
            i++;
         }
         fortprintf(fd, ", block %% mesh %% %s, &\n", lastdim->dim->name_in_code);
   
         if (is_derived_dim(lastdim->dim->name_in_code)) {
            split_derived_dim_string(lastdim->dim->name_in_code, &cp1, &cp2);
            fortprintf(fd, "                                send%sList, recv%sList)\n", cp1, cp1);
            free(cp1);
            free(cp2);
         }
         else
            fortprintf(fd, "                                send%sList, recv%sList)\n", lastdim->dim->name_in_code+1, lastdim->dim->name_in_code+1);


         /* Copy from super_ array to field */
         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            i = 1;
            dimlist_ptr = var_ptr->dimlist;
            while (i <= var_ptr->ndims) {
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, "      do i%i=1,block %% mesh %% %s\n", i, dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, "      do i%i=1,%s\n", i, dimlist_ptr->dim->name_in_code);
   
               i++;
               dimlist_ptr = dimlist_ptr->next;
            }
   
            if (var_ptr->timedim) 
               fortprintf(fd, "         block %% time_levs(1) %% state %% %s %% array(index_%s,", var_ptr->super_array, var_ptr->name_in_code);
            else
               fortprintf(fd, "         block %% mesh %% %s %% array(index_%s,", var_ptr->super_array, var_ptr->name_in_code);
            for(i=1; i<=var_ptr->ndims; i++) {
               fortprintf(fd, "i%i",i);
               if (i < var_ptr->ndims) fortprintf(fd, ",");
            }
            fortprintf(fd, ") = super_%s%id(", vtype, var_ptr->ndims);
            for(i=1; i<=var_ptr->ndims; i++) {
               fortprintf(fd, "i%i",i);
               if (i < var_ptr->ndims) fortprintf(fd, ",");
            }
            fortprintf(fd, ")\n");
   
            i = 1;
            while (i <= var_ptr->ndims) {
               fortprintf(fd, "      end do\n");
               i++;
            }
         }

         fortprintf(fd, "      deallocate(%s%id %% array)\n", vtype, var_ptr->ndims);
         if (strncmp(var_ptr->super_array, "-", 1024) != 0)
            fortprintf(fd, "      deallocate(super_%s%id)\n", vtype, var_ptr->ndims);
      }
      else {
         if (var_ptr->timedim) 
            fortprintf(fd, "      block %% time_levs(1) %% state %% %s %% scalar = %s%id %% scalar\n", var_ptr->name_in_code, vtype, var_ptr->ndims);
         else
            fortprintf(fd, "      block %% mesh %% %s %% scalar = %s%id %% scalar\n", var_ptr->name_in_code, vtype, var_ptr->ndims);
      }
     
      fortprintf(fd, "      end if\n\n");

      var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate NetCDF reads of dimension and variable IDs
    */
   fd = fopen("netcdf_read_ids.inc", "w");

   fortprintf(fd, "      nferr = nf_inq_unlimdim(input_obj %% rd_ncid, input_obj %% rdDimIDTime)\n");
   fortprintf(fd, "      nferr = nf_inq_dimlen(input_obj %% rd_ncid, input_obj %% rdDimIDTime, input_obj %% rdLocalTime)\n");
   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
         fortprintf(fd, "      nferr = nf_inq_dimid(input_obj %% rd_ncid, \'%s\', input_obj %% rdDimID%s)\n", dim_ptr->name_in_file, dim_ptr->name_in_file);
         fortprintf(fd, "      nferr = nf_inq_dimlen(input_obj %% rd_ncid, input_obj %% rdDimID%s, input_obj %% rdLocal%s)\n", dim_ptr->name_in_file, dim_ptr->name_in_file);
      }
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   var_ptr = vars;
   while (var_ptr) {
      fortprintf(fd, "      nferr = nf_inq_varid(input_obj %% rd_ncid, \'%s\', input_obj %% rdVarID%s)\n", var_ptr->name_in_file, var_ptr->name_in_file);
      var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate code to return dimension given its name
    */
   fd = fopen("get_dimension_by_name.inc", "w");

   dim_ptr = dims;
   while (dim_ptr->constant_value >= 0 || is_derived_dim(dim_ptr->name_in_code)) dim_ptr = dim_ptr->next;
   fortprintf(fd, "      if (trim(dimname) == \'%s\') then\n", dim_ptr->name_in_code);
   fortprintf(fd, "         dimsize = input_obj %% rdLocal%s\n", dim_ptr->name_in_file);
   dim_ptr = dim_ptr->next;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
         fortprintf(fd, "      else if (trim(dimname) == \'%s\') then\n", dim_ptr->name_in_code);
         fortprintf(fd, "         dimsize = input_obj %% rdLocal%s\n", dim_ptr->name_in_file);
      }
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "      end if\n");

   fclose(fd);
   
   
   /*
    *  Generate code to read 0d, 1d, 2d, 3d time-invariant fields
    */
   for(j=0; j<2; j++) {
      for(i=0; i<=3; i++) {
         if (j == 0) {
            sprintf(fname, "input_field%idinteger.inc", i);
            ivtype = INTEGER;
         }
         else {
            sprintf(fname, "input_field%idreal.inc", i);
            ivtype = REAL;
         }
         fd = fopen(fname, "w");
   
         var_ptr = vars;
         while (var_ptr && (var_ptr->ndims != i || var_ptr->vtype != ivtype || var_ptr->timedim)) var_ptr = var_ptr->next;
         if (var_ptr) {
            fortprintf(fd, "      if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
            fortprintf(fd, "         varID = input_obj %% rdVarID%s\n", var_ptr->name_in_file);
            var_ptr = var_ptr->next;
            while (var_ptr) {
               if (var_ptr->ndims == i && var_ptr->vtype == ivtype && !var_ptr->timedim) {
                  fortprintf(fd, "      else if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
                  fortprintf(fd, "         varID = input_obj %% rdVarID%s\n", var_ptr->name_in_file);
               }
               var_ptr = var_ptr->next;
            }
            fortprintf(fd, "      end if\n");
         }
      
         fclose(fd);
      } 
   } 
   
   
   /*
    *  Generate code to read 0d, 1d, 2d, 3d time-varying real fields
    */
   for(i=0; i<=3; i++) { 
      sprintf(fname, "input_field%idreal_time.inc", i);
      fd = fopen(fname, "w");
   
      var_ptr = vars;
      while (var_ptr && (var_ptr->ndims != i || var_ptr->vtype != REAL || !var_ptr->timedim)) var_ptr = var_ptr->next;
      if (var_ptr) {
         fortprintf(fd, "      if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
         fortprintf(fd, "         varID = input_obj %% rdVarID%s\n", var_ptr->name_in_file);
         var_ptr = var_ptr->next;
         while (var_ptr) {
            if (var_ptr->ndims == i && var_ptr->vtype == REAL && var_ptr->timedim) {
               fortprintf(fd, "      else if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
               fortprintf(fd, "         varID = input_obj %% rdVarID%s\n", var_ptr->name_in_file);
            }
            var_ptr = var_ptr->next;
         }
         fortprintf(fd, "      end if\n");
      }
   
      fclose(fd);
   } 
   
}


void gen_writes(struct variable * vars, struct dimension * dims, struct namelist * namelists)
{
   struct variable * var_ptr;
   struct dimension * dim_ptr;
   struct dimension_list * dimlist_ptr, * lastdim;
   struct dtable * dictionary;
   struct namelist * nl;
   FILE * fd;
   char vtype[5];
   char fname[32];
   char * cp1, * cp2;
   int i, j;
   int ivtype;
   
   
   /*
    *  Generate declarations of IDs belonging in io_output_object
    */
   fd = fopen("io_output_obj_decls.inc", "w");

   fortprintf(fd, "      integer :: wrDimIDTime\n");
   dim_ptr = dims;
   while (dim_ptr) {
      fortprintf(fd, "      integer :: wrDimID%s\n", dim_ptr->name_in_file);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   var_ptr = vars;
   while (var_ptr) {
      fortprintf(fd, "      integer :: wrVarID%s\n", var_ptr->name_in_file);
      var_ptr = var_ptr->next;
   }

   fclose(fd);


   /*
    *  Generate declarations of temporary dimension variables used for arguments
    */
   fd = fopen("output_dim_actual_decls.inc", "w");

   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      integer :: %sGlobal\n", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   fclose(fd);


   /*
    *  Generate initialization of temporary dimension variables used for arguments
    */
   fd = fopen("output_dim_inits.inc", "w");

   dim_ptr = dims;
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, "      %sGlobal = block_ptr %% mesh %% %s\n", dim_ptr->name_in_code, dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   fclose(fd);


   /*
    *  Generate actual dimension argument list
    */
   fd = fopen("output_dim_actual_args.inc", "w");
   dim_ptr = dims;
   if (dim_ptr && dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) {
      fortprintf(fd, "                            %sGlobal", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   while (dim_ptr) {
      if (dim_ptr->constant_value < 0 && !is_derived_dim(dim_ptr->name_in_code)) fortprintf(fd, ", %sGlobal", dim_ptr->name_in_code);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, " &\n");

   fclose(fd);


   /*
    *  Generate NetCDF calls to define dimensions, variables, and global attributes
    */
   fd = fopen("netcdf_def_dims_vars.inc", "w");

   fortprintf(fd, "      nferr = nf_def_dim(output_obj %% wr_ncid, \'Time\', NF_UNLIMITED, output_obj %% wrDimIDTime)\n");
   dim_ptr = dims;
   while (dim_ptr) {
      fortprintf(fd, "      nferr = nf_def_dim(output_obj %% wr_ncid, \'%s\', %s, output_obj %% wrDimID%s)\n", dim_ptr->name_in_file, dim_ptr->name_in_code, dim_ptr->name_in_file);
      dim_ptr = dim_ptr->next;
   }
   fortprintf(fd, "\n");

   var_ptr = vars;
   while (var_ptr) {
      fortprintf(fd, "      if (.false. &\n");
      if (var_ptr->iostreams & RESTART0) fortprintf(fd, "          .or. output_obj %% stream == RESTART &\n");
      if (var_ptr->iostreams & OUTPUT0)  fortprintf(fd, "          .or. output_obj %% stream == OUTPUT &\n");
      fortprintf(fd, "      ) then\n");
      dimlist_ptr = var_ptr->dimlist;
      i = 1;
      while(dimlist_ptr) {
         fortprintf(fd, "      dimlist(%i) = output_obj %% wrDimID%s\n", i++, dimlist_ptr->dim->name_in_file);
         dimlist_ptr = dimlist_ptr->next;
      }
      if (var_ptr->timedim) fortprintf(fd, "      dimlist(%i) = output_obj %% wrDimIDTime\n", i++);
      if (var_ptr->vtype == INTEGER)
         fortprintf(fd, "      nferr = nf_def_var(output_obj %% wr_ncid, \'%s\', NF_INT, %i, dimlist, output_obj %% wrVarID%s)\n", var_ptr->name_in_file, var_ptr->ndims + var_ptr->timedim, var_ptr->name_in_file);
      else if (var_ptr->vtype == REAL)
         fortprintf(fd, "      nferr = nf_def_var(output_obj %% wr_ncid, \'%s\', NF_DOUBLE, %i, dimlist, output_obj %% wrVarID%s)\n", var_ptr->name_in_file, var_ptr->ndims + var_ptr->timedim, var_ptr->name_in_file);

      fortprintf(fd, "      end if\n\n");

      var_ptr = var_ptr->next;
   }

   nl = namelists;
   while (nl) {
      if (nl->vtype == INTEGER)
         fortprintf(fd, "      nferr = nf_put_att_int(output_obj %% wr_ncid, NF_GLOBAL, \'%s\', NF_INT, 1, %s)\n", nl->name, nl->name);
      else if (nl->vtype == REAL) {
         fortprintf(fd, "      if (RKIND == 8) then\n", nl->name);
         fortprintf(fd, "         nferr = nf_put_att_double(output_obj %% wr_ncid, NF_GLOBAL, \'%s\', NF_DOUBLE, 1, %s)\n", nl->name, nl->name);
         fortprintf(fd, "      else if (RKIND == 4) then\n", nl->name);
         fortprintf(fd, "         nferr = nf_put_att_real(output_obj %% wr_ncid, NF_GLOBAL, \'%s\', NF_FLOAT, 1, %s)\n", nl->name, nl->name);
         fortprintf(fd, "      end if\n");
      }
      else if (nl->vtype == CHARACTER)
         fortprintf(fd, "      nferr = nf_put_att_text(output_obj %% wr_ncid, NF_GLOBAL, \'%s\', len_trim(%s), trim(%s))\n", nl->name, nl->name, nl->name);
      else if (nl->vtype == LOGICAL) {
         fortprintf(fd, "      if (%s) then\n", nl->name);
         fortprintf(fd, "         nferr = nf_put_att_text(output_obj %% wr_ncid, NF_GLOBAL, \'%s\', 1, \'T\')\n", nl->name);
         fortprintf(fd, "      else\n");
         fortprintf(fd, "         nferr = nf_put_att_text(output_obj %% wr_ncid, NF_GLOBAL, \'%s\', 1, \'F\')\n", nl->name);
         fortprintf(fd, "      end if\n");
      }
      nl = nl->next;
   }

   fclose(fd);   
   
   
   /*
    *  Generate collect and write code
    */
   fd = fopen("io_output_fields.inc", "w");

   var_ptr = vars;
   while (var_ptr) {
      i = 1;
      dimlist_ptr = var_ptr->dimlist;
      if (var_ptr->vtype == INTEGER) sprintf(vtype, "int"); 
      else if (var_ptr->vtype == REAL) sprintf(vtype, "real"); 

      if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
         if (var_ptr->timedim) {
            fortprintf(fd, "      if ((domain %% blocklist %% time_levs(1) %% state %% %s %% ioinfo %% output .and. output_obj %% stream == OUTPUT) .or. &\n", var_ptr->super_array);
            fortprintf(fd, "          (domain %% blocklist %% time_levs(1) %% state %% %s %% ioinfo %% restart .and. output_obj %% stream == RESTART)) then\n", var_ptr->super_array);
         }
         else {
            fortprintf(fd, "      if ((domain %% blocklist %% mesh %% %s %% ioinfo %% output .and. output_obj %% stream == OUTPUT) .or. &\n", var_ptr->super_array);
            fortprintf(fd, "          (domain %% blocklist %% mesh %% %s %% ioinfo %% restart .and. output_obj %% stream == RESTART)) then\n", var_ptr->super_array);
         }
      }
      else {
         if (var_ptr->timedim) {
            fortprintf(fd, "      if ((domain %% blocklist %% time_levs(1) %% state %% %s %% ioinfo %% output .and. output_obj %% stream == OUTPUT) .or. &\n", var_ptr->name_in_code);
            fortprintf(fd, "          (domain %% blocklist %% time_levs(1) %% state %% %s %% ioinfo %% restart .and. output_obj %% stream == RESTART)) then\n", var_ptr->name_in_code);
         }
         else {
            fortprintf(fd, "      if ((domain %% blocklist %% mesh %% %s %% ioinfo %% output .and. output_obj %% stream == OUTPUT) .or. &\n", var_ptr->name_in_code);
            fortprintf(fd, "          (domain %% blocklist %% mesh %% %s %% ioinfo %% restart .and. output_obj %% stream == RESTART)) then\n", var_ptr->name_in_code);
         }
      }

      if (var_ptr->ndims > 0) {
         while (dimlist_ptr) {
               if (i < var_ptr->ndims) {
                  fortprintf(fd, "      %s%id %% ioinfo %% start(%i) = 1\n", vtype, var_ptr->ndims, i);
                  if (dimlist_ptr->dim->constant_value < 0)
                     fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = domain %% blocklist %% mesh %% %s\n", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
                  else
                     fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = %s\n", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
               }
               else {
                  fortprintf(fd, "      %s%id %% ioinfo %% start(%i) = 1\n", vtype, var_ptr->ndims, i);
                  if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
                     split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
                     fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = n%sGlobal%s\n", vtype, var_ptr->ndims, i, cp1, cp2);
                     free(cp1);
                     free(cp2);
                  }
                  else
                     fortprintf(fd, "      %s%id %% ioinfo %% count(%i) = %sGlobal\n", vtype, var_ptr->ndims, i, dimlist_ptr->dim->name_in_code);
               }
            dimlist_ptr = dimlist_ptr->next;
            i++;
         }
   
         fortprintf(fd, "      allocate(%s%id %% array(", vtype, var_ptr->ndims);
         i = 1;
         dimlist_ptr = var_ptr->dimlist;
   
         if (i < var_ptr->ndims)
            if (dimlist_ptr->dim->constant_value < 0)
               fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
            else
               fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
         else {
            if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
               split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
               fortprintf(fd, "n%sGlobal%s", cp1, cp2);
               free(cp1);
               free(cp2);
            }
            else
               fortprintf(fd, "%sGlobal", dimlist_ptr->dim->name_in_code);
            lastdim = dimlist_ptr;
         }
         dimlist_ptr = dimlist_ptr->next;
         i++;
         while (dimlist_ptr) {
            if (i < var_ptr->ndims)
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
            else {
               if (is_derived_dim(dimlist_ptr->dim->name_in_code)) {
                  split_derived_dim_string(dimlist_ptr->dim->name_in_code, &cp1, &cp2);
                  fortprintf(fd, ", n%sGlobal%s", cp1, cp2);
                  free(cp1);
                  free(cp2);
               }
               else
                  fortprintf(fd, ", %sGlobal", dimlist_ptr->dim->name_in_code);
               lastdim = dimlist_ptr;
            }
            dimlist_ptr = dimlist_ptr->next;
            i++;
         }
         fortprintf(fd, "))\n\n");

         if (strncmp(var_ptr->super_array, "-", 1024) != 0) {
            if (var_ptr->ndims > 0) {
               fortprintf(fd, "      allocate(super_%s%id(", vtype, var_ptr->ndims);
               i = 1;
               dimlist_ptr = var_ptr->dimlist;
               while (dimlist_ptr) {
                  if (dimlist_ptr->dim->constant_value < 0)
                     fortprintf(fd, "domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
                  else
                     fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
   
                  if (i < var_ptr->ndims) fortprintf(fd, ", ");
      
                  dimlist_ptr = dimlist_ptr->next;
                  i++;
               }
               fortprintf(fd, "))\n\n");
            }

            /* Copy from field to super_ array */
            i = 1;
            dimlist_ptr = var_ptr->dimlist;
            while (i <= var_ptr->ndims) {
               if (dimlist_ptr->dim->constant_value < 0)
                  fortprintf(fd, "      do i%i=1,domain %% blocklist %% mesh %% %s\n", i, dimlist_ptr->dim->name_in_code);
               else
                  fortprintf(fd, "      do i%i=1,%s\n", i, dimlist_ptr->dim->name_in_code);

               i++;
               dimlist_ptr = dimlist_ptr->next;
            }

            fortprintf(fd, "         super_%s%id(", vtype, var_ptr->ndims);
            for(i=1; i<=var_ptr->ndims; i++) {
               fortprintf(fd, "i%i",i);
               if (i < var_ptr->ndims) fortprintf(fd, ",");
            }
            if (var_ptr->timedim) 
               fortprintf(fd, ") = domain %% blocklist %% time_levs(1) %% state %% %s %% array(", var_ptr->super_array);
            else
               fortprintf(fd, ") = domain %% blocklist %% mesh %% %s %% array(", var_ptr->super_array);
            fortprintf(fd, "index_%s", var_ptr->name_in_code);
            for(i=1; i<=var_ptr->ndims; i++) {
               fortprintf(fd, ",i%i",i);
            }
            fortprintf(fd, ")\n");

            i = 1;
            while (i <= var_ptr->ndims) {
               fortprintf(fd, "      end do\n");
               i++;
            }
         }

         fortprintf(fd, "      %s%id %% ioinfo %% fieldName = \'%s\'\n", vtype, var_ptr->ndims, var_ptr->name_in_file);
         fortprintf(fd, "      call dmpar_alltoall_field(domain %% dminfo, &\n");
         if (strncmp(var_ptr->super_array, "-", 1024) != 0)
            fortprintf(fd, "                                super_%s%id, %s%id %% array, &\n", vtype, var_ptr->ndims, vtype, var_ptr->ndims);
         else {
            if (var_ptr->timedim) 
               fortprintf(fd, "                                domain %% blocklist %% time_levs(1) %% state %% %s %% array, %s%id %% array, &\n", var_ptr->name_in_code, vtype, var_ptr->ndims);
            else
               fortprintf(fd, "                                domain %% blocklist %% mesh %% %s %% array, %s%id %% array, &\n", var_ptr->name_in_code, vtype, var_ptr->ndims);
         }
   
         i = 1;
         dimlist_ptr = var_ptr->dimlist;
         
         if (dimlist_ptr->dim->constant_value < 0)
            fortprintf(fd, "                                domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
         else
            fortprintf(fd, "                                %s", dimlist_ptr->dim->name_in_code);
    
         dimlist_ptr = dimlist_ptr->next;
         i++;
         while (dimlist_ptr) {
            if (dimlist_ptr->dim->constant_value < 0)
               fortprintf(fd, ", domain %% blocklist %% mesh %% %s", dimlist_ptr->dim->name_in_code);
            else
               fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
   
            dimlist_ptr = dimlist_ptr->next;
            i++;
         }     
   
         if (is_derived_dim(lastdim->dim->name_in_code)) {
            split_derived_dim_string(lastdim->dim->name_in_code, &cp1, &cp2);
            fortprintf(fd, ", n%sGlobal%s, &\n", cp1, cp2);
            fortprintf(fd, "                                output_obj %% send%sList, output_obj %% recv%sList)\n", cp1, cp1);
            free(cp1);
            free(cp2);
         }
         else {
            fortprintf(fd, ", %sGlobal, &\n", lastdim->dim->name_in_code);
            fortprintf(fd, "                                output_obj %% send%sList, output_obj %% recv%sList)\n", lastdim->dim->name_in_code+1, lastdim->dim->name_in_code+1);
         }
      }
      else {
         fortprintf(fd, "      %s%id %% ioinfo %% fieldName = \'%s\'\n", vtype, var_ptr->ndims, var_ptr->name_in_file);
         if (var_ptr->timedim) 
            fortprintf(fd, "      %s%id %% scalar = domain %% blocklist %% time_levs(1) %% state %% %s %% scalar\n", vtype, var_ptr->ndims, var_ptr->name_in_code);
         else
            fortprintf(fd, "      %s%id %% scalar = domain %% blocklist %% mesh %% %s %% scalar\n", vtype, var_ptr->ndims, var_ptr->name_in_code);
      }

      if (var_ptr->timedim)
         fortprintf(fd, "      if (domain %% dminfo %% my_proc_id == IO_NODE) call io_output_field_time(output_obj, %s%id)\n", vtype, var_ptr->ndims);
      else
         fortprintf(fd, "      if (domain %% dminfo %% my_proc_id == IO_NODE) call io_output_field(output_obj, %s%id)\n", vtype, var_ptr->ndims);
      if (var_ptr->ndims > 0) {
         fortprintf(fd, "      deallocate(%s%id %% array)\n", vtype, var_ptr->ndims);
         if (strncmp(var_ptr->super_array, "-", 1024) != 0)
            fortprintf(fd, "      deallocate(super_%s%id)\n", vtype, var_ptr->ndims);
      }
      fortprintf(fd, "      end if\n\n");

      var_ptr = var_ptr->next;
   }

   fclose(fd);
   
   
   /*
    *  Generate code to write 0d, 1d, 2d, 3d time-invariant fields
    */
   for(j=0; j<2; j++) {
      for(i=0; i<=3; i++) {
         if (j == 0) {
            sprintf(fname, "output_field%idinteger.inc", i);
            ivtype = INTEGER;
         }
         else {
            sprintf(fname, "output_field%idreal.inc", i);
            ivtype = REAL;
         }
         fd = fopen(fname, "w");
   
         var_ptr = vars;
         while (var_ptr && (var_ptr->ndims != i || var_ptr->vtype != ivtype || var_ptr->timedim)) var_ptr = var_ptr->next;
         if (var_ptr) {
            fortprintf(fd, "      if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
            fortprintf(fd, "         varID = output_obj %% wrVarID%s\n", var_ptr->name_in_file);
            var_ptr = var_ptr->next;
            while (var_ptr) {
               if (var_ptr->ndims == i && var_ptr->vtype == ivtype && !var_ptr->timedim) {
                  fortprintf(fd, "      else if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
                  fortprintf(fd, "         varID = output_obj %% wrVarID%s\n", var_ptr->name_in_file);
               }
               var_ptr = var_ptr->next;
            }
            fortprintf(fd, "      end if\n");
         }
      
         fclose(fd);
      } 
   } 

   
   /*
    *  Generate code to write 0d, 1d, 2d, 3d real time-varying fields
    */
   for(i=0; i<=3; i++) {
      sprintf(fname, "output_field%idreal_time.inc", i);
      fd = fopen(fname, "w");

      var_ptr = vars;
      while (var_ptr && (var_ptr->ndims != i || var_ptr->vtype != REAL || !var_ptr->timedim)) var_ptr = var_ptr->next;
      if (var_ptr) {
         fortprintf(fd, "      if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
         fortprintf(fd, "         varID = output_obj %% wrVarID%s\n", var_ptr->name_in_file);
         var_ptr = var_ptr->next;
         while (var_ptr) {
            if (var_ptr->ndims == i && var_ptr->vtype == REAL && var_ptr->timedim) {
               fortprintf(fd, "      else if (trim(field %% ioinfo %% fieldName) == \'%s\') then\n", var_ptr->name_in_file);
               fortprintf(fd, "         varID = output_obj %% wrVarID%s\n", var_ptr->name_in_file);
            }
            var_ptr = var_ptr->next;
         }
         fortprintf(fd, "      end if\n");
      }
   
      fclose(fd);
   }
   
}
-------------- next part --------------
module io_input

   use grid_types
   use dmpar
   use block_decomp
   use sort
   use configure

#ifdef HAVE_ZOLTAN
   use zoltan_interface
#endif

   type io_input_object
      character (len=1024) :: filename
      integer :: rd_ncid

      integer :: time

#include "io_input_obj_decls.inc"
   end type io_input_object


   interface io_input_field
      module procedure io_input_field0dReal
      module procedure io_input_field1dReal
      module procedure io_input_field2dReal
      module procedure io_input_field3dReal
      module procedure io_input_field1dInteger
      module procedure io_input_field2dInteger
   end interface io_input_field

   interface io_input_field_time
      module procedure io_input_field0dReal_time
      module procedure io_input_field1dReal_time
      module procedure io_input_field2dReal_time
      module procedure io_input_field3dReal_time
   end interface io_input_field_time
 

   contains


   subroutine input_state_for_domain(domain)
   
      implicit none
   
      type (domain_type), pointer :: domain
   
      integer :: i, j, k
      type (io_input_object) :: input_obj
#include "dim_decls.inc"
   
      integer :: readCellStart, readCellEnd, nReadCells
      integer :: readEdgeStart, readEdgeEnd, nReadEdges
      integer :: readVertexStart, readVertexEnd, nReadVertices
      integer :: readVertLevelStart, readVertLevelEnd, nReadVertLevels
   
      type (field1dInteger) :: indexToCellIDField
      type (field1dInteger) :: indexToEdgeIDField
      type (field1dInteger) :: indexToVertexIDField
      type (field1dInteger) :: nEdgesOnCellField
      type (field2dInteger) :: cellsOnCellField
      type (field2dInteger) :: edgesOnCellField
      type (field2dInteger) :: verticesOnCellField
      type (field2dInteger) :: cellsOnEdgeField
      type (field2dInteger) :: cellsOnVertexField

#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      type (field1dReal) :: xCellField,   yCellField,   zCellField
      type (field1dReal) :: xEdgeField,   yEdgeField,   zEdgeField
      type (field1dReal) :: xVertexField, yVertexField, zVertexField
#endif
#endif

      type (field1dReal) :: xtime
   
      integer, dimension(:), pointer :: indexToCellID_0Halo
      integer, dimension(:), pointer :: nEdgesOnCell_0Halo
      integer, dimension(:,:), pointer :: cellsOnCell_0Halo
   
      integer, dimension(:,:), pointer :: edgesOnCell_2Halo
      integer, dimension(:,:), pointer :: verticesOnCell_2Halo
      integer, dimension(:,:), pointer :: cellsOnEdge_2Halo
      integer, dimension(:,:), pointer :: cellsOnVertex_2Halo

      integer, dimension(:,:), pointer :: cellIDSorted
      integer, dimension(:,:), pointer :: edgeIDSorted
      integer, dimension(:,:), pointer :: vertexIDSorted

#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      real (kind=RKIND), dimension(:), pointer :: xCell,   yCell,   zCell
      real (kind=RKIND), dimension(:), pointer :: xEdge,   yEdge,   zEdge
      real (kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
#endif
#endif
   
      integer, dimension(:), pointer :: local_cell_list, local_edge_list, local_vertex_list
      integer, dimension(:), pointer :: local_vertlevel_list, needed_vertlevel_list
      integer :: nlocal_edges, nlocal_vertices
      type (exchange_list), pointer :: sendCellList, recvCellList
      type (exchange_list), pointer :: sendEdgeList, recvEdgeList
      type (exchange_list), pointer :: sendVertexList, recvVertexList
      type (exchange_list), pointer :: sendVertLevelList, recvVertLevelList
      type (exchange_list), pointer :: send1Halo, recv1Halo
      type (exchange_list), pointer :: send2Halo, recv2Halo
      type (graph) :: partial_global_graph_info
      type (graph) :: block_graph_0Halo, block_graph_1Halo, block_graph_2Halo
      integer :: ghostEdgeStart, ghostVertexStart

      real (kind=RKIND) :: tdiff
      integer :: idxTdiff

      if (config_do_restart) then
         input_obj % filename = trim(config_restart_name)
      else
         input_obj % filename = trim(config_input_name)
      end if
      call io_input_init(input_obj, domain % dminfo)
   

      !
      ! Read global number of cells/edges/vertices
      !
#include "read_dims.inc"
   
      !
      ! Determine the range of cells/edges/vertices that a processor will initially read
      !   from the input file
      !
      call dmpar_get_index_range(domain % dminfo, 1, nCells, readCellStart, readCellEnd)   
      nReadCells    = readCellEnd - readCellStart + 1
   
      call dmpar_get_index_range(domain % dminfo, 1, nEdges, readEdgeStart, readEdgeEnd)   
      nReadEdges    = readEdgeEnd - readEdgeStart + 1
   
      call dmpar_get_index_range(domain % dminfo, 1, nVertices, readVertexStart, readVertexEnd)   
      nReadVertices = readVertexEnd - readVertexStart + 1

      readVertLevelStart = 1
      readVertLevelEnd = nVertLevels
      nReadVertLevels = nVertLevels
   
   
      !
      ! Allocate and read fields that we will need in order to ultimately work out
      !   which cells/edges/vertices are owned by each block, and which are ghost
      !

      ! Global cell indices
      allocate(indexToCellIDField % ioinfo)
      indexToCellIDField % ioinfo % fieldName = 'indexToCellID'
      indexToCellIDField % ioinfo % start(1) = readCellStart
      indexToCellIDField % ioinfo % count(1) = nReadCells
      allocate(indexToCellIDField % array(nReadCells))
      call io_input_field(input_obj, indexToCellIDField)
   
#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      ! Cell x-coordinates (in 3d Cartesian space)
      allocate(xCellField % ioinfo)
      xCellField % ioinfo % fieldName = 'xCell'
      xCellField % ioinfo % start(1) = readCellStart
      xCellField % ioinfo % count(1) = nReadCells
      allocate(xCellField % array(nReadCells))
      call io_input_field(input_obj, xCellField)

      ! Cell y-coordinates (in 3d Cartesian space)
      allocate(yCellField % ioinfo)
      yCellField % ioinfo % fieldName = 'yCell'
      yCellField % ioinfo % start(1) = readCellStart
      yCellField % ioinfo % count(1) = nReadCells
      allocate(yCellField % array(nReadCells))
      call io_input_field(input_obj, yCellField)

      ! Cell z-coordinates (in 3d Cartesian space)
      allocate(zCellField % ioinfo)
      zCellField % ioinfo % fieldName = 'zCell'
      zCellField % ioinfo % start(1) = readCellStart
      zCellField % ioinfo % count(1) = nReadCells
      allocate(zCellField % array(nReadCells))
      call io_input_field(input_obj, zCellField)
#endif
#endif


      ! Global edge indices
      allocate(indexToEdgeIDField % ioinfo)
      indexToEdgeIDField % ioinfo % fieldName = 'indexToEdgeID'
      indexToEdgeIDField % ioinfo % start(1) = readEdgeStart
      indexToEdgeIDField % ioinfo % count(1) = nReadEdges
      allocate(indexToEdgeIDField % array(nReadEdges))
      call io_input_field(input_obj, indexToEdgeIDField)
   
#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      ! Edge x-coordinates (in 3d Cartesian space)
      allocate(xEdgeField % ioinfo)
      xEdgeField % ioinfo % fieldName = 'xEdge'
      xEdgeField % ioinfo % start(1) = readEdgeStart
      xEdgeField % ioinfo % count(1) = nReadEdges
      allocate(xEdgeField % array(nReadEdges))
      call io_input_field(input_obj, xEdgeField)

      ! Edge y-coordinates (in 3d Cartesian space)
      allocate(yEdgeField % ioinfo)
      yEdgeField % ioinfo % fieldName = 'yEdge'
      yEdgeField % ioinfo % start(1) = readEdgeStart
      yEdgeField % ioinfo % count(1) = nReadEdges
      allocate(yEdgeField % array(nReadEdges))
      call io_input_field(input_obj, yEdgeField)

      ! Edge z-coordinates (in 3d Cartesian space)
      allocate(zEdgeField % ioinfo)
      zEdgeField % ioinfo % fieldName = 'zEdge'
      zEdgeField % ioinfo % start(1) = readEdgeStart
      zEdgeField % ioinfo % count(1) = nReadEdges
      allocate(zEdgeField % array(nReadEdges))
      call io_input_field(input_obj, zEdgeField)
#endif
#endif

      ! Global vertex indices
      allocate(indexToVertexIDField % ioinfo)
      indexToVertexIDField % ioinfo % fieldName = 'indexToVertexID'
      indexToVertexIDField % ioinfo % start(1) = readVertexStart
      indexToVertexIDField % ioinfo % count(1) = nReadVertices
      allocate(indexToVertexIDField % array(nReadVertices))
      call io_input_field(input_obj, indexToVertexIDField)
   
#ifdef HAVE_ZOLTAN
#ifdef _MPI
      ! Vertex x-coordinates (in 3d Cartesian space)
      allocate(xVertexField % ioinfo)
      xVertexField % ioinfo % fieldName = 'xVertex'
      xVertexField % ioinfo % start(1) = readVertexStart
      xVertexField % ioinfo % count(1) = nReadVertices
      allocate(xVertexField % array(nReadVertices))
      call io_input_field(input_obj, xVertexField)

      ! Vertex y-coordinates (in 3d Cartesian space)
      allocate(yVertexField % ioinfo)
      yVertexField % ioinfo % fieldName = 'yVertex'
      yVertexField % ioinfo % start(1) = readVertexStart
      yVertexField % ioinfo % count(1) = nReadVertices
      allocate(yVertexField % array(nReadVertices))
      call io_input_field(input_obj, yVertexField)

      ! Vertex z-coordinates (in 3d Cartesian space)
      allocate(zVertexField % ioinfo)
      zVertexField % ioinfo % fieldName = 'zVertex'
      zVertexField % ioinfo % start(1) = readVertexStart
      zVertexField % ioinfo % count(1) = nReadVertices
      allocate(zVertexField % array(nReadVertices))
      call io_input_field(input_obj, zVertexField)
#endif
#endif

      ! Number of cell/edges/vertices adjacent to each cell
      allocate(nEdgesOnCellField % ioinfo)
      nEdgesOnCellField % ioinfo % fieldName = 'nEdgesOnCell'
      nEdgesOnCellField % ioinfo % start(1) = readCellStart
      nEdgesOnCellField % ioinfo % count(1) = nReadCells
      allocate(nEdgesOnCellField % array(nReadCells))
      call io_input_field(input_obj, nEdgesOnCellField)
   
      ! Global indices of cells adjacent to each cell
      allocate(cellsOnCellField % ioinfo)
      cellsOnCellField % ioinfo % fieldName = 'cellsOnCell'
      cellsOnCellField % ioinfo % start(1) = 1
      cellsOnCellField % ioinfo % start(2) = readCellStart
      cellsOnCellField % ioinfo % count(1) = maxEdges
      cellsOnCellField % ioinfo % count(2) = nReadCells
      allocate(cellsOnCellField % array(maxEdges,nReadCells))
      call io_input_field(input_obj, cellsOnCellField)
   
      ! Global indices of edges adjacent to each cell
      allocate(edgesOnCellField % ioinfo)
      edgesOnCellField % ioinfo % fieldName = 'edgesOnCell'
      edgesOnCellField % ioinfo % start(1) = 1
      edgesOnCellField % ioinfo % start(2) = readCellStart
      edgesOnCellField % ioinfo % count(1) = maxEdges
      edgesOnCellField % ioinfo % count(2) = nReadCells
      allocate(edgesOnCellField % array(maxEdges,nReadCells))
      call io_input_field(input_obj, edgesOnCellField)
   
      ! Global indices of vertices adjacent to each cell
      allocate(verticesOnCellField % ioinfo)
      verticesOnCellField % ioinfo % fieldName = 'verticesOnCell'
      verticesOnCellField % ioinfo % start(1) = 1
      verticesOnCellField % ioinfo % start(2) = readCellStart
      verticesOnCellField % ioinfo % count(1) = maxEdges
      verticesOnCellField % ioinfo % count(2) = nReadCells
      allocate(verticesOnCellField % array(maxEdges,nReadCells))
      call io_input_field(input_obj, verticesOnCellField)
   
      ! Global indices of cells adjacent to each edge
      !    used for determining which edges are owned by a block, where 
      !    iEdge is owned iff cellsOnEdge(1,iEdge) is an owned cell
      allocate(cellsOnEdgeField % ioinfo)
      cellsOnEdgeField % ioinfo % fieldName = 'cellsOnEdge'
      cellsOnEdgeField % ioinfo % start(1) = 1
      cellsOnEdgeField % ioinfo % start(2) = readEdgeStart
      cellsOnEdgeField % ioinfo % count(1) = 2
      cellsOnEdgeField % ioinfo % count(2) = nReadEdges
      allocate(cellsOnEdgeField % array(2,nReadEdges))
      call io_input_field(input_obj, cellsOnEdgeField)
   
      ! Global indices of cells adjacent to each vertex
      !    used for determining which vertices are owned by a block, where 
      !    iVtx is owned iff cellsOnVertex(1,iVtx) is an owned cell
      allocate(cellsOnVertexField % ioinfo)
      cellsOnVertexField % ioinfo % fieldName = 'cellsOnVertex'
      cellsOnVertexField % ioinfo % start(1) = 1
      cellsOnVertexField % ioinfo % start(2) = readVertexStart
      cellsOnVertexField % ioinfo % count(1) = vertexDegree
      cellsOnVertexField % ioinfo % count(2) = nReadVertices
      allocate(cellsOnVertexField % array(vertexDegree,nReadVertices))
      call io_input_field(input_obj, cellsOnVertexField)
   
   
      !
      ! Set up a graph derived data type describing the connectivity for the cells 
      !   that were read by this process
      ! A partial description is passed to the block decomp module by each process,
      !   and the block decomp module returns with a list of global cell indices
      !   that belong to the block on this process
      !
      partial_global_graph_info % nVertices = nReadCells
      partial_global_graph_info % nVerticesTotal = nCells
      partial_global_graph_info % maxDegree = maxEdges
      partial_global_graph_info % ghostStart = nVertices+1
      allocate(partial_global_graph_info % vertexID(nReadCells))
      allocate(partial_global_graph_info % nAdjacent(nReadCells))
      allocate(partial_global_graph_info % adjacencyList(maxEdges, nReadCells))
   
      partial_global_graph_info % vertexID(:) = indexToCellIDField % array(:)
      partial_global_graph_info % nAdjacent(:) = nEdgesOnCellField % array(:)
      partial_global_graph_info % adjacencyList(:,:) = cellsOnCellField % array(:,:)
      
   
      ! TODO: Ensure (by renaming or exchanging) that initial cell range on each proc is contiguous
      !       This situation may occur when reading a restart file with cells/edges/vertices written
      !       in a scrambled order
   

      ! Determine which cells are owned by this process
      call block_decomp_cells_for_proc(domain % dminfo, partial_global_graph_info, local_cell_list)

      deallocate(partial_global_graph_info % vertexID)
      deallocate(partial_global_graph_info % nAdjacent)
      deallocate(partial_global_graph_info % adjacencyList)
   
   
      allocate(indexToCellID_0Halo(size(local_cell_list)))
      allocate(nEdgesOnCell_0Halo(size(local_cell_list)))
      allocate(cellsOnCell_0Halo(maxEdges, size(local_cell_list)))
   
#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      allocate(xCell(size(local_cell_list)))
      allocate(yCell(size(local_cell_list)))
      allocate(zCell(size(local_cell_list)))
#endif
#endif
   
      !
      ! Now that each process has a list of cells that it owns, exchange cell connectivity 
      !   information between the processes that read info for a cell and those that own that cell
      !
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToCellIDField % array), size(local_cell_list), &
                                indexToCellIDField % array, local_cell_list, &
                                sendCellList, recvCellList)
   
      call dmpar_alltoall_field(domain % dminfo, indexToCellIDField % array, indexToCellID_0Halo, &
                                size(indexToCellIDField % array), size(local_cell_list), &
                                sendCellList, recvCellList)
   
      call dmpar_alltoall_field(domain % dminfo, nEdgesOnCellField % array, nEdgesOnCell_0Halo, &
                                size(indexToCellIDField % array), size(local_cell_list), &
                                sendCellList, recvCellList)
   
      call dmpar_alltoall_field(domain % dminfo, cellsOnCellField % array, cellsOnCell_0Halo, &
                                size(cellsOnCellField % array, 1), size(indexToCellIDField % array), size(local_cell_list), &
                                sendCellList, recvCellList)
   
#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      call dmpar_alltoall_field(domain % dminfo, xCellField % array, xCell, &
                                size(xCellField % array), size(local_cell_list), &
                                sendCellList, recvCellList)
   
      call dmpar_alltoall_field(domain % dminfo, yCellField % array, yCell, &
                                size(yCellField % array), size(local_cell_list), &
                                sendCellList, recvCellList)
   
      call dmpar_alltoall_field(domain % dminfo, zCellField % array, zCell, &
                                size(zCellField % array), size(local_cell_list), &
                                sendCellList, recvCellList)
#endif
#endif


      deallocate(sendCellList % list)
      deallocate(sendCellList)
      deallocate(recvCellList % list)
      deallocate(recvCellList)



      !
      ! Build a graph of cell connectivity based on cells owned by this process
      !
      block_graph_0Halo % nVerticesTotal = size(local_cell_list)
      block_graph_0Halo % nVertices = size(local_cell_list)
      block_graph_0Halo % maxDegree = maxEdges
      block_graph_0Halo % ghostStart = size(local_cell_list) + 1
      allocate(block_graph_0Halo % vertexID(size(local_cell_list)))
      allocate(block_graph_0Halo % nAdjacent(size(local_cell_list)))
      allocate(block_graph_0Halo % adjacencyList(maxEdges, size(local_cell_list)))
   
      block_graph_0Halo % vertexID(:) = indexToCellID_0Halo(:)
      block_graph_0Halo % nAdjacent(:) = nEdgesOnCell_0Halo(:)
      block_graph_0Halo % adjacencyList(:,:) = cellsOnCell_0Halo(:,:)
   
      ! Get back a graph describing the owned cells plus the cells in the 1-halo
      call block_decomp_add_halo(domain % dminfo, block_graph_0Halo, block_graph_1Halo)
   
   
      !
      ! Work out exchange lists for 1-halo and exchange cell information for 1-halo
      !
      call dmpar_get_owner_list(domain % dminfo, &
                                block_graph_0Halo % nVertices, block_graph_1Halo % nVerticesTotal, &
                                block_graph_0Halo % vertexID,  block_graph_1Halo % vertexID, &
                                send1Halo, recv1Halo)
   
      call dmpar_alltoall_field(domain % dminfo, block_graph_0Halo % vertexID, block_graph_1Halo % vertexID, &
                                block_graph_0Halo % nVertices, block_graph_1Halo % nVerticesTotal, &
                                send1Halo, recv1Halo)
   
      call dmpar_alltoall_field(domain % dminfo, block_graph_0Halo % nAdjacent, block_graph_1Halo % nAdjacent, &
                                block_graph_0Halo % nVertices, block_graph_1Halo % nVerticesTotal, &
                                send1Halo, recv1Halo)
   
      call dmpar_alltoall_field(domain % dminfo, block_graph_0Halo % adjacencyList, block_graph_1Halo % adjacencyList, &
                                block_graph_0Halo % maxDegree, block_graph_0Halo % nVertices, block_graph_1Halo % nVerticesTotal, &
                                send1Halo, recv1Halo)
   
   
      !
      ! Work out exchange lists for 2-halo and exchange cell information for 2-halo
      !
      block_graph_1Halo % nVertices = block_graph_1Halo % nVerticesTotal
      block_graph_1Halo % ghostStart = block_graph_1Halo % nVerticesTotal + 1
     
      ! Get back a graph describing the owned and 1-halo cells plus the cells in the 2-halo
      call block_decomp_add_halo(domain % dminfo, block_graph_1Halo, block_graph_2Halo)
   
      block_graph_2Halo % nVertices = block_graph_0Halo % nVertices
      block_graph_2Halo % ghostStart = block_graph_2Halo % nVertices + 1

#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      !! For now, only use Zoltan with MPI
      !! Zoltan initialization
      call zoltanStart()

      !! Zoltan hook for cells
      call zoltanOrderLocHSFC_Cells(block_graph_2Halo%nVertices,block_graph_2Halo%VertexID,3,xCell,yCell,zCell)
#endif
#endif

      call dmpar_get_owner_list(domain % dminfo, &
                                block_graph_0Halo % nVertices, block_graph_2Halo % nVerticesTotal, &
                                block_graph_0Halo % vertexID,  block_graph_2Halo % vertexID, &
                                send2Halo, recv2Halo)
   
      call dmpar_alltoall_field(domain % dminfo, block_graph_0Halo % vertexID, block_graph_2Halo % vertexID, &
                                block_graph_0Halo % nVertices, block_graph_2Halo % nVerticesTotal, &
                                send2Halo, recv2Halo)
   
      call dmpar_alltoall_field(domain % dminfo, block_graph_0Halo % nAdjacent, block_graph_2Halo % nAdjacent, &
                                block_graph_0Halo % nVertices, block_graph_2Halo % nVerticesTotal, &
                                send2Halo, recv2Halo)
   
      call dmpar_alltoall_field(domain % dminfo, block_graph_0Halo % adjacencyList, block_graph_2Halo % adjacencyList, &
                                block_graph_0Halo % maxDegree, block_graph_0Halo % nVertices, block_graph_2Halo % nVerticesTotal, &
                                send2Halo, recv2Halo)


   
      !
      ! Knowing which cells are in block and the 2-halo, we can exchange lists of which edges are
      !   on each cell and which vertices are on each cell from the processes that read these
      !   fields for each cell to the processes that own the cells
      !
      allocate(edgesOnCell_2Halo(maxEdges, block_graph_2Halo % nVerticesTotal))
      allocate(verticesOnCell_2Halo(maxEdges, block_graph_2Halo % nVerticesTotal))
   
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToCellIDField % array), block_graph_2Halo % nVerticesTotal, &
                                indexToCellIDField % array, block_graph_2Halo % vertexID, &
                                sendCellList, recvCellList)
   
      call dmpar_alltoall_field(domain % dminfo, edgesOnCellField % array, edgesOnCell_2Halo, &
                                maxEdges, nReadCells, block_graph_2Halo % nVerticesTotal, &
                                sendCellList, recvCellList)
   
      call dmpar_alltoall_field(domain % dminfo, verticesOnCellField % array, verticesOnCell_2Halo, &
                                maxEdges, nReadCells, block_graph_2Halo % nVerticesTotal, &
                                sendCellList, recvCellList)

   
      ! 
      ! Get a list of which edges and vertices are adjacent to cells (including halo cells) in block
      ! 
      call block_decomp_all_edges_in_block(maxEdges, block_graph_2Halo % nVerticesTotal, block_graph_2Halo % nAdjacent, &
                                           edgesOnCell_2Halo, nlocal_edges, local_edge_list)
      call block_decomp_all_edges_in_block(maxEdges, block_graph_2Halo % nVerticesTotal, block_graph_2Halo % nAdjacent, &
                                           verticesOnCell_2Halo, nlocal_vertices, local_vertex_list)
   
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToEdgeIDField % array), nlocal_edges, &
                                indexToEdgeIDField % array, local_edge_list, &
                                sendEdgeList, recvEdgeList)
   
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToVertexIDField % array), nlocal_vertices, &
                                indexToVertexIDField % array, local_vertex_list, &
                                sendVertexList, recvVertexList)
   
   
   
      ! 
      ! Work out which edges and vertices are owned by this process, and which are ghost
      ! 
      allocate(cellsOnEdge_2Halo(2,nlocal_edges))
      allocate(cellsOnVertex_2Halo(vertexDegree,nlocal_vertices))
   
      call dmpar_alltoall_field(domain % dminfo, cellsOnEdgeField % array, cellsOnEdge_2Halo, &
                                2, size(cellsOnEdgeField % array, 2), nlocal_edges, &
                                sendEdgeList, recvEdgeList)
   
      call dmpar_alltoall_field(domain % dminfo, cellsOnVertexField % array, cellsOnVertex_2Halo, &
                                vertexDegree, size(cellsOnVertexField % array, 2), nlocal_vertices, &
                                sendVertexList, recvVertexList)
   
   
      call block_decomp_partitioned_edge_list(block_graph_2Halo % nVertices, &
                                              block_graph_2Halo % vertexID(1:block_graph_2Halo % nVertices), &
                                              2, nlocal_edges, cellsOnEdge_2Halo, local_edge_list, ghostEdgeStart)
      call block_decomp_partitioned_edge_list(block_graph_2Halo % nVertices, &
                                              block_graph_2Halo % vertexID(1:block_graph_2Halo % nVertices), &
                                              vertexDegree, nlocal_vertices, cellsOnVertex_2Halo, local_vertex_list, ghostVertexStart)


      ! At this point, local_edge_list(1;ghostEdgeStart-1) contains all of the owned edges for this block
      !   and local_edge_list(ghostEdgeStart:nlocal_edges) contains all of the ghost edges

      ! At this point, local_vertex_list(1;ghostVertexStart-1) contains all of the owned vertices for this block
      !   and local_vertex_list(ghostVertexStart:nlocal_vertices) contains all of the ghost vertices

      ! Also, at this point, block_graph_2Halo % vertexID(1:block_graph_2Halo%nVertices) contains all of the owned
      !   cells for this block, and block_graph_2Halo % vertexID(block_graph_2Halo%nVertices+1:block_graph_2Halo%nVerticesTotal)
      !   contains all of the ghost cells


      deallocate(sendEdgeList % list)
      deallocate(sendEdgeList)
      deallocate(recvEdgeList % list)
      deallocate(recvEdgeList)
   
      deallocate(sendVertexList % list)
      deallocate(sendVertexList)
      deallocate(recvVertexList % list)
      deallocate(recvVertexList)
   
#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      allocate(xEdge(nlocal_edges))
      allocate(yEdge(nlocal_edges))
      allocate(zEdge(nlocal_edges))
      allocate(xVertex(nlocal_vertices))
      allocate(yVertex(nlocal_vertices))
      allocate(zVertex(nlocal_vertices))
#endif
#endif
    
      !
      ! Knowing which edges/vertices are owned by this block and which are actually read
      !   from the input or restart file, we can build exchange lists to perform 
      !   all-to-all field exchanges from process that reads a field to the processes that
      !   need them
      !
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToEdgeIDField % array), nlocal_edges, &
                                indexToEdgeIDField % array, local_edge_list, &
                                sendEdgeList, recvEdgeList)
   
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToVertexIDField % array), nlocal_vertices, &
                                indexToVertexIDField % array, local_vertex_list, &
                                sendVertexList, recvVertexList)

#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      call dmpar_alltoall_field(domain % dminfo, xEdgeField % array, xEdge, &
                                size(xEdgeField % array), nlocal_edges, &
                                sendEdgeList, recvEdgeList)
      call dmpar_alltoall_field(domain % dminfo, yEdgeField % array, yEdge, &
                                size(yEdgeField % array), nlocal_edges, &
                                sendEdgeList, recvEdgeList)
      call dmpar_alltoall_field(domain % dminfo, zEdgeField % array, zEdge, &
                                size(zEdgeField % array), nlocal_edges, &
                                sendEdgeList, recvEdgeList)

      call dmpar_alltoall_field(domain % dminfo, xVertexField % array, xVertex, &
                                size(xVertexField % array), nlocal_vertices, &
                                sendVertexList, recvVertexList)
      call dmpar_alltoall_field(domain % dminfo, yVertexField % array, yVertex, &
                                size(yVertexField % array), nlocal_vertices, &
                                sendVertexList, recvVertexList)
      call dmpar_alltoall_field(domain % dminfo, zVertexField % array, zVertex, &
                                size(zVertexField % array), nlocal_vertices, &
                                sendVertexList, recvVertexList)
      !!!!!!!!!!!!!!!!!!
      !! Reorder edges
      !!!!!!!!!!!!!!!!!!
      call zoltanOrderLocHSFC_Edges(ghostEdgeStart-1,local_edge_list,3,xEdge,yEdge,zEdge)
      !!!!!!!!!!!!!!!!!!

      !!!!!!!!!!!!!!!!!!
      !! Reorder vertices
      !!!!!!!!!!!!!!!!!!
      call zoltanOrderLocHSFC_Verts(ghostVertexStart-1,local_vertex_list,3,xVertex,yVertex,zVertex)
      !!!!!!!!!!!!!!!!!!

      deallocate(sendEdgeList % list)
      deallocate(sendEdgeList)
      deallocate(recvEdgeList % list)
      deallocate(recvEdgeList)
   
      deallocate(sendVertexList % list)
      deallocate(sendVertexList)
      deallocate(recvVertexList % list)
      deallocate(recvVertexList)
    
      !
      ! Knowing which edges/vertices are owned by this block and which are actually read
      !   from the input or restart file, we can build exchange lists to perform 
      !   all-to-all field exchanges from process that reads a field to the processes that
      !   need them
      !
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToEdgeIDField % array), nlocal_edges, &
                                indexToEdgeIDField % array, local_edge_list, &
                                sendEdgeList, recvEdgeList)
   
      call dmpar_get_owner_list(domain % dminfo, &
                                size(indexToVertexIDField % array), nlocal_vertices, &
                                indexToVertexIDField % array, local_vertex_list, &
                                sendVertexList, recvVertexList)

#endif
#endif

      ! 
      ! Build ownership and exchange lists for vertical levels
      ! Essentially, process 0 owns all vertical levels when reading and writing,
      ! and it distributes them or gathers them to/from all other processes
      ! 
      if (domain % dminfo % my_proc_id == 0) then
         allocate(local_vertlevel_list(nVertLevels))
         do i=1,nVertLevels
            local_vertlevel_list(i) = i
         end do
      else
         allocate(local_vertlevel_list(0))
      end if
      allocate(needed_vertlevel_list(nVertLevels))
      do i=1,nVertLevels
         needed_vertlevel_list(i) = i
      end do

      call dmpar_get_owner_list(domain % dminfo, &
                                size(local_vertlevel_list), size(needed_vertlevel_list), &
                                local_vertlevel_list, needed_vertlevel_list, &
                                sendVertLevelList, recvVertLevelList)

      deallocate(local_vertlevel_list)
      deallocate(needed_vertlevel_list)


      !
      ! Read and distribute all fields given ownership lists and exchange lists (maybe already in block?)
      !
      allocate(domain % blocklist)

      nCells = block_graph_2Halo % nVerticesTotal
      nEdges = nlocal_edges
      nVertices = nlocal_vertices

      call allocate_block(domain % blocklist, domain, &
#include "dim_dummy_args.inc"
                         )

      if (.not. config_do_restart) then
         input_obj % time = 1
      else
         input_obj % time = 1

         !
         ! If doing a restart, we need to decide which time slice to read from the 
         !   restart file
         !
         if (input_obj % rdLocalTime <= 0) then
            write(0,*) 'Error: Couldn''t find any times in restart file.'
            call dmpar_abort(domain % dminfo)
         end if
         if (domain % dminfo % my_proc_id == IO_NODE) then
            allocate(xtime % ioinfo)
            xtime % ioinfo % start(1) = 1
            xtime % ioinfo % count(1) = input_obj % rdLocalTime
            allocate(xtime % array(input_obj % rdLocalTime))

            xtime % ioinfo % fieldName = 'xtime'
            call io_input_field(input_obj, xtime)

            tdiff = 1.E20
            do i=1,input_obj % rdLocalTime
               if (abs(xtime % array(i) - config_restart_time) < tdiff) then
                  input_obj % time = i
                  tdiff = abs(xtime % array(i) - config_restart_time)
               end if
            end do

            tdiff = xtime % array(input_obj % time)

            deallocate(xtime % ioinfo)
            deallocate(xtime % array)
         end if

         call dmpar_bcast_int(domain % dminfo, input_obj % time)
         call dmpar_bcast_real(domain % dminfo, tdiff)

         write(0,*) 'Restarting model from time ', tdiff
    
      end if


      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      ! Do the actual work of reading all fields in from the input or restart file
      ! For each field:
      !   1) Each process reads a contiguous range of cell/edge/vertex indices, which
      !      may not correspond with the cells/edges/vertices that are owned by the
      !      process
      !   2) All processes then send the global indices that were read to the 
      !      processes that own those indices based on 
      !      {send,recv}{Cell,Edge,Vertex,VertLevel}List
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      call read_and_distribute_fields(domain % dminfo, input_obj, domain % blocklist, &
                                      readCellStart, nReadCells, readEdgeStart, nReadEdges, readVertexStart, nReadVertices, &
                                      readVertLevelStart, nReadVertLevels, &
                                      sendCellList, recvCellList, sendEdgeList, recvEdgeList, sendVertexList, recvVertexList, &
                                      sendVertLevelList, recvVertLevelList) 


      call io_input_finalize(input_obj, domain % dminfo)

   
      !
      ! Rename vertices in cellsOnCell, edgesOnCell, etc. to local indices
      !
      allocate(cellIDSorted(2,domain % blocklist % mesh % nCells))
      allocate(edgeIDSorted(2,domain % blocklist % mesh % nEdges))
      allocate(vertexIDSorted(2,domain % blocklist % mesh % nVertices))

      do i=1,domain % blocklist % mesh % nCells
         cellIDSorted(1,i) = domain % blocklist % mesh % indexToCellID % array(i)
         cellIDSorted(2,i) = i
      end do
      call quicksort(block_graph_2Halo % nVerticesTotal, cellIDSorted)

      do i=1,domain % blocklist % mesh % nEdges
         edgeIDSorted(1,i) = domain % blocklist % mesh % indexToEdgeID % array(i)
         edgeIDSorted(2,i) = i
      end do
      call quicksort(nlocal_edges, edgeIDSorted)

      do i=1,domain % blocklist % mesh % nVertices
         vertexIDSorted(1,i) = domain % blocklist % mesh % indexToVertexID % array(i)
         vertexIDSorted(2,i) = i
      end do
      call quicksort(nlocal_vertices, vertexIDSorted)


      do i=1,domain % blocklist % mesh % nCells
         do j=1,domain % blocklist % mesh % nEdgesOnCell % array(i)

            k = binary_search(cellIDSorted, 2, 1, domain % blocklist % mesh % nCells, &
                              domain % blocklist % mesh % cellsOnCell % array(j,i))
            if (k <= domain % blocklist % mesh % nCells) then
               domain % blocklist % mesh % cellsOnCell % array(j,i) = cellIDSorted(2,k)
            else
               domain % blocklist % mesh % cellsOnCell % array(j,i) = 0
            end if

            k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
                              domain % blocklist % mesh % edgesOnCell % array(j,i))
            if (k <= domain % blocklist % mesh % nEdges) then
               domain % blocklist % mesh % edgesOnCell % array(j,i) = edgeIDSorted(2,k)
            else
               domain % blocklist % mesh % edgesOnCell % array(j,i) = 0
            end if

            k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
                              domain % blocklist % mesh % verticesOnCell % array(j,i))
            if (k <= domain % blocklist % mesh % nVertices) then
               domain % blocklist % mesh % verticesOnCell % array(j,i) = vertexIDSorted(2,k)
            else
               domain % blocklist % mesh % verticesOnCell % array(j,i) = 0
            end if

         end do
      end do

      do i=1,domain % blocklist % mesh % nEdges
         do j=1,2

            k = binary_search(cellIDSorted, 2, 1, domain % blocklist % mesh % nCells, &
                              domain % blocklist % mesh % cellsOnEdge % array(j,i))
            if (k <= domain % blocklist % mesh % nCells) then
               domain % blocklist % mesh % cellsOnEdge % array(j,i) = cellIDSorted(2,k)
            else
               domain % blocklist % mesh % cellsOnEdge % array(j,i) = 0
            end if

            k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
                              domain % blocklist % mesh % verticesOnEdge % array(j,i))
            if (k <= domain % blocklist % mesh % nVertices) then
               domain % blocklist % mesh % verticesOnEdge % array(j,i) = vertexIDSorted(2,k)
            else
               domain % blocklist % mesh % verticesOnEdge % array(j,i) = 0
            end if

         end do

         do j=1,domain % blocklist % mesh % nEdgesOnEdge % array(i)

            k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
                              domain % blocklist % mesh % edgesOnEdge % array(j,i))
            if (k <= domain % blocklist % mesh % nEdges) then
               domain % blocklist % mesh % edgesOnEdge % array(j,i) = edgeIDSorted(2,k)
            else
               domain % blocklist % mesh % edgesOnEdge % array(j,i) = 0
            end if

         end do
      end do

      do i=1,domain % blocklist % mesh % nVertices
         do j=1,vertexDegree

            k = binary_search(cellIDSorted, 2, 1, domain % blocklist % mesh % nCells, &
                              domain % blocklist % mesh % cellsOnVertex % array(j,i))
            if (k <= domain % blocklist % mesh % nCells) then
               domain % blocklist % mesh % cellsOnVertex % array(j,i) = cellIDSorted(2,k)
            else
               domain % blocklist % mesh % cellsOnVertex % array(j,i) = 0
            end if

            k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
                              domain % blocklist % mesh % edgesOnVertex % array(j,i))
            if (k <= domain % blocklist % mesh % nEdges) then
               domain % blocklist % mesh % edgesOnVertex % array(j,i) = edgeIDSorted(2,k)
            else
               domain % blocklist % mesh % edgesOnVertex % array(j,i) = 0
            end if

         end do
      end do

      deallocate(cellIDSorted)
      deallocate(edgeIDSorted)
      deallocate(vertexIDSorted)


      !
      ! Work out halo exchange lists for cells, edges, and vertices
      !
      call dmpar_get_owner_list(domain % dminfo, &
                                block_graph_2Halo % nVertices, block_graph_2Halo % nVerticesTotal, &
                                block_graph_2Halo % vertexID(1:block_graph_2Halo % nVertices), block_graph_2Halo % vertexID, &
                                domain % blocklist % parinfo % cellsToSend, domain % blocklist % parinfo % cellsToRecv)

      call dmpar_get_owner_list(domain % dminfo, &
                                ghostEdgeStart-1, nlocal_edges, &
                                local_edge_list(1:ghostEdgeStart-1), local_edge_list, &
                                domain % blocklist % parinfo % edgesToSend, domain % blocklist % parinfo % edgesToRecv)

      call dmpar_get_owner_list(domain % dminfo, &
                                ghostVertexStart-1, nlocal_vertices, &
                                local_vertex_list(1:ghostVertexStart-1), local_vertex_list, &
                                domain % blocklist % parinfo % verticesToSend, domain % blocklist % parinfo % verticesToRecv)

      domain % blocklist % mesh % nCellsSolve = block_graph_2Halo % nVertices
      domain % blocklist % mesh % nEdgesSolve = ghostEdgeStart-1
      domain % blocklist % mesh % nVerticesSolve = ghostVertexStart-1
      domain % blocklist % mesh % nVertLevelsSolve = domain % blocklist % mesh % nVertLevels   ! No vertical decomp yet...

   
      !
      ! Deallocate fields, graphs, and other memory
      !
      deallocate(indexToCellIDField % ioinfo)
      deallocate(indexToCellIDField % array)
#ifdef HAVE_ZOLTAN
#ifdef _MPI 
      deallocate(xCellField % ioinfo)
      deallocate(xCellField % array)
      deallocate(yCellField % ioinfo)
      deallocate(yCellField % array)
      deallocate(zCellField % ioinfo)
      deallocate(zCellField % array)
#endif
#endif
      deallocate(indexToEdgeIDField % ioinfo)
      deallocate(indexToEdgeIDField % array)
      deallocate(indexToVertexIDField % ioinfo)
      deallocate(indexToVertexIDField % array)
      deallocate(cellsOnCellField % ioinfo)
      deallocate(cellsOnCellField % array)
      deallocate(edgesOnCellField % ioinfo)
      deallocate(edgesOnCellField % array)
      deallocate(verticesOnCellField % ioinfo)
      deallocate(verticesOnCellField % array)
      deallocate(cellsOnEdgeField % ioinfo)
      deallocate(cellsOnEdgeField % array)
      deallocate(cellsOnVertexField % ioinfo)
      deallocate(cellsOnVertexField % array)
      deallocate(cellsOnCell_0Halo)
      deallocate(nEdgesOnCell_0Halo)
      deallocate(indexToCellID_0Halo)
      deallocate(cellsOnEdge_2Halo)
      deallocate(cellsOnVertex_2Halo)
      deallocate(edgesOnCell_2Halo)
      deallocate(verticesOnCell_2Halo)
      deallocate(block_graph_0Halo % vertexID)
      deallocate(block_graph_0Halo % nAdjacent)
      deallocate(block_graph_0Halo % adjacencyList)
#ifdef HAVE_ZOLTAN
#ifdef _MPI
      deallocate(xCell)
      deallocate(yCell)
      deallocate(zCell)
#endif
#endif
   end subroutine input_state_for_domain


   subroutine read_and_distribute_fields(dminfo, input_obj, block, &
                                     readCellsStart, readCellsCount, &
                                     readEdgesStart, readEdgesCount, &
                                     readVerticesStart, readVerticesCount, &
                                     readVertLevelsStart, readVertLevelsCount, &
                                     sendCellsList, recvCellsList, &
                                     sendEdgesList, recvEdgesList, &
                                     sendVerticesList, recvVerticesList, &
                                     sendVertLevelsList, recvVertLevelsList) 
      
      implicit none

      type (dm_info), intent(in) :: dminfo
      type (io_input_object), intent(in) :: input_obj
      type (block_type), intent(inout) :: block
      integer, intent(in) :: readCellsStart, readCellsCount, readEdgesStart, readEdgesCount, readVerticesStart, readVerticesCount
      integer, intent(in) :: readVertLevelsStart, readVertLevelsCount
      type (exchange_list), pointer :: sendCellsList, recvCellsList
      type (exchange_list), pointer :: sendEdgesList, recvEdgesList
      type (exchange_list), pointer :: sendVerticesList, recvVerticesList
      type (exchange_list), pointer :: sendVertLevelsList, recvVertLevelsList

      type (field1dInteger) :: int1d
      type (field2dInteger) :: int2d
      type (field0dReal) :: real0d
      type (field1dReal) :: real1d
      type (field2dReal) :: real2d
      type (field3dReal) :: real3d

      integer :: i1, i2, i3, i4

      integer, dimension(:), pointer :: super_int1d
      integer, dimension(:,:), pointer :: super_int2d
      real (kind=RKIND) :: super_real0d
      real (kind=RKIND), dimension(:), pointer :: super_real1d
      real (kind=RKIND), dimension(:,:), pointer :: super_real2d
      real (kind=RKIND), dimension(:,:,:), pointer :: super_real3d

      integer :: k

      allocate(int1d % ioinfo)
      allocate(int2d % ioinfo)
      allocate(real0d % ioinfo)
      allocate(real1d % ioinfo)
      allocate(real2d % ioinfo)
      allocate(real3d % ioinfo)


#include "io_input_fields.inc"

   end subroutine read_and_distribute_fields



   subroutine io_input_init(input_obj, dminfo)
 
      implicit none

      type (io_input_object), intent(inout) :: input_obj
      type (dm_info), intent(in) :: dminfo
 
      include 'netcdf.inc'
 
      integer :: nferr
 
 
#ifdef OFFSET64BIT
      nferr = nf_open(trim(input_obj % filename), ior(NF_SHARE,NF_64BIT_OFFSET), input_obj % rd_ncid)
#else
      nferr = nf_open(trim(input_obj % filename), NF_SHARE, input_obj % rd_ncid)
#endif

      if (nferr /= NF_NOERR) then
         write(0,*) ' '
         if (config_do_restart) then
            write(0,*) 'Error opening restart file ''', trim(input_obj % filename), ''''
         else
            write(0,*) 'Error opening input file ''', trim(input_obj % filename), ''''
         end if
         write(0,*) ' '
         call dmpar_abort(dminfo)
      end if
 
#include "netcdf_read_ids.inc"

#ifdef EXPAND_LEVELS
      if (.not. config_do_restart) then
         input_obj % rdLocalnVertLevels = EXPAND_LEVELS
         write(0,*) 'Expanding nVertLevels to ',input_obj % rdLocalnVertLevels,' by duplicating the first level.'
      end if
#endif

   end subroutine io_input_init

  
   subroutine io_input_get_dimension(input_obj, dimname, dimsize)

      implicit none

      type (io_input_object), intent(in) :: input_obj
      character (len=*), intent(in) :: dimname
      integer, intent(out) :: dimsize

#include "get_dimension_by_name.inc"

   end subroutine io_input_get_dimension


   subroutine io_input_field0dReal(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field0dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1
 
      start1(1) = 1
      count1(1) = 1

#include "input_field0dreal.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % scalar)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
#endif
 
   end subroutine io_input_field0dReal


   subroutine io_input_field1dReal(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field1dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1
 
      start1(1) = field % ioinfo % start(1)
      count1(1) = field % ioinfo % count(1)

      !
      ! Special case: we may want to read the xtime variable across the
      !   time dimension as a 1d array.
      !
      if (trim(field % ioinfo % fieldName) == 'xtime') then
         varID = input_obj % rdVarIDxtime
      end if
 
#include "input_field1dreal.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % array)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % array)
#endif
 
   end subroutine io_input_field1dReal


   subroutine io_input_field2dReal(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field2dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(2) :: start2, count2
 
      start2(1) = field % ioinfo % start(1)
      start2(2) = field % ioinfo % start(2)
      count2(1) = field % ioinfo % count(1)
      count2(2) = field % ioinfo % count(2)
 
#include "input_field2dreal.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start2, count2, field % array)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start2, count2, field % array)
#endif

   end subroutine io_input_field2dReal


   subroutine io_input_field3dReal(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field3dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(3) :: start3, count3
 
      start3(1) = field % ioinfo % start(1)
      start3(2) = field % ioinfo % start(2)
      start3(3) = field % ioinfo % start(3)
      count3(1) = field % ioinfo % count(1)
      count3(2) = field % ioinfo % count(2)
      count3(3) = field % ioinfo % count(3)
 
#include "input_field3dreal.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start3, count3, field % array)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start3, count3, field % array)
#endif

   end subroutine io_input_field3dReal


   subroutine io_input_field0dReal_time(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field0dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1
 
      start1(1) = input_obj % time
      count1(1) = 1
 
#include "input_field0dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start1, count1, field % scalar)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start1, count1, field % scalar)
#endif

   end subroutine io_input_field0dReal_time


   subroutine io_input_field1dReal_time(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field1dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(2) :: start2, count2
 
      start2(1) = field % ioinfo % start(1)
      start2(2) = input_obj % time
      count2(1) = field % ioinfo % count(1)
      count2(2) = 1
 
#include "input_field1dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start2, count2, field % array)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start2, count2, field % array)
#endif

   end subroutine io_input_field1dReal_time


   subroutine io_input_field2dReal_time(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field2dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(3) :: start3, count3
 
      start3(1) = field % ioinfo % start(1)
      start3(2) = field % ioinfo % start(2)
      start3(3) = input_obj % time
      count3(1) = field % ioinfo % count(1)
      count3(2) = field % ioinfo % count(2)
      count3(3) = 1
 
#include "input_field2dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start3, count3, field % array)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start3, count3, field % array)
#endif

   end subroutine io_input_field2dReal_time


   subroutine io_input_field3dReal_time(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field3dReal), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(4) :: start4, count4
 
      start4(1) = field % ioinfo % start(1)
      start4(2) = field % ioinfo % start(2)
      start4(3) = field % ioinfo % start(3)
      start4(4) = input_obj % time
      count4(1) = field % ioinfo % count(1)
      count4(2) = field % ioinfo % count(2)
      count4(3) = field % ioinfo % count(3)
      count4(4) = 1
 
#include "input_field3dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_get_vara_double(input_obj % rd_ncid, varID, start4, count4, field % array)
#else
      nferr = nf_get_vara_real(input_obj % rd_ncid, varID, start4, count4, field % array)
#endif

   end subroutine io_input_field3dReal_time


   subroutine io_input_field1dInteger(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field1dInteger), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1
 
      start1(1) = field % ioinfo % start(1)
      count1(1) = field % ioinfo % count(1)
      
#include "input_field1dinteger.inc"

      nferr = nf_get_vara_int(input_obj % rd_ncid, varID, start1, count1, field % array)
 
   end subroutine io_input_field1dInteger


   subroutine io_input_field2dInteger(input_obj, field)
 
      implicit none

      type (io_input_object), intent(in) :: input_obj      
      type (field2dInteger), intent(inout) :: field
 
      include 'netcdf.inc'
 
      integer :: nferr
      integer :: varID
      integer, dimension(2) :: start2, count2
 
      start2(1) = field % ioinfo % start(1)
      start2(2) = field % ioinfo % start(2)
      count2(1) = field % ioinfo % count(1)
      count2(2) = field % ioinfo % count(2)

#include "input_field2dinteger.inc"

      nferr = nf_get_vara_int(input_obj % rd_ncid, varID, start2, count2, field % array)

   end subroutine io_input_field2dInteger


   subroutine io_input_finalize(input_obj, dminfo)
 
      implicit none
 
      type (io_input_object), intent(inout) :: input_obj
      type (dm_info), intent(in) :: dminfo

      include 'netcdf.inc'
 
      integer :: nferr
 
      nferr = nf_close(input_obj % rd_ncid)
 
   end subroutine io_input_finalize
 
end module io_input
-------------- next part --------------
module io_output

   use grid_types
   use dmpar
   use sort
   use configure

   integer, parameter :: OUTPUT = 1
   integer, parameter :: RESTART = 2
 
   type io_output_object
      integer :: wr_ncid
      character (len=1024) :: filename

      integer :: time

      integer :: stream

#include "io_output_obj_decls.inc"

      logical :: validExchangeLists
      type (exchange_list), pointer :: sendCellsList, recvCellsList
      type (exchange_list), pointer :: sendEdgesList, recvEdgesList
      type (exchange_list), pointer :: sendVerticesList, recvVerticesList
      type (exchange_list), pointer :: sendVertLevelsList, recvVertLevelsList
   end type io_output_object


   interface io_output_field
      module procedure io_output_field0dReal
      module procedure io_output_field1dReal
      module procedure io_output_field2dReal
      module procedure io_output_field3dReal
      module procedure io_output_field1dInteger
      module procedure io_output_field2dInteger
   end interface io_output_field

   interface io_output_field_time
      module procedure io_output_field0dReal_time
      module procedure io_output_field1dReal_time
      module procedure io_output_field2dReal_time
      module procedure io_output_field3dReal_time
   end interface io_output_field_time
 

   contains

 
   subroutine output_state_init(output_obj, domain, stream)

      implicit none

      type (io_output_object), intent(inout) :: output_obj
      type (domain_type), intent(in) :: domain
      character (len=*) :: stream

      type (block_type), pointer :: block_ptr
#include "output_dim_actual_decls.inc"

      block_ptr => domain % blocklist
      nullify(output_obj % sendCellsList)
      nullify(output_obj % recvCellsList)
      nullify(output_obj % sendEdgesList)
      nullify(output_obj % recvEdgesList)
      nullify(output_obj % sendVerticesList)
      nullify(output_obj % recvVerticesList)
      nullify(output_obj % sendVertLevelsList)
      nullify(output_obj % recvVertLevelsList)
      output_obj % validExchangeLists = .false.

#include "output_dim_inits.inc"

      call dmpar_sum_int(domain % dminfo, block_ptr % mesh % nCellsSolve, nCellsGlobal) 
      call dmpar_sum_int(domain % dminfo, block_ptr % mesh % nEdgesSolve, nEdgesGlobal) 
      call dmpar_sum_int(domain % dminfo, block_ptr % mesh % nVerticesSolve, nVerticesGlobal) 
      nVertLevelsGlobal = block_ptr % mesh % nVertLevels

      if (trim(stream) == 'OUTPUT') then
         output_obj % filename = trim(config_output_name)
         output_obj % stream = OUTPUT
      else if (trim(stream) == 'RESTART') then
         output_obj % filename = trim(config_restart_name)
         output_obj % stream = RESTART
      end if

      ! For now, we assume that a domain consists only of one block,
      !   although in future, work needs to be done to write model state
      !   from many distributed blocks
      call io_output_init(output_obj, domain % dminfo, &
#include "output_dim_actual_args.inc"
                         )
   
   end subroutine output_state_init


   subroutine output_state_for_domain(output_obj, domain, itime)
   
      implicit none
   
      type (io_output_object), intent(inout) :: output_obj
      type (domain_type), intent(inout) :: domain
      integer, intent(in) :: itime

      integer :: i, j
      integer :: nCellsGlobal
      integer :: nEdgesGlobal
      integer :: nVerticesGlobal
      integer :: nVertLevelsGlobal
      integer, dimension(:), pointer :: neededCellList
      integer, dimension(:), pointer :: neededEdgeList
      integer, dimension(:), pointer :: neededVertexList
      integer, dimension(:), pointer :: neededVertLevelList
      integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, verticesOnCell, &
                                          cellsOnEdge, verticesOnEdge, edgesOnEdge, cellsOnVertex, edgesOnVertex
      integer, dimension(:,:), pointer :: cellsOnCell_save, edgesOnCell_save, verticesOnCell_save, &
                                          cellsOnEdge_save, verticesOnEdge_save, edgesOnEdge_save, &
                                          cellsOnVertex_save, edgesOnVertex_save
      type (field1dInteger) :: int1d
      type (field2dInteger) :: int2d
      type (field0dReal) :: real0d
      type (field1dReal) :: real1d
      type (field2dReal) :: real2d
      type (field3dReal) :: real3d

      integer :: i1, i2, i3, i4

      integer, dimension(:), pointer :: super_int1d
      integer, dimension(:,:), pointer :: super_int2d
      real :: super_real0d
      real (kind=RKIND), dimension(:), pointer :: super_real1d
      real (kind=RKIND), dimension(:,:), pointer :: super_real2d
      real (kind=RKIND), dimension(:,:,:), pointer :: super_real3d

      output_obj % time = itime

      allocate(int1d % ioinfo)
      allocate(int2d % ioinfo)
      allocate(real0d % ioinfo)
      allocate(real1d % ioinfo)
      allocate(real2d % ioinfo)
      allocate(real3d % ioinfo)

      call dmpar_sum_int(domain % dminfo, domain % blocklist % mesh % nCellsSolve, nCellsGlobal)
      call dmpar_sum_int(domain % dminfo, domain % blocklist % mesh % nEdgesSolve, nEdgesGlobal)
      call dmpar_sum_int(domain % dminfo, domain % blocklist % mesh % nVerticesSolve, nVerticesGlobal)
      nVertLevelsGlobal = domain % blocklist % mesh % nVertLevels

      allocate(cellsOnCell(domain % blocklist % mesh % maxEdges, domain % blocklist % mesh % nCellsSolve))
      allocate(edgesOnCell(domain % blocklist % mesh % maxEdges, domain % blocklist % mesh % nCellsSolve))
      allocate(verticesOnCell(domain % blocklist % mesh % maxEdges, domain % blocklist % mesh % nCellsSolve))
      allocate(cellsOnEdge(2, domain % blocklist % mesh % nEdgesSolve))
      allocate(verticesOnEdge(2, domain % blocklist % mesh % nEdgesSolve))
      allocate(edgesOnEdge(2 * domain % blocklist % mesh % maxEdges, domain % blocklist % mesh % nEdgesSolve))
      allocate(cellsOnVertex(domain % blocklist % mesh % vertexDegree, domain % blocklist % mesh % nVerticesSolve))
      allocate(edgesOnVertex(domain % blocklist % mesh % vertexDegree, domain % blocklist % mesh % nVerticesSolve))


      !
      ! Convert connectivity information from local to global indices
      !
      do i=1,domain % blocklist % mesh % nCellsSolve
         do j=1,domain % blocklist % mesh % nEdgesOnCell % array(i)
            cellsOnCell(j,i) = domain % blocklist % mesh % indexToCellID % array( &
                                                                           domain % blocklist % mesh % cellsOnCell % array(j,i))
            edgesOnCell(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
                                                                           domain % blocklist % mesh % edgesOnCell % array(j,i))
            verticesOnCell(j,i) = domain % blocklist % mesh % indexToVertexID % array( &
                                                                           domain % blocklist % mesh % verticesOnCell % array(j,i))
         end do
         do j=domain % blocklist % mesh % nEdgesOnCell % array(i)+1,domain % blocklist % mesh % maxEdges
            cellsOnCell(j,i) = domain % blocklist % mesh % indexToCellID % array( &
                                                                           domain % blocklist % mesh % nEdgesOnCell % array(i))
            edgesOnCell(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
                                                                           domain % blocklist % mesh % nEdgesOnCell % array(i))
            verticesOnCell(j,i) = domain % blocklist % mesh % indexToVertexID % array( &
                                                                           domain % blocklist % mesh % nEdgesOnCell % array(i))
         end do
      end do
      do i=1,domain % blocklist % mesh % nEdgesSolve
         cellsOnEdge(1,i) = domain % blocklist % mesh % indexToCellID % array(domain % blocklist % mesh % cellsOnEdge % array(1,i))
         cellsOnEdge(2,i) = domain % blocklist % mesh % indexToCellID % array(domain % blocklist % mesh % cellsOnEdge % array(2,i))
         verticesOnEdge(1,i) = domain % blocklist % mesh % indexToVertexID % array( &
                                                                           domain % blocklist % mesh % verticesOnEdge % array(1,i))
         verticesOnEdge(2,i) = domain % blocklist % mesh % indexToVertexID % array( &
                                                                           domain % blocklist % mesh % verticesOnEdge % array(2,i))
         do j=1,domain % blocklist % mesh % nEdgesOnEdge % array(i)
            edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
                                                                           domain % blocklist % mesh % edgesOnEdge % array(j,i))
         end do
         do j=domain % blocklist % mesh % nEdgesOnEdge % array(i)+1,2*domain % blocklist % mesh % maxEdges
            edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
                                                                           domain % blocklist % mesh % nEdgesOnEdge % array(i))
         end do
      end do
      do i=1,domain % blocklist % mesh % nVerticesSolve
         do j=1,domain % blocklist % mesh % vertexDegree
            cellsOnVertex(j,i) = domain % blocklist % mesh % indexToCellID % array( &
                                                                           domain % blocklist % mesh % cellsOnVertex % array(j,i))
            edgesOnVertex(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
                                                                           domain % blocklist % mesh % edgesOnVertex % array(j,i))
         end do
      end do

      if (domain % dminfo % my_proc_id == 0) then
         allocate(neededCellList(nCellsGlobal))
         allocate(neededEdgeList(nEdgesGlobal))
         allocate(neededVertexList(nVerticesGlobal))
         allocate(neededVertLevelList(nVertLevelsGlobal))
         do i=1,nCellsGlobal
            neededCellList(i) = i
         end do
         do i=1,nEdgesGlobal
            neededEdgeList(i) = i
         end do
         do i=1,nVerticesGlobal
            neededVertexList(i) = i
         end do
         do i=1,nVertLevelsGlobal
            neededVertLevelList(i) = i
         end do
      else
         allocate(neededCellList(0))
         allocate(neededEdgeList(0))
         allocate(neededVertexList(0))
         allocate(neededVertLevelList(0))
      end if

      if (.not. output_obj % validExchangeLists) then
         call dmpar_get_owner_list(domain % dminfo, &
                                   domain % blocklist % mesh % nCellsSolve, size(neededCellList), &
                                   domain % blocklist % mesh % indexToCellID % array, neededCellList, &
                                   output_obj % sendCellsList, output_obj % recvCellsList)

         call dmpar_get_owner_list(domain % dminfo, &
                                   domain % blocklist % mesh % nEdgesSolve, size(neededEdgeList), &
                                   domain % blocklist % mesh % indexToEdgeID % array, neededEdgeList, &
                                   output_obj % sendEdgesList, output_obj % recvEdgesList)

         call dmpar_get_owner_list(domain % dminfo, &
                                   domain % blocklist % mesh % nVerticesSolve, size(neededVertexList), &
                                   domain % blocklist % mesh % indexToVertexID % array, neededVertexList, &
                                   output_obj % sendVerticesList, output_obj % recvVerticesList)

         call dmpar_get_owner_list(domain % dminfo, &
                                   size(neededVertLevelList), size(neededVertLevelList), &
                                   neededVertLevelList, neededVertLevelList, &
                                   output_obj % sendVertLevelsList, output_obj % recvVertLevelsList)

         output_obj % validExchangeLists = .true.
      end if

      deallocate(neededCellList)
      deallocate(neededEdgeList)
      deallocate(neededVertexList)

      cellsOnCell_save => domain % blocklist % mesh % cellsOnCell % array
      edgesOnCell_save => domain % blocklist % mesh % edgesOnCell % array
      verticesOnCell_save => domain % blocklist % mesh % verticesOnCell % array
      cellsOnEdge_save => domain % blocklist % mesh % cellsOnEdge % array
      verticesOnEdge_save => domain % blocklist % mesh % verticesOnEdge % array
      edgesOnEdge_save => domain % blocklist % mesh % edgesOnEdge % array
      cellsOnVertex_save => domain % blocklist % mesh % cellsOnVertex % array
      edgesOnVertex_save => domain % blocklist % mesh % edgesOnVertex % array

      domain % blocklist % mesh % cellsOnCell % array => cellsOnCell
      domain % blocklist % mesh % edgesOnCell % array => edgesOnCell
      domain % blocklist % mesh % verticesOnCell % array => verticesOnCell
      domain % blocklist % mesh % cellsOnEdge % array => cellsOnEdge
      domain % blocklist % mesh % verticesOnEdge % array => verticesOnEdge
      domain % blocklist % mesh % edgesOnEdge % array => edgesOnEdge
      domain % blocklist % mesh % cellsOnVertex % array => cellsOnVertex
      domain % blocklist % mesh % edgesOnVertex % array => edgesOnVertex

#include "io_output_fields.inc"

      domain % blocklist % mesh % cellsOnCell % array => cellsOnCell_save
      domain % blocklist % mesh % edgesOnCell % array => edgesOnCell_save
      domain % blocklist % mesh % verticesOnCell % array => verticesOnCell_save
      domain % blocklist % mesh % cellsOnEdge % array => cellsOnEdge_save
      domain % blocklist % mesh % verticesOnEdge % array => verticesOnEdge_save
      domain % blocklist % mesh % edgesOnEdge % array => edgesOnEdge_save
      domain % blocklist % mesh % cellsOnVertex % array => cellsOnVertex_save
      domain % blocklist % mesh % edgesOnVertex % array => edgesOnVertex_save

      deallocate(cellsOnCell)
      deallocate(edgesOnCell)
      deallocate(verticesOnCell)
      deallocate(cellsOnEdge)
      deallocate(verticesOnEdge)
      deallocate(edgesOnEdge)
      deallocate(cellsOnVertex)
      deallocate(edgesOnVertex)

   end subroutine output_state_for_domain


   subroutine output_state_finalize(output_obj, dminfo)

      implicit none

      type (io_output_object), intent(inout) :: output_obj
      type (dm_info), intent(in) :: dminfo

      call io_output_finalize(output_obj, dminfo)

   end subroutine output_state_finalize


   subroutine io_output_init( output_obj, &
                              dminfo, &
#include "dim_dummy_args.inc"
                            )
 
      implicit none
 
      include 'netcdf.inc'
 
      type (io_output_object), intent(inout) :: output_obj
      type (dm_info), intent(in) :: dminfo
#include "dim_dummy_decls.inc"
 
      integer :: nferr
      integer, dimension(10) :: dimlist
 
      if (dminfo % my_proc_id == 0) then
#ifdef OFFSET64BIT
      nferr = nf_create(trim(output_obj % filename), ior(NF_CLOBBER,NF_64BIT_OFFSET), output_obj % wr_ncid)
#else
      nferr = nf_create(trim(output_obj % filename), NF_CLOBBER, output_obj % wr_ncid)
#endif
 
#include "netcdf_def_dims_vars.inc"
 
      nferr = nf_enddef(output_obj % wr_ncid)
      end if
 
   end subroutine io_output_init


   subroutine io_output_field0dReal(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field0dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1

      start1(1) = 1
      count1(1) = 1

#include "output_field0dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start1, count1, field % scalar)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field0dReal


   subroutine io_output_field1dReal(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field1dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1

      start1(1) = field % ioinfo % start(1)
      count1(1) = field % ioinfo % count(1)

#include "output_field1dreal.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, VarID, start1, count1, field % array)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, VarID, start1, count1, field % array)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field1dReal
 
 
   subroutine io_output_field2dReal(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field2dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(2) :: start2, count2

      start2(1) = field % ioinfo % start(1)
      start2(2) = field % ioinfo % start(2)
      count2(1) = field % ioinfo % count(1)
      count2(2) = field % ioinfo % count(2)

#include "output_field2dreal.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start2, count2, field % array)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start2, count2, field % array)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field2dReal
 
 
   subroutine io_output_field3dReal(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field3dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(3) :: start3, count3

      start3(1) = field % ioinfo % start(1)
      start3(2) = field % ioinfo % start(2)
      start3(3) = field % ioinfo % start(3)
      count3(1) = field % ioinfo % count(1)
      count3(2) = field % ioinfo % count(2)
      count3(3) = field % ioinfo % count(3)

#include "output_field3dreal.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start3, count3, field % array)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start3, count3, field % array)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field3dReal


   subroutine io_output_field0dReal_time(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field0dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1

      start1(1) = output_obj % time
      count1(1) = 1

#include "output_field0dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start1, count1, field % scalar)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start1, count1, field % scalar)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field0dReal_time


   subroutine io_output_field1dReal_time(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field1dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(2) :: start2, count2

      start2(1) = field % ioinfo % start(1)
      start2(2) = output_obj % time
      count2(1) = field % ioinfo % count(1)
      count2(2) = 1

#include "output_field1dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start2, count2, field % array)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start2, count2, field % array)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field1dReal_time


   subroutine io_output_field2dReal_time(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field2dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(3) :: start3, count3

      start3(1) = field % ioinfo % start(1)
      start3(2) = field % ioinfo % start(2)
      start3(3) = output_obj % time
      count3(1) = field % ioinfo % count(1)
      count3(2) = field % ioinfo % count(2)
      count3(3) = 1

#include "output_field2dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start3, count3, field % array)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start3, count3, field % array)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field2dReal_time


   subroutine io_output_field3dReal_time(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field3dReal), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(4) :: start4, count4

      start4(1) = field % ioinfo % start(1)
      start4(2) = field % ioinfo % start(2)
      start4(3) = field % ioinfo % start(3)
      start4(4) = output_obj % time
      count4(1) = field % ioinfo % count(1)
      count4(2) = field % ioinfo % count(2)
      count4(3) = field % ioinfo % count(3)
      count4(4) = 1

#include "output_field3dreal_time.inc"

#if (RKIND == 8)
      nferr = nf_put_vara_double(output_obj % wr_ncid, varID, start4, count4, field % array)
#else
      nferr = nf_put_vara_real(output_obj % wr_ncid, varID, start4, count4, field % array)
#endif
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field3dReal_time


   subroutine io_output_field1dInteger(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field1dInteger), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(1) :: start1, count1

      start1(1) = field % ioinfo % start(1)
      count1(1) = field % ioinfo % count(1)

#include "output_field1dinteger.inc"

      nferr = nf_put_vara_int(output_obj % wr_ncid, varID, start1, count1, field % array)
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field1dInteger


   subroutine io_output_field2dInteger(output_obj, field)

      implicit none

      type (io_output_object), intent(in) :: output_obj
      type (field2dInteger), intent(inout) :: field

      include 'netcdf.inc'

      integer :: nferr
      integer :: varID
      integer, dimension(2) :: start2, count2

      start2(1) = field % ioinfo % start(1)
      start2(2) = field % ioinfo % start(2)
      count2(1) = field % ioinfo % count(1)
      count2(2) = field % ioinfo % count(2)

#include "output_field2dinteger.inc"

      nferr = nf_put_vara_int(output_obj % wr_ncid, varID, start2, count2, field % array)
 
      nferr = nf_sync(output_obj % wr_ncid)

   end subroutine io_output_field2dInteger


   subroutine io_output_finalize(output_obj, dminfo)
 
      implicit none
 
      include 'netcdf.inc'

      type (io_output_object), intent(inout) :: output_obj
      type (dm_info), intent(in) :: dminfo
 
      integer :: nferr
 
      if (dminfo % my_proc_id == 0) then
      nferr = nf_close(output_obj % wr_ncid)
      end if
 
   end subroutine io_output_finalize
 
end module io_output
 


More information about the mpas-developers mailing list