[usol,idmn,zn,zm,pertrb]=ortho4(usol,idmn,zn,zm,pertrb);

function [usol,idmn,zn,zm,pertrb]=ortho4(usol,idmn,zn,zm,pertrb);
persistent ete i ifnl ii istr j jfnl jj jstr ute ; 

global spl4_5; if isempty(spl4_5), spl4_5=0; end;
global spl4_6; if isempty(spl4_6), spl4_6=0; end;
global spl4_7; if isempty(spl4_7), spl4_7=0; end;
global spl4_8; if isempty(spl4_8), spl4_8=0; end;
global spl4_15; if isempty(spl4_15), spl4_15=0; end;
global spl4_19; if isempty(spl4_19), spl4_19=0; end;
global spl4_16; if isempty(spl4_16), spl4_16=0; end;
global spl4_20; if isempty(spl4_20), spl4_20=0; end;
if isempty(ete), ete=0; end;
global spl4_17; if isempty(spl4_17), spl4_17=0; end;
global spl4_18; if isempty(spl4_18), spl4_18=0; end;
if isempty(ute), ute=0; end;
if isempty(i), i=0; end;
if isempty(ifnl), ifnl=0; end;
if isempty(ii), ii=0; end;
global spl4_11; if isempty(spl4_11), spl4_11=0; end;
if isempty(istr), istr=0; end;
if isempty(j), j=0; end;
if isempty(jfnl), jfnl=0; end;
if isempty(jj), jj=0; end;
global spl4_13; if isempty(spl4_13), spl4_13=0; end;
if isempty(jstr), jstr=0; end;
global spl4_3; if isempty(spl4_3), spl4_3=0; end;
global spl4_1; if isempty(spl4_1), spl4_1=0; end;
global spl4_2; if isempty(spl4_2), spl4_2=0; end;
global spl4_4; if isempty(spl4_4), spl4_4=0; end;
global spl4_9; if isempty(spl4_9), spl4_9=0; end;
global spl4_12; if isempty(spl4_12), spl4_12=0; end;
global spl4_10; if isempty(spl4_10), spl4_10=0; end;
global spl4_14; if isempty(spl4_14), spl4_14=0; end;
%***BEGIN PROLOGUE  ORTHO4
%***SUBSIDIARY
%***PURPOSE  Subsidiary to SEPX4
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (ORTHO4-S)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     This subroutine orthogonalizes the array USOL with respect to
%     the constant array in a weighted least squares norm.
%
%***SEE ALSO  SEPX4
%***ROUTINES CALLED  (NONE)
%***COMMON BLOCKS    SPL4
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  ORTHO4
%
% common :: ;
%% common /spl4  / kswx , kswy , k , l , ait , bit , cit , dit ,mit , nit , is , ms , js , ns , dlx , dly ,tdlx3 , tdly3 , dlx4 , dly4;
%% common /spl4  / spl4_1 , spl4_2 , spl4_3 , spl4_4 , spl4_5 , spl4_6 , spl4_7 , spl4_8 ,spl4_9 , spl4_10 , spl4_11 , spl4_12 , spl4_13 , spl4_14 , spl4_15 , spl4_16 ,spl4_17 , spl4_18 , spl4_19 , spl4_20;
usol_shape=size(usol);usol=reshape([usol(:).',zeros(1,ceil(numel(usol)./prod([idmn])).*prod([idmn])-numel(usol))],idmn,[]);
zn_shape=size(zn);zn=reshape(zn,1,[]);
zm_shape=size(zm);zm=reshape(zm,1,[]);
%***FIRST EXECUTABLE STATEMENT  ORTHO4
istr = fix(spl4_11);
ifnl = fix(spl4_12);
jstr = fix(spl4_13);
jfnl = fix(spl4_14);
%
%     COMPUTE WEIGHTED INNER PRODUCTS
%
ute = 0.0;
ete = 0.0;
for i = spl4_11 : spl4_12;
ii = fix(i - spl4_11 + 1);
for j = spl4_13 : spl4_14;
jj = fix(j - spl4_13 + 1);
ete = ete + zm(ii).*zn(jj);
ute = ute + usol(i,j).*zm(ii).*zn(jj);
end; j = fix(spl4_14+1);
end; i = fix(spl4_12+1);
%
%     SET PERTURBATION PARAMETER
%
pertrb = ute./ete;
%
%     SUBTRACT OFF CONSTANT PERTRB
%
for i = istr : ifnl;
for j = jstr : jfnl;
usol(i,j) = usol(i,j) - pertrb;
end; j = fix(jfnl+1);
end; i = fix(ifnl+1);
usol_shape=zeros(usol_shape);usol_shape(:)=usol(1:numel(usol_shape));usol=usol_shape;
zn_shape=zeros(zn_shape);zn_shape(:)=zn(1:numel(zn_shape));zn=zn_shape;
zm_shape=zeros(zm_shape);zm_shape(:)=zm(1:numel(zm_shape));zm=zm_shape;
end
%DECK ORTHOG