What is on Your Mind?

/* =============================================================================== */
/*     SUBROUTINE ONSET_OF_IMPERMEABILITY */
/*     Purpose:  This subroutine uses the Bisection Method to find the ambient */
/*     pressure and gas tension at the onset of impermeability for a given */
/*     compartment.  Source:  "Numerical Recipes in Fortran 77", */
/*     Cambridge University Press, 1992. */
/* =============================================================================== */

int onset_of_impermeability(real *starting_ambient_pressure, 
                              real *ending_ambient_pressure, 
                              real *rate,
                              integer *i)
{
    /* Local variables */
    static real time, last_diff_change, mid_range_nitrogen_pressure;
    static integer j;
    static real gas_tension_at_mid_range, 
           initial_inspired_n2_pressure, 
	       gradient_onset_of_imperm, 
           starting_gas_tension, 
           low_bound, 
	       initial_inspired_he_pressure, 
           high_bound_nitrogen_pressure, 
	       nitrogen_rate, 
           function_at_mid_range, 
           function_at_low_bound,
	       high_bound, 
           mid_range_helium_pressure, 
           mid_range_time, 
	       ending_gas_tension, 
           function_at_high_bound;

    static real mid_range_ambient_pressure, 
                high_bound_helium_pressure, 
	            helium_rate, 
                differential_change;

/* loop */
/* =============================================================================== */
/*     CALCULATIONS */
/*     First convert the Gradient for Onset of Impermeability to the diving */
/*     pressure units that are being used */
/* =============================================================================== */

    gradient_onset_of_imperm = gradient_onset_of_imperm_atm * units_factor;

/* =============================================================================== */
/*     ESTABLISH THE BOUNDS FOR THE ROOT SEARCH USING THE BISECTION METHOD */
/*     In this case, we are solving for time - the time when the ambient pressure */
/*     minus the gas tension will be equal to the Gradient for Onset of */
/*     Impermeabliity.  The low bound for time is set at zero and the high */
/*     bound is set at the elapsed time (segment time) it took to go from the */
/*     starting ambient pressure to the ending ambient pressure.  The desired */
/*     ambient pressure and gas tension at the onset of impermeability will */
/*     be found somewhere between these endpoints.  The algorithm checks to */
/*     make sure that the solution lies in between these bounds by first */
/*     computing the low bound and high bound function values. */
/* =============================================================================== */

    initial_inspired_he_pressure = 
        (*starting_ambient_pressure - water_vapor_pressure) * fraction_helium[mix_number - 1];
    initial_inspired_n2_pressure = 
        (*starting_ambient_pressure - water_vapor_pressure) * fraction_nitrogen[mix_number - 1];
    helium_rate = *rate * fraction_helium[mix_number - 1];
    nitrogen_rate = *rate * fraction_nitrogen[mix_number - 1];
    low_bound = 0.;
    high_bound = (*ending_ambient_pressure - *starting_ambient_pressure) / *rate;
    starting_gas_tension = 
        initial_helium_pressure[*i - 1] +
	    initial_nitrogen_pressure[*i - 1] + 
	    constant_pressure_other_gases;
    function_at_low_bound = 
        *starting_ambient_pressure - 
	    starting_gas_tension - 
        gradient_onset_of_imperm;
    high_bound_helium_pressure = 
        schreiner_equation__(&initial_inspired_he_pressure, 
                             &helium_rate, 
                             &high_bound, 
                             &helium_time_constant[*i - 1], 
                             &initial_helium_pressure[*i - 1]);
    high_bound_nitrogen_pressure = 
        schreiner_equation__(&initial_inspired_n2_pressure, 
                             &nitrogen_rate, 
                             &high_bound, 
                             &nitrogen_time_constant[*i - 1], 
                             &initial_nitrogen_pressure[*i - 1]);
    ending_gas_tension = 
        high_bound_helium_pressure + 
	    high_bound_nitrogen_pressure + 
	    constant_pressure_other_gases;
    function_at_high_bound = 
        *ending_ambient_pressure - 
	    ending_gas_tension - 
        gradient_onset_of_imperm;
    if (function_at_high_bound * function_at_low_bound >= 0.) {
	    printf("\nERROR! ROOT IS NOT WITHIN BRACKETS");
        pause();
    }

/* =============================================================================== */
/*     APPLY THE BISECTION METHOD IN SEVERAL ITERATIONS UNTIL A SOLUTION WITH */
/*     THE DESIRED ACCURACY IS FOUND */
/*     Note: the program allows for up to 100 iterations.  Normally an exit will */
/*     be made from the loop well before that number.  If, for some reason, the */
/*     program exceeds 100 iterations, there will be a pause to alert the user. */
/* =============================================================================== */

    if (function_at_low_bound < 0.) {
	    time = low_bound;
	    differential_change = high_bound - low_bound;
    } else {
	    time = high_bound;
	    differential_change = low_bound - high_bound;
    }
    for (j = 1; j <= 100; ++j) {
	    last_diff_change = differential_change;
	    differential_change = last_diff_change * .5;
	    mid_range_time = time + differential_change;
	    mid_range_ambient_pressure = *starting_ambient_pressure + *rate * mid_range_time;
	    mid_range_helium_pressure = 
            schreiner_equation__(&initial_inspired_he_pressure, 
                                 &helium_rate, 
                                 &mid_range_time, 
                                 &helium_time_constant[*i - 1], 
                                 &initial_helium_pressure[*i - 1]);
	    mid_range_nitrogen_pressure = 
            schreiner_equation__(&initial_inspired_n2_pressure, 
                                 &nitrogen_rate, 
                                 &mid_range_time, 
                                 &nitrogen_time_constant[*i - 1], 
                                 &initial_nitrogen_pressure[*i - 1]);
	    gas_tension_at_mid_range = 
            mid_range_helium_pressure + 
		    mid_range_nitrogen_pressure + 
		    constant_pressure_other_gases;
	    function_at_mid_range = 
            mid_range_ambient_pressure - 
		    gas_tension_at_mid_range - 
            gradient_onset_of_imperm;
	    if (function_at_mid_range <= 0.) {
	        time = mid_range_time;
	    }
	    if (fabs(differential_change) < .001 || 
		    function_at_mid_range == 0.) {
	        goto L100;
	    }
    } 
	printf("\nERROR! ROOT SEARCH EXCEEDED MAXIMUM ITERATIONS");
    pause();

/* =============================================================================== */
/*     When a solution with the desired accuracy is found, the program jumps out */
/*     of the loop to Line 100 and assigns the solution values for ambient */
/*     pressure and gas tension at the onset of impermeability. */
/* =============================================================================== */

L100:
    amb_pressure_onset_of_imperm[*i - 1] = mid_range_ambient_pressure;
    gas_tension_onset_of_imperm[*i - 1] = gas_tension_at_mid_range;
    return 0;
} /* onset_of_impermeability */

/* =============================================================================== */
/*     Varying Permeability Model (VPM) Decompression Program in FORTRAN */
/* 
/*     Author:  Erik C. Baker */
/* 
/*     "DISTRIBUTE FREELY - CREDIT THE AUTHORS" */
/* 
/*     This program extends the 1986 VPM algorithm (Yount & Hoffman) to include */
/*     mixed gas, repetitive, and altitude diving.  Developments to the algorithm */
/*     were made by David E. Yount, Eric B. Maiken, and Erik C. Baker over a */
/*     period from 1999 to 2001.  This work is dedicated in remembrance of */
/*     Professor David E. Yount who passed away on April 27, 2000. */
/* 
/*     Notes: */
/*     1.  This program uses the sixteen (16) half-time compartments of the */
/*         Buhlmann ZH-L16 model.  The optional Compartment 1b is used here with */
/*         half-times of 1.88 minutes for helium and 5.0 minutes for nitrogen. */
/* 
/*     2.  This program uses various DEC, IBM, and Microsoft extensions which */
/*         may not be supported by all FORTRAN compilers.  Comments are made with */
/*         a capital "C" in the first column or an exclamation point "!" placed */
/*         in a line after code.  An asterisk "*" in column 6 is a continuation */
/*         of the previous line.  All code, except for line numbers, starts in */
/*         column 7. */
/* 
/*     3.  Comments and suggestions for improvements are welcome.  Please */
/*         respond by e-mail to:  EBaker@se.aeieng.com */
/*
/*     Acknowledgment:  Thanks to Kurt Spaugh for recommendations on how to clean */
/*     up the code. */
/* ===============================================================================
/*     Converted to vpmdeco.c using f2c; R.McGinnis (CABER Swe) 5/01
/* =============================================================================== */

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <iostream.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include "support.h"

