sub _equation_air {
        my ( $args_ref, $level_iii_ref, $y ) = @_;
#    my $_equation_air = sub {
#        my $y = shift;
#    local *_equation_air = sub {
#        my $y = shift;
...
[download]

package Levels::Level_IV;

#--------------------------------------------------------------

#Heading:       Level_VI.pl

#Author:        Anonymous

#Purpose:       Level_VI calculations.

#Usage:         It's a library. Use only level_iv subroutine.

#References:    Don Mackay. Multimedia Environmental Models
#               The Fugacity Approach Second Edition

#File formats:  No files needed or printed.

#Restrictions:  Reasonal input.

#Notes:         This module is free software; you can redistribute it
#               and/or modify it under the same terms as Perl itself.

#--------------------------------------------------------------

use strict;
use warnings;
#use diagnostics;
use Carp;
use Exporter;
#use List::Util;
use List::MoreUtils qw(zip);
use Scalar::Util qw(reftype);
#use Math::ODE; #Solve N-th Order (O)rdinary (D)ifferential (E)quation
+s
#use Math::Calculus::Differentiate; #Algebraic Differentiation Engine
use Math::RungeKutta; #Runge-Kutta numeric algo. for Differential Equa
+tions
use Term::Clui;
use FileHandle;
use Smart::Comments;
use Contextual::Return;
use Levels::Level_III;

our ( @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS, $VERSION );

$VERSION = 0.9;
@ISA = qw(Exporter);

@EXPORT = qw(level_iv);

#non equilibrium unsteady-state flow system
sub level_iv {
    my ( $args_ref, $R, $factor, $level_iii_ref ) = @_;
    no warnings qw(closure); #'variable will not stay shared' -> switc
+hed off
    #see end of file for possible solution

    #call level_iii subroutine to generate difusivities and fugacities
    $level_iii_ref = level_iii($args_ref, $R, $factor) unless $level_i
+ii_ref;

    #choose differentiation algorithm
    print "\n";
    my $routine 
        = choose( 'Algorithm: (use arrow keys and <Return>, or q to qu
+it)',
        'rk4', 'rk4_auto' );

    #initialize y's for the different mediums
    my @y_med;
    INPUTY:
    while (1) {
        my @try = split m/\s*,\s*/, 
            ask( "Initialize ya, yw, ys and yse: (separate y's with ',
+')" );
        if ( scalar @try < 4 ) {
            print "Four (4) arguments needed. Your input was: '@try'.\
+n";
            redo INPUTY;
        }
        elsif ( "$try[0]$try[1]$try[2]$try[3]" !~ m/ \A \d+ \z /xms ) 
+{
            print "Arguments are not numeric. Your input was: '@try'.\
+n";
            redo INPUTY;            
        }
        else {
            @y_med = @try;
            last INPUTY;
        }
    }

    #choose time step, starting and end time
    my $dt = choose( 'Time step: (use arrow keys and <Return>, or q to
+ quit)',
        '0.005', '0.01', '0.02', '0.05', '0.1', '0.2' );
    my $t = 0;
    my $t_final = 1 / $dt;  # from 0+step to 1 in that many steps

    #open file for printing runge-kutta-merson results
    open my $fh, '>', 'Level_IVa_Caculations.txt' or croak $!;
    $fh->autoflush(1);

    #Runge-Kutta-Merson 4th-order method
    my @iterations;
    my @output_merson;

    if ( $routine eq 'rk4_auto' ) {
        #choose error for the steps in auto mode.
        my $eps 
            = choose( 'epsilon: (use arrow keys and <Return>, or q to 
+quit)', 
            '0.000001', '0.00001', '0.0001' );

        #invoke runge kutta algorithm and store calculations
        RK4_AUTO:
        for my $i ( $t..$t_final ) {    ### Differentiating:...   % do
+ne
            #print initialisation 
            if ( $i == 0 ) {
                printf {$fh} "i = %3d  t = %16s  dt = %16s  ya = %20s 
+ "
                    . "yw = %20s  ys = %20s  yse = %20s\n", $i, $t, $d
+t,  
                    $y_med[0], $y_med[1], $y_med[2], $y_med[3];
                push @iterations, $i;
                push @output_merson, 
                    [ $t, $dt, $y_med[0], $y_med[1], $y_med[2], $y_med
+[3] ];
                next RK4_AUTO;
            }

            #begin calculations after initial print
            ( $t, $dt, @y_med ) = rk4_auto( \@y_med, \&_dydt, $t, $dt,
+ $eps );
            printf {$fh} "i = %3d  t = %16s  dt = %16s  ya = %20s  yw 
+= %20s"
                    . "  ys = %20s  yse = %20s\n", $i, $t, $dt, $y_med
+[0], 
                    $y_med[1], $y_med[2], $y_med[3];
            push @iterations, $i;
            push @output_merson, 
                [ $t, $dt, $y_med[0], $y_med[1], $y_med[2], $y_med[3] 
+];
        }
    }
    else {
        #invoke runge kutta algorithm and store calculations
        RK4:
        for my $i ( $t..$t_final ) {    ### Differentiating:...   % do
+ne
            #print initialisation 
            if ( $i == 0 ) {
                printf {$fh} "i = %3d  t = %.3f  dt = %.3f  ya = %20s 
+ "
                    . "yw = %20s  ys = %20s  yse = %20s\n", $i, $t, $d
+t,  
                    $y_med[0], $y_med[1], $y_med[2], $y_med[3];
                push @iterations, $i;
                push @output_merson, 
                    [ $t, $dt, $y_med[0], $y_med[1], $y_med[2], $y_med
+[3] ];
                next RK4;
            }

            #begin calculations after initial print
            ( $t, @y_med ) = rk4( \@y_med, \&_dydt, $t, $dt );
            printf {$fh} "i = %3d  t = %.3f  dt = %.3f  ya = %20s  yw 
+= %20s"
                    . "  ys = %20s  yse = %20s\n", $i, $t, $dt, $y_med
+[0], 
                    $y_med[1], $y_med[2], $y_med[3];
            push @iterations, $i;
            push @output_merson, 
                [ $t, $dt, $y_med[0], $y_med[1], $y_med[2], $y_med[3] 
+];
        }
    }

    print "File Level_IVa_Caculations.txt written.\n";
    close $fh or croak $!;
    
    #_equation and _dydt_* have to be declared in the main routine or 
