    
    input/output:

         sens, evap, tau_u, tau_v -- real, dimension(:,:)
             Surface fluxes (positive upwards) 
             shape must conform to first two dimensions of input fields

             Units:  sens   : W/(m**2) = J /(m**2 s)
                     evap   : Kg/(m**2 s)
                     tau_u  : Kg/(m s**2) = (Kg m/s)/(m**2 s) 
                               = (Kg/(m**3))((m/s)**2) = Pascals
                     tau_v  : as above

             Flux should be evaluated with atm