/* 'c' file I/O */

FILE    *fd;
int     err;

/* temp vars to simplify UNFMTLISTs */
real    fO2, fHe, fN2;
real    dc, rc, ssc;
integer mc;

/* Common Block Declarations */

real    water_vapor_pressure;
real    surface_tension_gamma, 
        skin_compression_gammac;
real    crit_volume_parameter_lambda;
real    minimum_deco_stop_time;
real    regeneration_time_constant;
real    constant_pressure_other_gases;
real    gradient_onset_of_imperm_atm;
real    run_time;
integer segment_number;
real    segment_time;
real    ending_ambient_pressure;
integer mix_number;
real    barometric_pressure;
logical units_equal_fsw, 
        units_equal_msw;
real    units_factor;
real    helium_time_constant[16];
real    nitrogen_time_constant[16];
real    helium_pressure[16], 
        nitrogen_pressure[16];
real    fraction_helium[10], 
        fraction_nitrogen[10];
real    initial_critical_radius_he[16], 
        initial_critical_radius_n2[16];
real    adjusted_critical_radius_he[16], 
        adjusted_critical_radius_n2[16];
real    max_crushing_pressure_he[16], 
        max_crushing_pressure_n2[16];
real    surface_phase_volume_time[16];
real    max_actual_gradient[16];
real    amb_pressure_onset_of_imperm[16], 
        gas_tension_onset_of_imperm[16];
real    initial_helium_pressure[16], 
        initial_nitrogen_pressure[16];
real    regenerated_radius_he[16], 
        regenerated_radius_n2[16];
real    adjusted_crushing_pressure_he[16], 
        adjusted_crushing_pressure_n2[16];
real    allowable_gradient_he[16], 
        allowable_gradient_n2[16];
real    initial_allowable_gradient_he[16], 
        initial_allowable_gradient_n2[16];

/* Forward references */

extern /* Subroutine */ int calc_barometric_pressure();
extern /* Subroutine */ int calc_deco_ceiling();
extern /* Subroutine */ int onset_of_impermeability();
extern /* Subroutine */ int radius_root_finder();
extern /* Subroutine */ int nuclear_regeneration(), clock_();
extern /* Subroutine */ int calc_deco_ceiling();
extern /* Subroutine */ int calc_barometric_pressure();
extern /* Subroutine */ int vpm_repetitive_algorithm();
extern /* Subroutine */ int gas_loadings_surface_interval();
extern /* Subroutine */ int critical_volume();	    ;
extern /* Subroutine */ int calc_crushing_pressure(), 
                            calc_surface_phase_volume_time(), 
                            decompression_stop();
extern /* Subroutine */ int calc_start_of_deco_zone();
extern /* Subroutine */ int calc_initial_allowable_gradient(), 
                            gas_loadings_ascent_descen();
extern /* Subroutine */ int gas_loadings_constant_depth();
extern /* Subroutine */ int vpm_altitude_dive_algorithm();
extern /* Subroutine */ int calc_max_actual_gradient(), 
                            projected_ascent();

extern /* Subroutine */ doublereal schreiner_equation__();
extern /* Subroutine */ doublereal haldane_equation__();

/* =============================================================================== */
/*     Begin MAIN proc (see endofpgm for alternate entry) */
/* =============================================================================== */

