createEmptyFile.c

/*********************************************************************

    File:           cdfdistiller.c

    Created:        14/3/2006

    Author:         Beth Fulton
                    CSIRO Division of Marine Resaech

    Purpose:        Main routine for distilling end of old output into
					input for new atlantis mse box model

    Arguments:      See Usage() routine below for explanation
					of command line arguments

    Returns:        int - 0 if sucessful

    Revisions:


*********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <signal.h>
#include <math.h>
#include <string.h>
#include <sjwlib.h>
#include <netcdf.h>
#include "cdfdistiller.h"

/****************************************************************************
Main run routine
*/
int createEmptyFile(char *dummyFileName, char *outputFileName, char *inputFileName)
{
	//char str[1000];
    cdfBoxModel bm;
    int fid, dummyfid, b;
    //, totnz;

    strcpy(bm.ncIfname, inputFileName);
	//strcpy(bm.dummyFileName,dummyFileName);
	strcpy(bm.ncOfname, outputFileName);


    /* Clear the box model - this assumes that it doesn't
     * already contain valid data. If so, allocated memory should
     * be freed first by whoever is calling this routine.
     */
    memset(&bm, 0, sizeof(cdfBoxModel));

    /* Write arguments to parameter string */
    sprintf(bm.params,"cdfDistiller -i %s -o %s ", inputFileName, outputFileName);

	/* Open input netCDF file */
    printf("Opening the input file %s\n", bm.ncIfname);
    fid = ncopen(inputFileName,NC_NOWRITE);

    /* Read the model geometry */
    ncattget(fid,NC_GLOBAL,"geometry",bm.geomIfname);

    /* This reads in the number of boxes and layers */
    readcdfBoxModelGeom(bm.geomIfname,&bm);

    /* Open the dummy netcdf file */
    dummyfid = ncopen(dummyFileName,NC_NOWRITE);

    /* Read the model variable info. Physical info must be first,
     * as it reads the number of water column and sediment layers
     * from the input netCDF file
     */

    /* Note the different files here.
     * The physical info is read in from the atlantis output file - this
     * sets the number of water column layers and boxes.
     *
     * Then the tracer information is read in from the dummy file. This is
     * used to ensure that there are tracers for all groups even those that
     * were turned 'off' in the original file.
     */
    readBMphysInfo(fid,&bm, 1);
	readBMTracerInfo(dummyfid,&bm, 1);
	readBMEpiInfo(dummyfid,&bm, 1);

	/* Allocate memory for variables. Note that the casts (FPTYPE ***)
     * are only needed to supress compiler warnings. Note also that the
     * allocation routines used here have parameter ordering opposite
     * to the way that the arrays are accessed. Thus if we want a
     * 3d array X, accessed X[k][j][i] (i varies fastest), then we say
     * X = ...alloc3d(ni,nj,nk), where ni is the number of elements
     * in the i dimension etc.
     */
    switch( sizeof(FPTYPE) ) {
	case sizeof(float):
	    bm.nom_dz = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.dz = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
		bm.wctr = (FPTYPE ***)f_alloc3d(bm.ntracer,bm.wcnz,bm.nbox);
		bm.sedtr = (FPTYPE ***)f_alloc3d(bm.ntracer,bm.sednz,bm.nbox);
	    bm.epi = (FPTYPE **)f_alloc2d(bm.nepi,bm.nbox);
	    bm.vol = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.por = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.hddrain = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.hdpool = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.efluxes = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.vfluxes = (FPTYPE **)f_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    break;
	case sizeof(double):
	    bm.nom_dz = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.dz = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
		bm.wctr = (FPTYPE ***)d_alloc3d(bm.ntracer,bm.wcnz,bm.nbox);
		bm.sedtr = (FPTYPE ***)d_alloc3d(bm.ntracer,bm.sednz,bm.nbox);
	    bm.epi = (FPTYPE **)d_alloc2d(bm.nepi,bm.nbox);
	    bm.vol = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.por = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.hddrain = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.hdpool = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.efluxes = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    bm.vfluxes = (FPTYPE **)d_alloc2d(bm.wcnz+bm.sednz,bm.nbox);
	    break;
	default:
	    quit("setupcdfBoxModel: Unknown size for FPTYPE\n");
		break;
    }

    /* Set up pointers to data from each box */
    for(b=0; b<bm.nbox; b++) {
		int k;
		bm.boxes[b].tr = bm.wctr[b];
		bm.boxes[b].sm.tr = bm.sedtr[b];
		bm.boxes[b].epi = bm.epi[b];
		bm.boxes[b].volume = bm.vol[b];
		bm.boxes[b].sm.volume = bm.vol[b]+bm.wcnz;
		bm.boxes[b].dz = bm.dz[b];
		bm.boxes[b].sm.dz = bm.dz[b]+bm.wcnz;
		bm.boxes[b].hdsource = bm.hddrain[b];
		bm.boxes[b].sm.hdsource = bm.hddrain[b]+bm.wcnz;
		bm.boxes[b].hdsink = bm.hdpool[b];
		bm.boxes[b].sm.hdsink = bm.hdpool[b]+bm.wcnz;
		bm.boxes[b].eflux = bm.efluxes[b];
		bm.boxes[b].vflux = bm.vfluxes[b];
		bm.boxes[b].sm.eflux = bm.efluxes[b]+bm.wcnz;
		bm.boxes[b].sm.vflux = bm.vfluxes[b]+bm.wcnz;
		for(k=0; k<bm.wcnz; k++){
			/* Set water column porosity values to 1 */
			bm.por[b][k] = 1.0;
		}
		/* Set sediment porosity */
		k = bm.boxes[b].sm.topk;
		bm.boxes[b].sm.porosity = bm.por[b]+bm.wcnz;
	}
	//totnz = bm.wcnz + bm.sednz;


	/* Test to see if can open output files */
    /* Open the general tracer output file */
	printf("createEmptyFile - Create general output file\n");
    ncopts = 0;
    bm.ncOfid = ncopen(outputFileName,NC_WRITE);
    if( bm.ncOfid == -1 ) {
    	/* Create file */
    	bm.ncOfid = createBMDataFile(outputFileName,&bm);
    	bm.ncOfid = ncopen(outputFileName,NC_WRITE);
		bm.ncOfdump = 0;

    } else {
		/* File does exist and don't want to replace or re-use it */
		quit("file with name %s already exists, either delete, rename or set reuse param in run.prm to 1\n", outputFileName);
	}

  //  bm.ncOfid = fid;
	ncopts = NC_VERBOSE | NC_FATAL;

#if 0
    /* Set signal to allow clean up on SIGTERM */
    signal(SIGTERM,setkill);
#endif

    /* Set output dump to zero (as want it to be first entry in new input file) */
   	bm.ncOfdump = 0;

    /* Write new input dump and close the input files */
    writeBMphysData(bm.ncOfid,0,&bm);
    writeBMTracerData(bm.ncOfid,0,&bm);
    writeBMEpiData(bm.ncOfid,0,&bm);

    ncclose(bm.ncOfid);

    ncclose(fid);
    ncclose(dummyfid);
    /* Free up all the memory */

    d_free2d(bm.nom_dz);
    d_free2d(bm.dz);
    d_free3d(bm.wctr);
    d_free3d(bm.sedtr);
    d_free2d(bm.epi);
    d_free2d(bm.vol);
    d_free2d(bm.por);
    d_free2d(bm.hddrain);
    d_free2d(bm.hdpool);
    d_free2d(bm.efluxes);
    d_free2d(bm.vfluxes);

    freeBMEpiInfo(&bm);
    freeBMTracerInfo(&bm);
    freeBMphysInfo(&bm);
    freecdfBoxModelGeom(&bm);

    return 0;

}