+they 
    #would not have access to the $args_ref and $level_iii_ref (closur
+es)

    #level_iv '_equation' for mass is used as helper equation in dydt 
+function 
    sub _equation_air {
        my ( $args_ref, $level_iii_ref, $y ) = @_;
#    my $_equation_air = sub {
#        my $y = shift;
#    local *_equation_air = sub {
#        my $y = shift;

        #for easy calculations make copies of the parameters (not alia
+ses)
        my ( $Ei, $Gi, $Cbi, $Vi, $Zi );
        if (    exists $args_ref->{'-Ea'}
            and exists $args_ref->{'-Ga'}
            and exists $args_ref->{'-Cba'}
            and exists $args_ref->{'-Va'}
            and exists $args_ref->{'-Za'} ) 
        {
            $Ei  = $args_ref->{'-Ea'};
            $Gi  = $args_ref->{'-Ga'};
            $Cbi = $args_ref->{'-Cba'};
            $Vi  = $args_ref->{'-Va'};
            $Zi  = $args_ref->{'-Za'};
        }
        else {
            croak "No parameter '-Ea', '-Ga', '-Cba', '-Va' or '-Za' f
+ound. "
                . "Differentiation not possible";
        }
        my $f_j1    = $level_iii_ref->[0]{'air'}{'fugacity'};
        my $D_j1_i1 = 0;
        my $D_j1_i2 = $level_iii_ref->[1]{'diffusivity'}{'air-water'};
        my $D_j1_i3 = $level_iii_ref->[1]{'diffusivity'}{'air-soil'};
        my $D_j1_i4 = 0;
        my $f_j2    = $level_iii_ref->[0]{'water'}{'fugacity'};
        my $D_j2_i1 = $level_iii_ref->[1]{'diffusivity'}{'water-air'};
        my $D_j2_i2 = 0;
        my $D_j2_i3 = $level_iii_ref->[1]{'diffusivity'}{'water-soil'}