int main()
{

/* =============================================================================== */
/*     ASSIGN HALF-TIME VALUES TO BUHLMANN COMPARTMENT ARRAYS */
/* =============================================================================== */
    /* Initialized data */

    real helium_half_time[16] = 
        {  1.88,   3.02,   4.72,   6.99,
          10.21,  14.48,  20.53,  29.11,
          41.2,   55.19,  70.69,  90.34,
         115.29, 147.42, 188.24, 240.03 };
    real nitrogen_half_time[16] = 
        { 5.,    8.,  12.5, 18.5,
          27.,  38.3, 54.3, 77.,
          109.,146., 187.,  239.,
          305.,390., 498.,  635. };

    /* Format strings */
    char fmt_900[] = "\n";
    // was "(\002 \002)";
    char fmt_901[] = "ERROR! UNITS MUST BE FSW OR MSW\n";
    // was: "(\0020ERROR! UNITS MUST BE FSW OR MSW\002)";
    char fmt_902[] = "ERROR! ALTITUDE DIVE ALGORITHM MUST BE ON OR OFF\n";
    // was "(\0020ERROR! ALTITUDE DIVE ALGORITHM MUST BE ON OR OFF\002)";
    char fmt_903[] = "ERROR! RADIUS MUST BE BETWEEN 0.2 AND 1.35 MICRONS\n";
    // was "(\0020ERROR! RADIUS MUST BE BETWEEN 0.2 AND 1.35 MICRONS\002)";
    char fmt_904[] = "ERROR! CRITICAL VOLUME ALGORITHM MUST BE ON OR OFF";
    // was "(\0020ERROR! CRITICAL VOLUME ALGORITHM MUST BE ON OR OFF\002)";
    char fmt_800[] = "UNITS = FEET OF SEAWATER (FSW)\n";
    // was "(\0020UNITS = FEET OF SEAWATER (FSW)\002)";
    char fmt_801[] = "UNITS = METERS OF SEAWATER (MSW)\n";
    // was "(\0020UNITS = METERS OF SEAWATER (MSW)\002)";
    char fmt_802[] = "ALTITUDE =  %7.1lf\nBAROMETRIC PRESSURE = %6.3lf\n\n";
    // was "(\0020ALTITUDE = \002,1x,f7.1,4x,\002BAROMETRIC PRESSURE = \002,f6.3)";
    char fmt_805[] = "%s\n";
    // was "(a70)";
    char fmt_810[] = "(\033E\033&a10L\033&l80F\033&l8D\033(s0p16.67h8.5\002)"; /* tbd */
    // was "(\002\033E\033&a10L\033&l80F\033&l8D\033(s0p16.67h8.5\002)";
    char fmt_811[] = "                          DECOMPRESSION CALCULATION PROGRAM\n";
    // was "(26x,\002DECOMPRESSION CALCULATION PROGRAM\002)";
    char fmt_812[] = "                        Developed in FORTRAN by Erik C. Baker\n";
    // was "(24x,\002Developed in FORTRAN by Erik C. Baker\002)";
    char fmt_813[] = "\n";
    // was "(\002 \002)";
    char fmt_814[] = "Program Run:    %2.2d-%2.2d-%4d at %2.2d:%2.2d %1sm                       Model: VPM 2001\n";
    // was "(\002Program Run:\002,4x,i2.2,\002-\002,i2.2,\002-\002,i4,1x,\002at\002,1x,i2.2,\002:\002,i2.2,1x,a1,\002m\002,23x,\002Model: VPM 2001\002)";
    char fmt_815[] = "Description:    %70s\n";
    // was "(\002Description:\002,4x,a70)";
    char fmt_906[] = "ERROR IN INPUT FILE (GASMIX DATA)\n";
    // was "(\0020ERROR IN INPUT FILE (GASMIX DATA)\002)";
    char fmt_820[] = "Gasmix Summary:                        FO2    FHe    FN2\n";
    // was "(\002Gasmix Summary:\002,24x,\002FO2\002,4x,\002FHe\002,4x,\002FN2\002)";
    char fmt_821[] = "                          Gasmix #%2d  %5.3lf  %5.3lf  %5.3lf\n";
    // was "(26x,\002Gasmix #\002,i2,2x,f5.3,2x,f5.3,2x,f5.3)";
    char fmt_830[] = "                                    DIVE PROFILE\n";
    // was "(36x,\002DIVE PROFILE\002)";
    char fmt_831[] = "Seg-  Segm.  Run   | Gasmix | Ascent    From     To      Rate    | Constant\n";
    // was "(\002Seg-\002,2x,\002Segm.\002,2x,\002Run\002,3x,\002|\002,1x,\002Gasmix\002,1x,\002|\002,1x,\002Ascent\002,4x,\002From\002,5x,\002To\002,6x,\002Rate\002,4x,\002|\002,1x,\002Constant\002)";
    char fmt_832[] = "ment  Time   Time  |  Used  |   or     Depth   Depth    +Dn/-Up  |  Depth\n";
    // was "(\002ment\002,2x,\002Time\002,3x,\002Time\002,2x,\002|\002,2x,\002Used\002,2x,\002|\002,3x,\002or\002,5x,\002Depth\002,3x,\002Depth\002,4x,\002+Dn/-Up\002,2x,\002|\002,2x,\002Depth\002)";
    char fmt_833[] = "  #   (min)  (min) |    #   | Descent  (%4s)  (%4s)  (%7s) |  (%4s)\n";
    // was "(2x,\002#\002,3x,\002(min)\002,2x,\002(min)\002,1x,\002|\002,4x,\002#\002,3x,\002|\002,1x,\002Descent\002,2x,\002(\002,a4,\002)\002,2x,\002(\002,a4,\002)\002,2x,\002(\002,a7,\002)\002,1x,\002|\002,2x,\002(\002,a4,\002)\002)";
    char fmt_834[] = "----- -----  ----- | ------ | -------  ------  ------  --------- | --------\n";
    // was "(\002-----\002,1x,\002-----\002,2x,\002-----\002,1x,\002|\002,1x,\002------\002,1x,\002|\002,1x,\002-------\002,2x,\002------\002,2x,\002------\002,2x,\002---------\002,1x,\002|\002,1x,\002--------\002)";
    char fmt_840[] = "%3d   %5.1lf %6.1lf |   %2d   | %7s%7.0lf %7.0lf   %7.1lf   |\n";
    // was "(i3,3x,f5.1,1x,f6.1,1x,\002|\002,3x,i2,3x,\002|\002,1x,a7,f7.0,1x,f7.0,3x,f7.1,3x,\002|\002)";
    char fmt_845[] = "%3d   %5.1lf %6.1lf |   %2d   |                                    | %7.0lf\n";
    // was "(i3,3x,f5.1,1x,f6.1,1x,\002|\002,3x,i2,3x,\002|\002,36x,\002|\002,f7.0)";
    char fmt_907[] = "ERROR IN INPUT FILE (PROFILE CODE)\n";
    // was "(\0020ERROR IN INPUT FILE (PROFILE CODE)\002)";
    char fmt_850[] = "                               DECOMPRESSION PROFILE\n";
    // was "(31x,\002DECOMPRESSION PROFILE\002)";
    char fmt_857[] = "          Leading compartment enters the decompression zone at%7.1lf %4s\n";
    // was "(10x,\002Leading compartment enters the decompression zone\002,1x,\002at\002,f7.1,1x,a4\n";
    char fmt_858[] = "                 Deepest possible decompression stop is%7.1lf %4s\n";
    // was "(17x,\002Deepest possible decompression stop is\002,f7.1,1x,a4)";
    char fmt_851[] = "Seg-  Segm.  Run   | Gasmix | Ascent   Ascent   Col   |  DECO   STOP   RUN\n";
    // was "(\002Seg-\002,2x,\002Segm.\002,2x,\002Run\002,3x,\002|\002,1x,\002Gasmix\002,1x,\002|\002,1x,\002Ascent\002,3x,\002Ascent\002,3x,\002Col\002,3x,\002|\002,2x,\002DECO\002,3x,\002STOP\002,3x,\002RUN\002)";
    char fmt_852[] = "ment  Time   Time  |  Used  |   To      Rate    Not   |  STOP   TIME   TIME\n";
    // was "(\002ment\002,2x,\002Time\002,3x,\002Time\002,2x,\002|\002,2x,\002Used\002,2x,\002|\002,3x,\002To\002,6x,\002Rate\002,4x,\002Not\002,3x,\002|\002,2x,\002STOP\002,3x,\002TIME\002,3x,\002TIME\002)";
    char fmt_853[] = "  #   (min)  (min) |    #   | (%4s) (%7s)  Used  | (%4s)  (min)  (min)\n";
    // was "(2x,\002#\002,3x,\002(min)\002,2x,\002(min)\002,1x,\002|\002,4x,\002#\002,3x,\002|\002,1x,\002(\002,a4,\002)\002,1x,\002(\002,a7,\002)\002,2x,\002Used\002,2x,\002|\002,1x,\002(\002,a4,\002)\002,2x,\002(min)\002,2x,\002(min)\002)";
    char fmt_854[] = "----- -----  ----- | ------ | ------ --------- ------ | ------ -----  -----\n";
    // was "(\002-----\002,1x,\002-----\002,2x,\002-----\002,1x,\002|\002,1x,\002------\002,1x,\002|\002,1x,\002------\002,1x,\002---------\002,1x,\002------\002,1x,\002|\002,1x,\002------\002,2x,\002-----\002,2x,\002-----\002)";
    char fmt_905[] = "ERROR! STEP SIZE IS TOO LARGE TO DECOMPRESS\n";
    // was "(\0020ERROR! STEP SIZE IS TOO LARGE TO DECOMPRESS\002)";
    char fmt_860[] = "%3d   %5.1lf %6.1lf |   %2d   |  %4.0lf   %6.1lf          |\n";
    // was "(i3,3x,f5.1,1x,f6.1,1x,\002|\002,3x,i2,3x,\002|\002,2x,f4.0,3x,f6.1,10x,\002|\002)";
    char fmt_862[] = "%3d   %5.1lf %6.1lf |   %2d   |                         |  %4d   %4d  %5d\n";
    // was "(i3,3x,f5.1,1x,f6.1,1x,\002|\002,3x,i2,3x,\002|\002,25x,\002|\002,2x,i4,3x,i4,2x,i5)";
    char fmt_863[] = "%3d   %5.1lf %6.1lf |   %2d   |                         |  %5.0lf %6.1lf %7.1lf\n";
    // was "(i3,3x,f5.1,1x,f6.1,1x,\002|\002,3x,i2,3x,\002|\002,25x,\002|\002,2x,f5.0,1x,f6.1,1x,f7.1)";
    char fmt_880[] = "\n";
    // was "(\002 \002)";
    char fmt_890[] = "REPETITIVE DIVE:\n";
    // was "(\002REPETITIVE DIVE:\002)";
    char fmt_908[] = "ERROR IN INPUT FILE (REPETITIVE DIVE CODE)\n";
    // was "(\0020ERROR IN INPUT FILE (REPETITIVE DIVE CODE)\002)";
    char fmt_871[] = " PROGRAM CALCULATIONS COMPLETE\n";
    // was "(\002 PROGRAM CALCULATIONS COMPLETE\002)";
    char fmt_872[] = "Output data is located in the file VPMDECO.OUT\n";
    // was "(\0020Output data is located in the file VPMDECO.OUT\002)";

    /* System generated locals */
    integer i1;
    real    r1;

    /* Local variables */
    real    fraction_oxygen[10];
    real    run_time_end_of_segment, 
            rate, 
	        n2_pressure_start_of_ascent[16];
    shortint year;
    char    altitude_dive_algorithm[4];                         /* allow null stopper */
    real    depth_change[10];
    char    word[8];                                            /* allow null stopper */
    logical altitude_dive_algorithm_off;
    real    ending_depth;
    integer profile_code;
    real    run_time_start_of_deco_zone;
    char    line1[71];                                          /* allow null stopper */
    real    deepest_possible_stop_depth, 
            he_pressure_start_of_ascent[16], 
            altitude_of_dive, 
            step_size_change[10];
    integer i, 
            j, 
            repetitive_dive_flag;
    char    m[2];                                               /* allow null stopper */
    real    sum_of_fractions;
    real    first_stop_depth;
    real    depth, depth_start_of_deco_zone;
    real    sum_check;
    shortint month;
    real    run_time_start_of_ascent;
    char    units[4];                                           /* allow null stopper */
    integer number_of_changes;
    real    stop_time;
    real    step_size, 
            next_stop, 
            last_run_time, 
	        phase_volume_time[16], 
            surface_interval_time;
    integer mix_change[10];

    shortint minute;
    char    critical_volume_algorithm[4];                       /* allow null stopper */                
    real    deco_ceiling_depth;
    shortint clock_hour;
    logical critical_volume_algorithm_off;
    real    pressure_other_gases_mmhg;
    logical schedule_converged;
    real    critical_radius_n2_microns, 
            starting_depth, 
	        deco_phase_volume_time;
    shortint day;
    real    rate_change[10], 
            critical_radius_he_microns, 
	        last_phase_volume_time[16], 
            n2_pressure_start_of_deco_zone[16];
    char    units_word1[5], units_word2[8];                     /* allow null stopper */
    real    critical_volume_comparison, deco_stop_depth;
    integer segment_number_start_of_ascent;
    real    rounding_operation, 
            rounding_operation2, 
	        he_pressure_start_of_deco_zone[16];
    integer number_of_mixes;

/* =============================================================================== */
/*     NAMELIST FOR PROGRAM SETTINGS (READ IN FROM ASCII TEXT FILE) */
/* =============================================================================== */
    /* NAMELIST stuff */

    VALDESC altitude_dive_algorithm_dv =    
        { "ALTITUDE_DIVE_ALGORITHM", altitude_dive_algorithm, CHARTYPE(altitude_dive_algorithm) };
    VALDESC skin_compression_gammac_dv =    
        { "SKIN_COMPRESSION_GAMMAC", (char *)&skin_compression_gammac, REALTYPE };
    VALDESC crit_volume_parameter_lambda_dv = 
        { "CRIT_VOLUME_PARAMETER_LAMBDA", (char *)&crit_volume_parameter_lambda, REALTYPE };
    VALDESC gradient_onset_of_imperm_atm_dv = 
        { "GRADIENT_ONSET_OF_IMPERM_ATM", (char *)&gradient_onset_of_imperm_atm, REALTYPE };
    VALDESC units_dv =                      
        { "UNITS", units, CHARTYPE(units) };
    VALDESC surface_tension_gamma_dv =      
        { "SURFACE_TENSION_GAMMA", (char *)&surface_tension_gamma, REALTYPE };
    VALDESC critical_volume_algorithm_dv =  
        { "CRITICAL_VOLUME_ALGORITHM", critical_volume_algorithm, CHARTYPE(critical_volume_algorithm) };
    VALDESC pressure_other_gases_mmhg_dv =  
        { "PRESSURE_OTHER_GASES_MMHG", (char *)&pressure_other_gases_mmhg, REALTYPE };
    VALDESC critical_radius_n2_microns_dv = 
        { "CRITICAL_RADIUS_N2_MICRONS", (char *)&critical_radius_n2_microns, REALTYPE };
    VALDESC critical_radius_he_microns_dv =     
        { "CRITICAL_RADIUS_HE_MICRONS", (char *)&critical_radius_he_microns, REALTYPE };
    VALDESC minimum_deco_stop_time_dv =     
        { "MINIMUM_DECO_STOP_TIME", (char *)&minimum_deco_stop_time, REALTYPE };
    VALDESC regeneration_time_constant_dv = 
        { "REGENERATION_TIME_CONSTANT", (char *)&regeneration_time_constant, REALTYPE };

    VALDESC *program_settings_vl[] = { 
        &units_dv, 
        &altitude_dive_algorithm_dv, 
        &minimum_deco_stop_time_dv, 
        &critical_radius_n2_microns_dv, 
        &critical_radius_he_microns_dv,
	    &critical_volume_algorithm_dv, 
        &crit_volume_parameter_lambda_dv, 
        &gradient_onset_of_imperm_atm_dv, 
        &surface_tension_gamma_dv, 
        &skin_compression_gammac_dv, 
        &regeneration_time_constant_dv, 
        &pressure_other_gases_mmhg_dv };

    NAMELIST program_settings = { 
        "PROGRAM_SETTINGS", 
	    program_settings_vl, 
        sizeof(program_settings_vl)/sizeof(VALDESC*)};

/* =============================================================================== */
/*     UNFMTLISTS FOR INPUT (ASCII TEXT FILE) */
/* =============================================================================== */
    /* UNFMTLIST stuff */

    
    VALDESC line1_dv =    
        {"", line1, CHARTYPE(line1) };
    VALDESC *line1_vl[] = {&line1_dv}; 
    UNFMTLIST line1_ul =  { 
	    line1_vl, 
        sizeof(line1_vl)/sizeof(VALDESC*)};

    VALDESC line2_dv =    {"", &number_of_mixes, INTTYPE };
    VALDESC *line2_vl[] = {&line2_dv}; 
    UNFMTLIST line2_ul =  { 
	    line2_vl, 
        sizeof(line2_vl)/sizeof(VALDESC*)};

    VALDESC fO2_dv =                            /* note: uses temp vars */
        {"", &fO2, REALTYPE };
    VALDESC fHe_dv =    
        {"", &fHe, REALTYPE };
    VALDESC fN2_dv =  
        {"", &fN2, REALTYPE };
    VALDESC *line3_vl[] = {&fO2_dv,
                           &fHe_dv,
                           &fN2_dv}; 
    UNFMTLIST line3_ul = { 
	    line3_vl, 
        sizeof(line3_vl)/sizeof(VALDESC*)};

    VALDESC line4_dv =    
        {"", &profile_code, INTTYPE };
    VALDESC *line4_vl[] = {&line4_dv}; 
    UNFMTLIST line4_ul =  { 
        line4_vl, 
        sizeof(line4_vl)/sizeof(VALDESC*)};

    VALDESC starting_depth_dv = 
        {"", &starting_depth, REALTYPE };
    VALDESC ending_depth_dv =   
        {"", &ending_depth, REALTYPE };
    VALDESC rate_dv =           
        {"", &rate, REALTYPE };
    VALDESC mix_number_dv =     
        {"", &mix_number, INTTYPE };
    VALDESC *line5_vl[] = {&starting_depth_dv,
                           &ending_depth_dv,
                           &rate_dv,
                           &mix_number_dv}; 
    UNFMTLIST line5_ul = { 
	    line5_vl, 
        sizeof(line5_vl)/sizeof(VALDESC*)};

    VALDESC depth_dv =                      
        {"", &depth, REALTYPE };
    VALDESC run_time_end_of_segment_dv =    
        {"", &run_time_end_of_segment, REALTYPE };
    VALDESC *line6_vl[] = {&depth_dv,
                           &run_time_end_of_segment_dv,
                           &mix_number_dv}; 
    UNFMTLIST line6_ul = { line6_vl, 
        sizeof(line6_vl)/sizeof(VALDESC*)};

    VALDESC line7_dv =    
        {"", &number_of_changes, INTTYPE };
    VALDESC *line7_vl[] = 
        {&line7_dv}; 
    UNFMTLIST line7_ul =  { line7_vl, 
        sizeof(line7_vl)/sizeof(VALDESC*)};

    VALDESC dc_dv =                             /* note: uses temp vars */
        {"", &dc, REALTYPE };
    VALDESC mc_dv =         
        {"", &mc, INTTYPE };
    VALDESC rc_dv =        
        {"", &rc, REALTYPE };
    VALDESC ssc_dv =   
        {"", &ssc, REALTYPE };
    VALDESC *line8_vl[] = {&dc_dv,
                           &mc_dv,
                           &rc_dv,
                           &ssc_dv}; 
    UNFMTLIST line8_ul = { 
	    line8_vl, 
        sizeof(line8_vl)/sizeof(VALDESC*)};

    VALDESC line9_dv =    
        {"", &repetitive_dive_flag, INTTYPE };
    VALDESC *line9_vl[] = {&line9_dv}; 
    UNFMTLIST line9_ul =  { 
	    line9_vl, 
        sizeof(line9_vl)/sizeof(VALDESC*)};

    VALDESC line10_dv =    
        {"", &surface_interval_time, REALTYPE };
    VALDESC *line10_vl[] = {&line10_dv}; 
    UNFMTLIST line10_ul =  { 
	    line10_vl, 
        sizeof(line10_vl)/sizeof(VALDESC*)};

/* loop */
/* =============================================================================== */
/*     BEGIN PROGRAM EXECUTION WITH OUTPUT MESSAGE TO SCREEN */
/* =============================================================================== */
    printf("\t\n");                     /* Pass "clear screen" to MS operating sys */
    printf("\n");
    printf("PROGRAM VPMDECO\n");
    printf("\n");

/* =============================================================================== */
/*     READ IN PROGRAM SETTINGS AND CHECK FOR ERRORS */
/*     IF THERE ARE ERRORS, WRITE AN ERROR MESSAGE AND TERMINATE PROGRAM */
/* =============================================================================== */

    /* open a file for reading (LUN7 in the fortran program) */
    fd = fopen("VPMDECO.IN", "rt");
    if (!fd) {
        printf("\nCannot open VPMDECO.IN");
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }

    // READ (10,Program_Settings)
    if (err = ReadNameList("VPMDECO.SET", &program_settings)) {
        printf("\nError %d returned reading VPMDECO.SET", err);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }

    if (strcmp(units, "fsw") == 0 || strcmp(units, "FSW") == 0) {
	    units_equal_fsw = TRUE_;
	    units_equal_msw = FALSE_;
    } else if (strcmp(units, "msw") == 0 || strcmp(units, "MSW") == 0) {
	    units_equal_fsw = FALSE_;
	    units_equal_msw = TRUE_;
    } else {
	    printf("\t\n");
        printf(fmt_901);
        printf(fmt_900);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }
    if (strcmp(altitude_dive_algorithm, "ON") == 0 || 
	    strcmp(altitude_dive_algorithm, "on") == 0)
	     {
	    altitude_dive_algorithm_off = FALSE_;
    } else if (strcmp(altitude_dive_algorithm, "OFF") == 0 || 
               strcmp(altitude_dive_algorithm, "off") == 0) {
	    altitude_dive_algorithm_off = TRUE_;
    } else {
	    printf("\t\n");
        printf(fmt_902);
        printf(fmt_900);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }
    if (critical_radius_n2_microns < .2 || 
	    critical_radius_n2_microns > 1.35) {
	    printf("\t\n");
        printf(fmt_903);
        printf(fmt_900);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }
    if (critical_radius_he_microns < .2 || 
	    critical_radius_he_microns > 1.35) {
	    printf("\t\n");
        printf(fmt_903);
        printf(fmt_900);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }
    if (strcmp(critical_volume_algorithm, "ON") == 0 ||
	        strcmp(critical_volume_algorithm, "on") == 0) {
	    critical_volume_algorithm_off = FALSE_;
    } else if (strcmp(critical_volume_algorithm, "OFF") == 0 || 
	           strcmp(critical_volume_algorithm, "off") == 0) {
	    critical_volume_algorithm_off = TRUE_;
    } else {
	    printf("\t\n");
        printf(fmt_904);
        printf(fmt_900);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }

/* =============================================================================== */
/*     INITIALIZE CONSTANTS/VARIABLES BASED ON SELECTION OF UNITS - FSW OR MSW */
/*     fsw = feet of seawater, a unit of pressure */
/*     msw = meters of seawater, a unit of pressure */
/* =============================================================================== */

    if (units_equal_fsw) {
    	printf(fmt_800);
	    strcpy(units_word1, "fswg");
	    strcpy(units_word2, "fsw/min");
	    units_factor = 33.;
	    water_vapor_pressure = 1.607;/* (Schreiner value)  based on respiratory quotient */
    }
    if (units_equal_msw) {
	    printf(fmt_801);
	    strcpy(units_word1, "mswg");
	    strcpy(units_word2, "msw/min");
	    units_factor = 10.1325;
	    water_vapor_pressure = .493;/* (Schreiner value)  based on respiratory quotien */
    }

/* =============================================================================== */
/*     INITIALIZE CONSTANTS/VARIABLES */
/* =============================================================================== */

    constant_pressure_other_gases = 
        pressure_other_gases_mmhg / 760. * units_factor;
    run_time = 0.;
    segment_number = 0;
    for (i = 1; i <= 16; ++i) {
	    helium_time_constant[i - 1] = log(2.) / helium_half_time[i - 1];
	    nitrogen_time_constant[i - 1] = log(2.) / nitrogen_half_time[i - 1];
	    max_crushing_pressure_he[i - 1] = 0.;
	    max_crushing_pressure_n2[i - 1] = 0.;
	    max_actual_gradient[i - 1] = 0.;
	    surface_phase_volume_time[i - 1] = 0.;
	    amb_pressure_onset_of_imperm[i - 1] = 0.;
	    gas_tension_onset_of_imperm[i - 1] = 0.;
	    initial_critical_radius_n2[i - 1] = 
	    	critical_radius_n2_microns * 1e-6;
	    initial_critical_radius_he[i - 1] = 
		    critical_radius_he_microns * 1e-6;
    }

/* =============================================================================== */
/*     INITIALIZE VARIABLES FOR SEA LEVEL OR ALTITUDE DIVE */
/*     See subroutines for explanation of altitude calculations.  Purposes are */
/*     1) to determine barometric pressure and 2) set or adjust the VPM critical */
/*     radius variables and gas loadings, as applicable, based on altitude, */
/*     ascent to altitude before the dive, and time at altitude before the dive */
/* =============================================================================== */

    if (altitude_dive_algorithm_off) {
	    altitude_of_dive = 0.;
	    calc_barometric_pressure(&altitude_of_dive); 
	    printf(fmt_802, altitude_of_dive, 
                        barometric_pressure);
	    for (i = 1; i <= 16; ++i) {
	        adjusted_critical_radius_n2[i - 1] = 
		        initial_critical_radius_n2[i - 1];
	        adjusted_critical_radius_he[i - 1] = 
		        initial_critical_radius_he[i - 1];
	        helium_pressure[i - 1] = 0.;
	        nitrogen_pressure[i - 1] = 
		        (barometric_pressure - water_vapor_pressure) * ..79;
	    }
    } else {
	    vpm_altitude_dive_algorithm();  
    }

/* =============================================================================== */
/*     START OF REPETITIVE DIVE LOOP */
/*     This is the largest loop in the main program and operates between Lines */
/*     30 and 330.  If there is one or more repetitive dives, the program will */
/*     return to this point to process each repetitive dive. */
/* =============================================================================== */
/* L30: */
    while(TRUE_) {                      /* until there is an break statement */
                                        /* loop will run continuous */

/* =============================================================================== */
/*     INPUT DIVE DESCRIPTION AND GAS MIX DATA FROM ASCII TEXT INPUT FILE */
/*     BEGIN WRITING HEADINGS/OUTPUT TO ASCII TEXT OUTPUT FILE */
/*     See separate explanation of format for input file. */
/* =============================================================================== */

	    // READ (7,805) Line1  
        if (err = ReadUnformatted(fd, &line1_ul)) {
            printf("\nError %d returned reading VPMDECO.IN", err);
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
        }       
	    clock_(&year, &month, &day, &clock_hour, &minute, m); 
//	    printf(fmt_810);        /* tbd */
        printf(fmt_811);
        printf(fmt_812);
        printf(fmt_813);
        printf(fmt_814, month, 
                        day, 
                        year, 
                        clock_hour, 
                        minute, 
                        m);
	    printf(fmt_813);
        printf(fmt_815, line1);
	    printf(fmt_813);

	    // READ (7,*) Number_of_Mixes                   !check for errors in gasmixes 
        if (err = ReadUnformatted(fd, &line2_ul)) {
            printf("\nError %d returned reading VPMDECO.IN", err);
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
        }       
	    i1 = number_of_mixes;
	    for (i = 1; i <= i1; ++i) {
	        // READ (7,*) Fraction_Oxygen(I), Fraction_Helium(I), Fraction_Nitrogen(I) 
            if (err = ReadUnformatted(fd, &line3_ul)) {
                printf("\nError %d returned reading VPMDECO.IN", err);
                printf("\nPROGRAM TERMINATED\n");
                exit(1);
            }  
            fraction_oxygen[i - 1] = fO2;
            fraction_helium[i - 1] = fHe;
            fraction_nitrogen[i - 1] = fN2;

	        sum_of_fractions = 
                fraction_oxygen[i - 1] + 
		        fraction_helium[i - 1] + 
		        fraction_nitrogen[i - 1];
	        sum_check = sum_of_fractions;
	        if (sum_check != 1.) {
	            printf("\t\n");
                printf(fmt_906);
                printf(fmt_900);
                printf("\nPROGRAM TERMINATED\n");
                exit(1);
	        }
	    }
	    printf(fmt_820);
	    i1 = number_of_mixes;
	    for (j = 1; j <= i1; ++j) {
	        printf(fmt_821, j, 
                            fraction_oxygen[j - 1], 
                            fraction_helium[j - 1], 
                            fraction_nitrogen[j - 1]);
	    }
        printf(fmt_813);
        printf(fmt_813);
        printf(fmt_830);
        printf(fmt_813);
        printf(fmt_831);
        printf(fmt_832);
        printf(fmt_833, units_word1, 
                        units_word1, 
                        units_word2, 
                        units_word1);
        printf(fmt_834);

/* =============================================================================== */
/*     DIVE PROFILE LOOP - INPUT DIVE PROFILE DATA FROM ASCII TEXT INPUT FILE */
/*     AND PROCESS DIVE AS A SERIES OF ASCENT/DESCENT AND CONSTANT DEPTH */
/*     SEGMENTS.  THIS ALLOWS FOR MULTI-LEVEL DIVES AND UNUSUAL PROFILES.  UPDATE */
/*     GAS LOADINGS FOR EACH SEGMENT.  IF IT IS A DESCENT SEGMENT, CALC CRUSHING */
/*     PRESSURE ON CRITICAL RADII IN EACH COMPARTMENT. */
/*     "Instantaneous" descents are not used in the VPM.  All ascent/descent */
/*     segments must have a realistic rate of ascent/descent.  Unlike Haldanian */
/*     models, the VPM is actually more conservative when the descent rate is */
/*     slower becuase the effective crushing pressure is reduced.  Also, a */
/*     realistic actual supersaturation gradient must be calculated during */
/*     ascents as this affects critical radii adjustments for repetitive dives. */
/*     Profile codes: 1 = Ascent/Descent, 2 = Constant Depth, 99 = Decompress */
/* =============================================================================== */

	while(TRUE_) {          /* until there is an exit statement loop will run continuous */
	    // READ (7,*) Profile_Code 
        if (err = ReadUnformatted(fd, &line4_ul)) {
            printf("\nError %d returned reading VPMDECO.IN", err);
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
        }       
	    if (profile_code == 1) {                        /* ascent or descent */
		    // READ (7,*) Starting_Depth, Ending_Depth, Rate, Mix_Number  
            if (err = ReadUnformatted(fd, &line5_ul)) {
                printf("\nError %d returned reading VPMDECO.IN", err);
                printf("\nPROGRAM TERMINATED\n");
                exit(1);
            }       
		    gas_loadings_ascent_descen(&starting_depth, 
                                       &ending_depth, 
                                       &rate);           
		    if (ending_depth > starting_depth) {
		        calc_crushing_pressure(&starting_depth, 
                                    &ending_depth, &rate);       
		    }
		    if (ending_depth > starting_depth) {
		        strcpy(word, "Descent");
		    } else if (starting_depth > ending_depth) {
		        strcpy(word, "Ascent ");
		    } else {
		        strcpy(word, "ERROR");
		    }
		    printf(fmt_840, segment_number, 
                            segment_time, 
                            run_time, 
                            mix_number, 
                            word,
                            starting_depth, 
                            ending_depth, 
                            rate);
	    } else if (profile_code == 2) {                 /* constant depth */
		    // READ (7,*) Depth, Run_Time_End_of_Segment, Mix_Number  
            if (err = ReadUnformatted(fd, &line6_ul)) {
                printf("\nError %d returned reading VPMDECO.IN", err);
                printf("\nPROGRAM TERMINATED\n");
                exit(1);
            }       
		    gas_loadings_constant_depth(&depth, 
                                        &run_time_end_of_segment);       
		    printf(fmt_845, segment_number, 
                            segment_time, 
                            run_time, 
                            mix_number, 
                            depth);
	    } else if (profile_code == 99) {                /* begin decompression */
	    	break;                                      /* .. get out of while loop */
	    } else {                                        /* unknown profile_code */
	        printf("\t\n");
            printf(fmt_907);
            printf(fmt_900);
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
	    }
	}

/* =============================================================================== */
/*     BEGIN PROCESS OF ASCENT AND DECOMPRESSION */
/*     First, calculate the regeneration of critical radii that takes place over */
/*     the dive time.  The regeneration time constant has a time scale of weeks */
/*     so this will have very little impact on dives of normal length, but will */
/*     have major impact for saturation dives. */
/* =============================================================================== */

	nuclear_regeneration(&run_time);

/* =============================================================================== */
/*     CALCULATE INITIAL ALLOWABLE GRADIENTS FOR ASCENT */
/*     This is based on the maximum effective crushing pressure on critical radii */
/*     in each compartment achieved during the dive profile. */
/* =============================================================================== */

	calc_initial_allowable_gradient();

/* =============================================================================== */
/*     SAVE VARIABLES AT START OF ASCENT (END OF BOTTOM TIME) SINCE THESE WILL */
/*     BE USED LATER TO COMPUTE THE FINAL ASCENT PROFILE THAT IS WRITTEN TO THE */
/*     OUTPUT FILE. */
/*     The VPM uses an iterative process to compute decompression schedules so */
/*     there will be more than one pass through the decompression loop. */
/* =============================================================================== */

	for (i = 1; i <= 16; ++i) {
	    he_pressure_start_of_ascent[i - 1] = 
		    helium_pressure[i - 1];
	    n2_pressure_start_of_ascent[i - 1] = 
		    nitrogen_pressure[i - 1];
	}
	run_time_start_of_ascent = run_time;
	segment_number_start_of_ascent = segment_number;

/* =============================================================================== */
/*     INPUT PARAMETERS TO BE USED FOR STAGED DECOMPRESSION AND SAVE IN ARRAYS. */
/*     ASSIGN INITAL PARAMETERS TO BE USED AT START OF ASCENT */
/*     The user has the ability to change mix, ascent rate, and step size in any */
/*     combination at any depth during the ascent. */
/* =============================================================================== */

	//  READ (7,*) Number_of_Changes
    if (err = ReadUnformatted(fd, &line7_ul)) {
        printf("\nError %d returned reading VPMDECO.IN", err);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }       
	i1 = number_of_changes;
	for (i = 1; i <= i1; ++i) {
	    // READ (7,*) Depth_Change(I), Mix_Change(I), Rate_Change(I), Step_Size_Change(I)
        if (err = ReadUnformatted(fd, &line8_ul)) {
            printf("\nError %d returned reading VPMDECO.IN", err);
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
        }              
	    depth_change[i - 1] = dc;
	    mix_change[i - 1] = mc;
	    rate_change[i - 1] = rc;
	    step_size_change[i - 1] = ssc;
	}
	starting_depth = depth_change[0];
	mix_number = mix_change[0];
	rate = rate_change[0];
	step_size = step_size_change[0];

/* =============================================================================== */
/*     CALCULATE THE DEPTH WHERE THE DECOMPRESSION ZONE BEGINS FOR THIS PROFILE */
/*     BASED ON THE INITIAL ASCENT PARAMETERS AND WRITE THE DEEPEST POSSIBLE */
/*     DECOMPRESSION STOP DEPTH TO THE OUTPUT FILE */
/*     Knowing where the decompression zone starts is very important.  Below */
/*     that depth there is no possibility for bubble formation because there */
/*     will be no supersaturation gradients.  Deco stops should never start */
/*     below the deco zone.  The deepest possible stop deco stop depth is */
/*     defined as the next "standard" stop depth above the point where the */
/*     leading compartment enters the deco zone.  Thus, the program will not */
/*     base this calculation on step sizes larger than 10 fsw or 3 msw.  The */
/*     deepest possible stop depth is not used in the program, per se, rather */
/*     it is information to tell the diver where to start putting on the brakes */
/*     during ascent.  This should be prominently displayed by any deco program. */
/* =============================================================================== */

	calc_start_of_deco_zone(&starting_depth, &rate, &depth_start_of_deco_zone);  
	if (units_equal_fsw) {
	    if (step_size < 10.) {
		    rounding_operation = depth_start_of_deco_zone / step_size - .5;
		    deepest_possible_stop_depth = r_nint(&rounding_operation) * step_size;
	    } else {
		    rounding_operation = depth_start_of_deco_zone /  10. - .5;
		    deepest_possible_stop_depth = r_nint(&rounding_operation) * 10.;
	    }
	}
	if (units_equal_msw) {
	    if (step_size < 3.) {
		    rounding_operation = depth_start_of_deco_zone / step_size - .5;
		    deepest_possible_stop_depth = r_nint(&rounding_operation) * step_size;
	    } else {
		    rounding_operation = depth_start_of_deco_zone /  3. - .5;
		    deepest_possible_stop_depth = r_nint(&rounding_operation) * 3.;
	    }
	}
    printf(fmt_813);
    printf(fmt_813);
    printf(fmt_850);
    printf(fmt_813);
    printf(fmt_857, depth_start_of_deco_zone, 
                    units_word1);
    printf(fmt_858, deepest_possible_stop_depth, 
                    units_word1);
    printf(fmt_813);
    printf(fmt_851);
    printf(fmt_852);
    printf(fmt_853, units_word1, 
                    units_word2, 
                    units_word1);
    printf(fmt_854);

/* =============================================================================== */
/*     TEMPORARILY ASCEND PROFILE TO THE START OF THE DECOMPRESSION ZONE, SAVE */
/*     VARIABLES AT THIS POINT, AND INITIALIZE VARIABLES FOR CRITICAL VOLUME LOOP */
/*     The iterative process of the VPM Critical Volume Algorithm will operate */
/*     only in the decompression zone since it deals with excess gas volume */
/*     released as a result of supersaturation gradients (not possible below the */
/*     decompression zone). */
/* =============================================================================== */

	gas_loadings_ascent_descen(&starting_depth, &depth_start_of_deco_zone, &rate);
	run_time_start_of_deco_zone = run_time;
	deco_phase_volume_time = 0.;
	last_run_time = 0.;
	schedule_converged = FALSE_;
	for (i = 1; i <= 16; ++i) {
	    last_phase_volume_time[i - 1] = 0.;
	    he_pressure_start_of_deco_zone[i - 1] = 
		    helium_pressure[i - 1];
	    n2_pressure_start_of_deco_zone[i - 1] = 
		    nitrogen_pressure[i - 1];
	    max_actual_gradient[i - 1] = 0.;
	}
/* =============================================================================== */
/*     START OF CRITICAL VOLUME LOOP */
/*     This loop operates between Lines 50 and 100.  If the Critical Volume */
/*     Algorithm is toggled "off" in the program settings, there will only be */
/*     one pass through this loop.  Otherwise, there will be two or more passes */
/*     through this loop until the deco schedule is "converged" - that is when a */
/*     comparison between the phase volume time of the present iteration and the */
/*     last iteration is less than or equal to one minute.  This implies that */
/*     the volume of released gas in the most recent iteration differs from the */
/*     "critical" volume limit by an acceptably small amount.  The critical */
/*     volume limit is set by the Critical Volume Parameter Lambda in the program */
/*     settings (default setting is 7500 fsw-min with adjustability range from */
/*     from 6500 to 8300 fsw-min according to Bruce Wienke). */
/* =============================================================================== */
/* L50: */

	while(TRUE_) {          /* loop will run continuous there is an exit stateme */

/* =============================================================================== */
/*     CALCULATE CURRENT DECO CEILING BASED ON ALLOWABLE SUPERSATURATION */
/*     GRADIENTS AND SET FIRST DECO STOP.  CHECK TO MAKE SURE THAT SELECTED STEP */
/*     SIZE WILL NOT ROUND UP FIRST STOP TO A DEPTH THAT IS BELOW THE DECO ZONE. */
/* =============================================================================== */

	    calc_deco_ceiling(&deco_ceiling_depth);
	    if (deco_ceiling_depth <= 0.) {
		    deco_stop_depth = 0.;
	    } else {
		    rounding_operation2 = 
                deco_ceiling_depth / step_size + ( float).5;
		    deco_stop_depth = 
                r_nint(&rounding_operation2) * step_size;
	    }
	    if (deco_stop_depth > depth_start_of_deco_zone) {
	        printf("\t\n");
            printf(fmt_905);
            printf(fmt_900);    
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
	    }

/* =============================================================================== */
/*     PERFORM A SEPARATE "PROJECTED ASCENT" OUTSIDE OF THE MAIN PROGRAM TO MAKE */
/*     SURE THAT AN INCREASE IN GAS LOADINGS DURING ASCENT TO THE FIRST STOP WILL */
/*     NOT CAUSE A VIOLATION OF THE DECO CEILING.  IF SO, ADJUST THE FIRST STOP */
/*     DEEPER BASED ON STEP SIZE UNTIL A SAFE ASCENT CAN BE MADE. */
/*     Note: this situation is a possibility when ascending from extremely deep */
/*     dives or due to an unusual gas mix selection. */
/*     CHECK AGAIN TO MAKE SURE THAT ADJUSTED FIRST STOP WILL NOT BE BELOW THE */
/*     DECO ZONE. */
/* =============================================================================== */

	    projected_ascent(&depth_start_of_deco_zone, &rate, &deco_stop_depth, &step_size);
	    if (deco_stop_depth > depth_start_of_deco_zone) {
	        printf("\t\n");
            printf(fmt_905);
            printf(fmt_900);    
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
	    }

/* =============================================================================== */
/*     HANDLE THE SPECIAL CASE WHEN NO DECO STOPS ARE REQUIRED - ASCENT CAN BE */
/*     MADE DIRECTLY TO THE SURFACE */
/*     Write ascent data to output file and exit the Critical Volume Loop. */
/* =============================================================================== */

	    if (deco_stop_depth == 0.) {
		    for (i = 1; i <= 16; ++i) {
		        helium_pressure[i - 1] = 
			        he_pressure_start_of_ascent[i - 1];
		        nitrogen_pressure[i - 1] = 
			        n2_pressure_start_of_ascent[i - 1];
		    }
		    run_time = run_time_start_of_ascent;
		    segment_number = 
			    segment_number_start_of_ascent;
		    starting_depth = depth_change[0];
		    ending_depth = 0.;
		    gas_loadings_ascent_descen(&starting_depth, 
                                       &ending_depth, &rate);

		    printf(fmt_860, segment_number, 
                            segment_time, 
                            run_time, 
                            mix_number, 
                            deco_stop_depth,
                            rate);
		    break;          /* exit the critical volume l */
	    }

/* =============================================================================== */
/*     ASSIGN VARIABLES FOR ASCENT FROM START OF DECO ZONE TO FIRST STOP.  SAVE */
/*     FIRST STOP DEPTH FOR LATER USE WHEN COMPUTING THE FINAL ASCENT PROFILE */
/* =============================================================================== */

	    starting_depth = depth_start_of_deco_zone;
	    first_stop_depth = deco_stop_depth;

/* =============================================================================== */
/*     DECO STOP LOOP BLOCK WITHIN CRITICAL VOLUME LOOP */
/*     This loop computes a decompression schedule to the surface during each */
/*     iteration of the critical volume loop.  No output is written from this */
/*     loop, rather it computes a schedule from which the in-water portion of the */
/*     total phase volume time (Deco_Phase_Volume_Time) can be extracted.  Also, */
/*     the gas loadings computed at the end of this loop are used the subroutine */
/*     which computes the out-of-water portion of the total phase volume time */
/*     (Surface_Phase_Volume_Time) for that schedule. */

/*     Note that exit is made from the loop after last ascent is made to a deco */
/*     stop depth that is less than or equal to zero.  A final deco stop less */
/*     than zero can happen when the user makes an odd step size change during */
/*     ascent - such as specifying a 5 msw step size change at the 3 msw stop! */
/* =============================================================================== */

	    while(TRUE_) {          /* loop will run continuous there is an break statement */
		    gas_loadings_ascent_descen(&starting_depth, &deco_stop_depth, &rate);
		    if (deco_stop_depth <= 0.) {
	    	    break;
		    }
		    if (number_of_changes > 1) {
		        i1 = number_of_changes;
		        for (i = 2; i <= i1; ++i) {
			        if (depth_change[i - 1] >= deco_stop_depth) {
			            mix_number = mix_change[i - 1];
			            rate = rate_change[i - 1];
			            step_size = step_size_change[i - 1];
			        }
		        }
		    }
		    decompression_stop(&deco_stop_depth, &step_size);
		    starting_depth = deco_stop_depth;
		    next_stop = deco_stop_depth - step_size;
		    deco_stop_depth = next_stop;
		    last_run_time = run_time;
/* L60: */
	    }

/* =============================================================================== */
/*     COMPUTE TOTAL PHASE VOLUME TIME AND MAKE CRITICAL VOLUME COMPARISON */
/*     The deco phase volume time is computed from the run time.  The surface */
/*     phase volume time is computed in a subroutine based on the surfacing gas */
/*     loadings from previous deco loop block.  Next the total phase volume time */
/*     (in-water + surface) for each compartment is compared against the previous */
/*     total phase volume time.  The schedule is converged when the difference is */
/*     less than or equal to 1 minute in any one of the 16 compartments. */

/*     Note:  the "phase volume time" is somewhat of a mathematical concept. */
/*     It is the time divided out of a total integration of supersaturation */
/*     gradient x time (in-water and surface).  This integration is multiplied */
/*     by the excess bubble number to represent the amount of free-gas released */
/*     as a result of allowing a certain number of excess bubbles to form. */
/* =============================================================================== */
/* end of deco stop loop */

	    deco_phase_volume_time = run_time - run_time_start_of_deco_zone;
	    calc_surface_phase_volume_time();

	    for (i = 1; i <= 16; ++i) {
		    phase_volume_time[i - 1] = 
                deco_phase_volume_time + surface_phase_volume_time[i - 1];
		    critical_volume_comparison = 
                (r1 = phase_volume_time[i - 1] - last_phase_volume_time[i - 1], fabs(r1));
			
		    if (critical_volume_comparison <= 1.) {
		        schedule_converged = TRUE_;
		    }
	    }

/* =============================================================================== */
/*     CRITICAL VOLUME DECISION TREE BETWEEN LINES 70 AND 99 */
/*     There are two options here.  If the Critical Volume Agorithm setting is */
/*     "on" and the schedule is converged, or the Critical Volume Algorithm */
/*     setting was "off" in the first place, the program will re-assign variables */
/*     to their values at the start of ascent (end of bottom time) and process */
/*     a complete decompression schedule once again using all the same ascent */
/*     parameters and first stop depth.  This decompression schedule will match */
/*     the last iteration of the Critical Volume Loop and the program will write */
/*     the final deco schedule to the output file. */

/*     Note: if the Critical Volume Agorithm setting was "off", the final deco */
/*     schedule will be based on "Initial Allowable Supersaturation Gradients." */
/*     If it was "on", the final schedule will be based on "Adjusted Allowable */
/*     Supersaturation Gradients" (gradients that are "relaxed" as a result of */
/*     the Critical Volume Algorithm). */

/*     If the Critical Volume Agorithm setting is "on" and the schedule is not */
/*     converged, the program will re-assign variables to their values at the */
/*     start of the deco zone and process another trial decompression schedule. */
/* =============================================================================== */
/* L70: */

	    if (schedule_converged || critical_volume_algorithm_off) {
		    for (i = 1; i <= 16; ++i) {
		        helium_pressure[i - 1] = 
			        he_pressure_start_of_ascent[i - 1];
		        nitrogen_pressure[i - 1] = 
			        n2_pressure_start_of_ascent[i - 1];
		    }
		    run_time = run_time_start_of_ascent;
		    segment_number = 
			    segment_number_start_of_ascent;
		    starting_depth = depth_change[0];
		    mix_number = mix_change[0];
		    rate = rate_change[0];
		    step_size = step_size_change[0];
		    deco_stop_depth = first_stop_depth;
		    last_run_time = 0.;

/* =============================================================================== */
/*     DECO STOP LOOP BLOCK FOR FINAL DECOMPRESSION SCHEDULE */
/* =============================================================================== */

		while(TRUE_) {      /* loop will run continuous until there is an break statement */
		    gas_loadings_ascent_descen(&starting_depth, 
                                       &deco_stop_depth, &rate);

/* =============================================================================== */
/*     DURING FINAL DECOMPRESSION SCHEDULE PROCESS, COMPUTE MAXIMUM ACTUAL */
/*     SUPERSATURATION GRADIENT RESULTING IN EACH COMPARTMENT */
/*     If there is a repetitive dive, this will be used later in the VPM */
/*     Repetitive Algorithm to adjust the values for critical radii. */
/* =============================================================================== */

		    calc_max_actual_gradient(&deco_stop_depth);

		    printf(fmt_860, segment_number, 
                            segment_time, 
                            run_time, 
                            mix_number, 
                            deco_stop_depth, 
                            rate);
		    if (deco_stop_depth <= 0.) {
		        break;
		    }
		    if (number_of_changes > 1) {
			    i1 = number_of_changes;
			    for (i = 2; i <= i1; ++i) {
			        if (depth_change[i - 1] >= deco_stop_depth) {
			        	mix_number = mix_change[i - 1];
				        rate = rate_change[i - 1];
				        step_size = step_size_change[i - 1];
			        }
			    }
		    }
		    decompression_stop(&deco_stop_depth, 
                               &step_size);

/* =============================================================================== */
/*     This next bit justs rounds up the stop time at the first stop to be in */
/*     whole increments of the minimum stop time (to make for a nice deco table). */
/* =============================================================================== */

		    if (last_run_time == 0.) {
			    r1 = segment_time / minimum_deco_stop_time + ..5;
			    stop_time = r_nint(&r1) * minimum_deco_stop_time;
		    } else {
			    stop_time = run_time - last_run_time;
		    }

/* =============================================================================== */
/*     DURING FINAL DECOMPRESSION SCHEDULE, IF MINIMUM STOP TIME PARAMETER IS A */
/*     WHOLE NUMBER (i.e. 1 minute) THEN WRITE DECO SCHEDULE USING INTEGER */
/*     NUMBERS (looks nicer).  OTHERWISE, USE DECIMAL NUMBERS. */
/*     Note: per the request of a noted exploration diver(!), program now allows */
/*     a minimum stop time of less than one minute so that total ascent time can */
/*     be minimized on very long dives.  In fact, with step size set at 1 fsw or */
/*     0.2 msw and minimum stop time set at 0.1 minute (6 seconds), a near */
/*     continuous decompression schedule can be computed. */
/* =============================================================================== */

		    if (r_int(&minimum_deco_stop_time) == minimum_deco_stop_time) {
			    printf(fmt_862,    segment_number, 
                                   segment_time, 
                                   run_time, 
                                   mix_number,  
                                   (int)deco_stop_depth,
			                       (int)stop_time,
			                       (int)run_time);
		    } else {
			    printf(fmt_863,    segment_number, 
                                   segment_time, 
                                   run_time, 
                                   mix_number, 
                                   deco_stop_depth, 
                                   stop_time, 
                                   run_time);
		    }
		    starting_depth = deco_stop_depth;
		    next_stop = deco_stop_depth - step_size;
		    deco_stop_depth = next_stop;
		    last_run_time = run_time;
/* L80: */
		}
/* for final deco sche */
/* end of deco stop lo */
		break;
/* final deco schedule */
/* exit critical volume l */

/* =============================================================================== */
/*     IF SCHEDULE NOT CONVERGED, COMPUTE RELAXED ALLOWABLE SUPERSATURATION */
/*     GRADIENTS WITH VPM CRITICAL VOLUME ALGORITHM AND PROCESS ANOTHER */
/*     ITERATION OF THE CRITICAL VOLUME LOOP */
/* =============================================================================== */

	    } else {
		    critical_volume(&deco_phase_volume_time);

		    deco_phase_volume_time = 0.;
		    run_time = run_time_start_of_deco_zone;
		    starting_depth = depth_start_of_deco_zone;
		    mix_number = mix_change[0];
		    rate = rate_change[0];
		    step_size = step_size_change[0];
		    for (i = 1; i <= 16; ++i) {
		        last_phase_volume_time[i - 1] = phase_volume_time[i - 1];
		        helium_pressure[i - 1] = 
			        he_pressure_start_of_deco_zone[i - 1];
		        nitrogen_pressure[i - 1] = 
			        n2_pressure_start_of_deco_zone[i - 1];
		    }
		    continue;
/* (Line 50) to process another it */
/* Return to start of criti */
/* L99: */
	    }
/* end of critical volume decision */
/* L100: */
    }                                       /* end of critical vol loop */

/* =============================================================================== */
/*     PROCESSING OF DIVE COMPLETE.  READ INPUT FILE TO DETERMINE IF THERE IS A */
/*     REPETITIVE DIVE.  IF NONE, THEN EXIT REPETITIVE LOOP. */
/* =============================================================================== */

    //  READ (7,*) Repetitive_Dive_Flag
    if (err = ReadUnformatted(fd, &line9_ul)) {
        printf("\nError %d returned reading VPMDECO.IN", err);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
    }       
	if (repetitive_dive_flag == 0) {
		break;                              /* at Line 330 exit repet */

/* =============================================================================== */
/*     IF THERE IS A REPETITIVE DIVE, COMPUTE GAS LOADINGS (OFF-GASSING) DURING */
/*     SURFACE INTERVAL TIME.  ADJUST CRITICAL RADII USING VPM REPETITIVE */
/*     ALGORITHM.  RE-INITIALIZE SELECTED VARIABLES AND RETURN TO START OF */
/*     REPETITIVE LOOP AT LINE 30. */
/* =============================================================================== */

	} else if (repetitive_dive_flag == 1) {
        // READ (7,*) Surface_Interval_Time
        if (err = ReadUnformatted(fd, &line10_ul)) {
            printf("\nError %d returned reading VPMDECO.IN", err);
            printf("\nPROGRAM TERMINATED\n");
            exit(1);
        }       

	    gas_loadings_surface_interval(&surface_interval_time);

	    vpm_repetitive_algorithm(&surface_interval_time);

	    for (i = 1; i <= 16; ++i) {
		    max_crushing_pressure_he[i - 1] = 0.;
		    max_crushing_pressure_n2[i - 1] = 0.;
		    max_actual_gradient[i - 1] = 0.;
	    }
	    run_time = 0.;
	    segment_number = 0;
        printf(fmt_880);
        printf(fmt_890);
        printf(fmt_813);
        continue;

/* =============================================================================== */
/*     WRITE ERROR MESSAGE AND TERMINATE PROGRAM IF THERE IS AN ERROR IN THE */
/*     INPUT FILE FOR THE REPETITIVE DIVE FLAG */
/* =============================================================================== */

/* Return to start of repetitive loop to proce */
	} else {
	    printf("\t\n");
        printf(fmt_908);
        printf(fmt_900);
        printf("\nPROGRAM TERMINATED\n");
        exit(1);
	}
/* L330: */
    }

/* =============================================================================== */
/*     FINAL WRITES TO OUTPUT AND CLOSE PROGRAM FILES */
/* =============================================================================== */
/* End of repetit */

    printf(fmt_813);
    printf(fmt_871);
    printf(fmt_872);
    printf(fmt_813);
    printf(fmt_880);
    return 0;
} /* MAIN__ */