+;
        my $D_j2_i4 = $level_iii_ref->[1]{'diffusivity'}{'water-sedime
+nt'};
        my $f_j3    = $level_iii_ref->[0]{'soil'}{'fugacity'};
        my $D_j3_i1 = $level_iii_ref->[1]{'diffusivity'}{'soil-air'};
        my $D_j3_i2 = $level_iii_ref->[1]{'diffusivity'}{'soil-water'}
+;
        my $D_j3_i3 = 0;
        my $D_j3_i4 = 0;
        my $f_j4    = $level_iii_ref->[0]{'sediment'}{'fugacity'};
        my $D_j4_i1 = 0;
        my $D_j4_i2 = $level_iii_ref->[1]{'diffusivity'}{'sediment-wat
+er'};
        my $D_j4_i3 = 0;
        my $D_j4_i4 = 0;
        my $Dti     = $level_iii_ref->[0]{'air'}{'D_total'};
    
        #differential equation for mass
        my $dydt 
            = ( ( $Ei + $Gi * $Cbi ) 
                + ( $D_j1_i1 * $f_j1 ) 
                + ( $D_j1_i2 * $f_j1 ) 
                + ( $D_j1_i3 * $f_j1 ) 
                + ( $D_j1_i4 * $f_j1 ) 
                + ( $D_j2_i1 * $f_j2 )  
                + ( $D_j2_i2 * $f_j2 ) 
                + ( $D_j2_i3 * $f_j2 ) 
                + ( $D_j2_i4 * $f_j2 ) 
                + ( $D_j3_i1 * $f_j3 ) 
                + ( $D_j3_i2 * $f_j3 ) 
                + ( $D_j3_i3 * $f_j3 ) 
                + ( $D_j3_i4 * $f_j3 ) 
                + ( $D_j4_i1 * $f_j4 ) 
                + ( $D_j4_i2 * $f_j4 ) 
                + ( $D_j4_i3 * $f_j4 ) 
                + ( $D_j4_i4 * $f_j4 )
                - ( $Dti * $y ) ) / ( $Vi * $Zi );

        return $dydt;
    };

    #helper _equation for water medium/compartment
    sub _equation_water {
        my ( $args_ref, $level_iii_ref, $y ) = @_;

        #for easy calculations make copies of the parameters (not alia
+ses)
        my ( $Ei, $Gi, $Cbi, $Vi, $Zi );
        if (    exists $args_ref->{'-Ew'}
            and exists $args_ref->{'-Gw'}
            and exists $args_ref->{'-Cbw'}
            and exists $args_ref->{'-Vw'}
            and exists $args_ref->{'-Zw'} ) 
        {
            $Ei  = $args_ref->{'-Ew'};
            $Gi  = $args_ref->{'-Gw'};
            $Cbi = $args_ref->{'-Cbw'};
            $Vi  = $args_ref->{'-Vw'};
            $Zi  = $args_ref->{'-Zw'};
        }
        else {
            croak "No parameter '-Ew', '-Gw', '-Cbw', '-Vw' or '-Zw' f
+ound. "
                . "Differentiation not possible";
        }
        my $f_j1    = $level_iii_ref->[0]{'air'}{'fugacity'};
        my $D_j1_i1 = 0;
        my $D_j1_i2 = $level_iii_ref->[1]{'diffusivity'}{'air-water'};
        my $D_j1_i3 = $level_iii_ref->[1]{'diffusivity'}{'air-soil'};
        my $D_j1_i4 = 0;
        my $f_j2    = $level_iii_ref->[0]{'water'}{'fugacity'};
        my $D_j2_i1 = $level_iii_ref->[1]{'diffusivity'}{'water-air'};
        my $D_j2_i2 = 0;
        my $D_j2_i3 = $level_iii_ref->[1]{'diffusivity'}{'water-soil'}
+;
        my $D_j2_i4 = $level_iii_ref->[1]{'diffusivity'}{'water-sedime
+nt'};
        my $f_j3    = $level_iii_ref->[0]{'soil'}{'fugacity'};
        my $D_j3_i1 = $level_iii_ref->[1]{'diffusivity'}{'soil-air'};
        my $D_j3_i2 = $level_iii_ref->[1]{'diffusivity'}{'soil-water'}
+;
        my $D_j3_i3 = 0;
        my $D_j3_i4 = 0;
        my $f_j4    = $level_iii_ref->[0]{'sediment'}{'fugacity'};
        my $D_j4_i1 = 0;
        my $D_j4_i2 = $level_iii_ref->[1]{'diffusivity'}{'sediment-wat
+er'};
        my $D_j4_i3 = 0;
        my $D_j4_i4 = 0;
        my $Dti     = $level_iii_ref->[0]{'water'}{'D_total'};
    
        #differential equation for mass
        my $dydt 
            = ( ( $Ei + $Gi * $Cbi ) 
                + ( $D_j1_i1 * $f_j1 ) 
                + ( $D_j1_i2 * $f_j1 ) 
                + ( $D_j1_i3 * $f_j1 ) 
                + ( $D_j1_i4 * $f_j1 ) 
                + ( $D_j2_i1 * $f_j2 )  
                + ( $D_j2_i2 * $f_j2 ) 
                + ( $D_j2_i3 * $f_j2 ) 
                + ( $D_j2_i4 * $f_j2 ) 
                + ( $D_j3_i1 * $f_j3 ) 
                + ( $D_j3_i2 * $f_j3 ) 
                + ( $D_j3_i3 * $f_j3 ) 
                + ( $D_j3_i4 * $f_j3 ) 
                + ( $D_j4_i1 * $f_j4 ) 
                + ( $D_j4_i2 * $f_j4 ) 
                + ( $D_j4_i3 * $f_j4 ) 
                + ( $D_j4_i4 * $f_j4 )
                - ( $Dti * $y ) ) / ( $Vi * $Zi );

        return $dydt;
    }

    #helper _equation for soil medium/compartment
    sub _equation_soil {
        my ( $args_ref, $level_iii_ref, $y ) = @_;

        #for easy calculations make copies of the parameters (not alia
+ses)
        my ( $Ei, $Gi, $Cbi, $Vi, $Zi );
        if (    exists $args_ref->{'-Es'}
            and exists $args_ref->{'-Gs'}
            and exists $args_ref->{'-Cbs'}
            and exists $args_ref->{'-Vs'}
            and exists $args_ref->{'-Zs'} ) 
        {
            $Ei  = $args_ref->{'-Es'};
            $Gi  = $args_ref->{'-Gs'};
            $Cbi = $args_ref->{'-Cbs'};
            $Vi  = $args_ref->{'-Vs'};
            $Zi  = $args_ref->{'-Zs'};
        }
        else {
            croak "No parameter '-Es', '-Gs', '-Cbs', '-Vs' or '-Zs' f
+ound. "
                . "Differentiation not possible";
        }
        my $f_j1    = $level_iii_ref->[0]{'air'}{'fugacity'};
        my $D_j1_i1 = 0;
        my $D_j1_i2 = $level_iii_ref->[1]{'diffusivity'}{'air-water'};
        my $D_j1_i3 = $level_iii_ref->[1]{'diffusivity'}{'air-soil'};
        my $D_j1_i4 = 0;
        my $f_j2    = $level_iii_ref->[0]{'water'}{'fugacity'};
        my $D_j2_i1 = $level_iii_ref->[1]{'diffusivity'}{'water-air'};
        my $D_j2_i2 = 0;
        my $D_j2_i3 = $level_iii_ref->[1]{'diffusivity'}{'water-soil'}
+;
        my $D_j2_i4 = $level_iii_ref->[1]{'diffusivity'}{'water-sedime
+nt'};
        my $f_j3    = $level_iii_ref->[0]{'soil'}{'fugacity'};
        my $D_j3_i1 = $level_iii_ref->[1]{'diffusivity'}{'soil-air'};
        my $D_j3_i2 = $level_iii_ref->[1]{'diffusivity'}{'soil-water'}
+;
        my $D_j3_i3 = 0;
        my $D_j3_i4 = 0;
        my $f_j4    = $level_iii_ref->[0]{'sediment'}{'fugacity'};
        my $D_j4_i1 = 0;
        my $D_j4_i2 = $level_iii_ref->[1]{'diffusivity'}{'sediment-wat
+er'};
        my $D_j4_i3 = 0;
        my $D_j4_i4 = 0;
        my $Dti     = $level_iii_ref->[0]{'soil'}{'D_total'};
    
        #differential equation for mass
        my $dydt 
            = ( ( $Ei + $Gi * $Cbi ) 
                + ( $D_j1_i1 * $f_j1 ) 
                + ( $D_j1_i2 * $f_j1 ) 
                + ( $D_j1_i3 * $f_j1 ) 
                + ( $D_j1_i4 * $f_j1 ) 
                + ( $D_j2_i1 * $f_j2 )  
                + ( $D_j2_i2 * $f_j2 ) 
                + ( $D_j2_i3 * $f_j2 ) 
                + ( $D_j2_i4 * $f_j2 ) 
                + ( $D_j3_i1 * $f_j3 ) 
                + ( $D_j3_i2 * $f_j3 ) 
                + ( $D_j3_i3 * $f_j3 ) 
                + ( $D_j3_i4 * $f_j3 ) 
                + ( $D_j4_i1 * $f_j4 ) 
                + ( $D_j4_i2 * $f_j4 ) 
                + ( $D_j4_i3 * $f_j4 ) 
                + ( $D_j4_i4 * $f_j4 )
                - ( $Dti * $y ) ) / ( $Vi * $Zi );

        return $dydt;
    }

    #helper _equation for sediment medium/compartment
    sub _equation_sediment {
        my ( $args_ref, $level_iii_ref, $y ) = @_;

        #for easy calculations make copies of the parameters (not alia
+ses)
        my ( $Ei, $Gi, $Cbi, $Vi, $Zi );
        if (    exists $args_ref->{'-Ese'}
            and exists $args_ref->{'-Gse'}
            and exists $args_ref->{'-Cbse'}
            and exists $args_ref->{'-Vse'}
            and exists $args_ref->{'-Zse'} ) 
        {
            $Ei  = $args_ref->{'-Ese'};
            $Gi  = $args_ref->{'-Gse'};
            $Cbi = $args_ref->{'-Cbse'};
            $Vi  = $args_ref->{'-Vse'};
            $Zi  = $args_ref->{'-Zse'};
        }
        else {
            croak "No parameter '-Ese', '-Ges', '-Cbse', '-Vse' or '-Z
+se' "
                . "found. Differentiation not possible";
        }
        my $f_j1    = $level_iii_ref->[0]{'air'}{'fugacity'};
        my $D_j1_i1 = 0;
        my $D_j1_i2 = $level_iii_ref->[1]{'diffusivity'}{'air-water'};
        my $D_j1_i3 = $level_iii_ref->[1]{'diffusivity'}{'air-soil'};
        my $D_j1_i4 = 0;
        my $f_j2    = $level_iii_ref->[0]{'water'}{'fugacity'};
        my $D_j2_i1 = $level_iii_ref->[1]{'diffusivity'}{'water-air'};
        my $D_j2_i2 = 0;
        my $D_j2_i3 = $level_iii_ref->[1]{'diffusivity'}{'water-soil'}
+;
        my $D_j2_i4 = $level_iii_ref->[1]{'diffusivity'}{'water-sedime
+nt'};
        my $f_j3    = $level_iii_ref->[0]{'soil'}{'fugacity'};
        my $D_j3_i1 = $level_iii_ref->[1]{'diffusivity'}{'soil-air'};
        my $D_j3_i2 = $level_iii_ref->[1]{'diffusivity'}{'soil-water'}
+;
        my $D_j3_i3 = 0;
        my $D_j3_i4 = 0;
        my $f_j4    = $level_iii_ref->[0]{'sediment'}{'fugacity'};
        my $D_j4_i1 = 0;
        my $D_j4_i2 = $level_iii_ref->[1]{'diffusivity'}{'sediment-wat
+er'};
        my $D_j4_i3 = 0;
        my $D_j4_i4 = 0;
        my $Dti     = $level_iii_ref->[0]{'sediment'}{'D_total'};
    
        #differential equation for mass
        my $dydt 
            = ( ( $Ei + $Gi * $Cbi ) 
                + ( $D_j1_i1 * $f_j1 ) 
                + ( $D_j1_i2 * $f_j1 ) 
                + ( $D_j1_i3 * $f_j1 ) 
                + ( $D_j1_i4 * $f_j1 ) 
                + ( $D_j2_i1 * $f_j2 )  
                + ( $D_j2_i2 * $f_j2 ) 
                + ( $D_j2_i3 * $f_j2 ) 
                + ( $D_j2_i4 * $f_j2 ) 
                + ( $D_j3_i1 * $f_j3 ) 
                + ( $D_j3_i2 * $f_j3 ) 
                + ( $D_j3_i3 * $f_j3 ) 
                + ( $D_j3_i4 * $f_j3 ) 
                + ( $D_j4_i1 * $f_j4 ) 
                + ( $D_j4_i2 * $f_j4 ) 
                + ( $D_j4_i3 * $f_j4 ) 
                + ( $D_j4_i4 * $f_j4 )
                - ( $Dti * $y ) ) / ( $Vi * $Zi );

        return $dydt;
    }

    #differential equation for fugacity at start and end times
    sub _dydt {
        my ( $t, @y ) = @_;
        my @dydt; #e.g.: dy/dt = 4  ->  y'(t) = 4  ->  y(t) = 4t

        #two different stating points are accepted in this implementat
+ion 
        #so two different starting scenarios can be calculated in orde
+r to
        #show how the system runs against a boundary and balances itse
+lf
        
        #start differentiation at point y0 (air)
        $dydt[0] = _equation_air( $args_ref, $level_iii_ref, $y[0] );
        #$dydt[0] = _equation_air($y[0]);   #local
        #$dydt[0] = $_equation_air->($y[0]);#lexical

        #start again at a different point y1 (water)
        $dydt[1] = _equation_water( $args_ref, $level_iii_ref, $y[1] )
+;

        #start again at a different point y2 (soil)
        $dydt[2] = _equation_soil( $args_ref, $level_iii_ref, $y[2] );

        #start again at a different point y3 (sediment)
        $dydt[3] = _equation_sediment( $args_ref, $level_iii_ref, $y[3
+] );

        return @dydt;
    }

    #build resulthahes
    my @hash_pairs = zip @iterations, @output_merson ;    
    my %res_hash = ( 
        air      => { @hash_pairs }, 
        water    => { @hash_pairs },  
        soil     => { @hash_pairs },  
        sediment => { @hash_pairs },  
    );

    #my %res_hash_parameter = ( );

#    return wantarray 
#        ? ( \%res_hash, $level_iii_ref )
#        : [ \%res_hash, $level_iii_ref ];

    return
        SCALAR { [ \%res_hash, $level_iii_ref ] }
        BOOL   { scalar keys %res_hash > 0 }
        NUM    { keys %res_hash }
        STR    { reftype [] }

        LIST   { \%res_hash, $level_iii_ref }

        ARRAYREF  { [ \%res_hash, $level_iii_ref ] }

        CODEREF   { croak " > Don't use this result as code! < "; }
    ;

}

1;

#not correct
#sub aus {
#    my $x = $_[0] + 35;
#    sub in { return $x * 19 }  #beachte ; bei local
#    return $x + in();
#}

#correct
#sub aussen {
#    my $x = $_[0] + 35;
#    local *innen = sub { return $x * 19 };
#    return $x + innen();
#}
[download]