C
C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/SOSE/BoxAdj/code_ad/Attic/ctrl_init_obcs_variables.F,v 1.8 2011/04/19 23:33:00 mmazloff dead $
C $Name:  $

#include "CTRL_CPPOPTIONS.h"

      subroutine ctrl_init_obcs_variables( mythid )

c     ==================================================================
c     SUBROUTINE ctrl_init_obcs_variables
c     ==================================================================
c
c     o Set parts of the vector of control variables and initialize the
c       rest to zero.
c
c     started: heimbach@mit.edu 25-Mar-2002
c
c     ==================================================================
c     SUBROUTINE ctrl_init_obcs_variables
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"

#include "ctrl.h"
#ifdef ALLOW_OBCS_CONTROL
# include "OBCS.h"
#endif

c     == routine arguments ==

      integer mythid

c     == local variables ==

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer ntmp
      integer ivarindex
      integer iobcs

c     == end of interface ==

      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)
      jmin = 1-oly
      jmax = sny+oly
      imin = 1-olx
      imax = snx+olx

#ifdef ALLOW_OBCSN_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do i = imin,imax
                xx_obcsn0(i,k,bi,bj,iobcs) = 0. _d 0
                xx_obcsn1(i,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCSS_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do i = imin,imax
                xx_obcss0(i,k,bi,bj,iobcs) = 0. _d 0
                xx_obcss1(i,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCSW_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do j = jmin,jmax
                xx_obcsw0(j,k,bi,bj,iobcs) = 0. _d 0
                xx_obcsw1(j,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCSE_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do j = jmin,jmax
                xx_obcse0(j,k,bi,bj,iobcs) = 0. _d 0
                xx_obcse1(j,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCS_CONTROL
cih   Get matrices for reconstruction from barotropic-barclinic modes
CMM have not put in a good way to use/not use modes -- better to have
c    runtime option (e.g. defining filename=baro_invmodes.bin) 
c    for now hardcode with ECCO_CPPOPTION
#ifdef ALLOW_OBCS_CONTROL_MODES
CMM -- i want to use mdsreadvector but cannot as it assumes some
c       grid info and will increment record accordingly
c       so for now use this -- it works -- and fix in future....
c       should fix rel too - not universal!
        open(117, file='baro_invmodes.bin', access='direct',
     &     recl=2*WORDLENGTH*Nr*Nr, status='old')
        do j = 1,Nr
           read(117,rec=j) ((modesv(k,i,j), k=1,Nr), i=1,Nr)
        end do
#endif
CMM  AND HERE IS README PART OF CODE:
CMM  modesv should be orthonormal modes -- ideally solution of
c    planetary vertical mode equation using model mean dRho/dz
c    modesv is nr X nr X nr
c       dim one is z-space    
c       dim two is mode space 
c       dim three is the total depth for which this set of modes applies
c       so for example modesv(:,2,nr) will be the second mode
c       in z-space for the full model depth
c    It is important to note that the modes must be orthogonal
c       when weighted by dz!
c    i.e. if f_i(z) = mode i, sum_j(f_i(z_j)*f_j(z_j)*dz_j = delta_ij
c    first mode should also be constant in depth...barotropic
c    here is matlab code to ensure modes are orthonormal
c    ________________________________________________________
C     for k = 2:NZ;
C       iz = 1:k;% iz = vector of depth levels
C       NZ-k    % print out a progress number
C       % have to weight by dz!  Make a vector of fractional dz
C       zwt = DRF(iz)/sum(DRF(iz),1);
C       %ensure first mode is barotropic (constant in depth)
C       avm1=sum(modesv(iz,1,k).*zwt,1);
C       modesv(iz,1,k)=avm1;
C       %make all modes orthogonal weighted by delta z
C       % use gram-schmidt leaving first one the same
C       for mds = 1:k-1
C           R = sum((modesv(iz,mds,k).^2).*zwt,1);
C           R2 =(transpose(modesv(iz,mds,k).*zwt))*modesv(iz,mds+1:k,k)/R;
C           modesv(iz,mds+1:k,k) = modesv(iz,mds+1:k,k) - ...
C               modesv(iz,mds,k)*R2;
C         end
C       %All now orthognal, now normalize
C       for mds = 1:k
C          R = sqrt(sum((modesv(iz,mds,k).^2).*zwt,1));
C          if R < 1e-8
C            fprintf('Small R for mds = %2i and k = %2i\n',mds,k);
C            R = inf;
C          end
C          modesv(iz,mds,k) = modesv(iz,mds,k)./R;
C        end
C     end;% end of loop over depth level k
C     fid = fopen('baro_invmodes.bin','w','b');
C     fwrite(fid,modesv,'single' ); %or double depending on control prec
C     fclose(fid)
C     if 1 %plot first 5 modes for deepest case
C       figure
C       clf;plot(modesv(:,[1:5],NZ),RC(:));
C     end
C     if 1  % test orthogonality
C       % do whole matrix (need to include dz!)
C       k = NZ;
C       cmm=(transpose(squeeze(modesv(iz,iz,k)).*repmat(zwt,[1 k])))* ...
C             squeeze(modesv(iz,iz,k));
C       figure;imagesc(cmm);colorbar
C     end
C     save baro_invmodes modesv
c    ________________________________________________________

#endif


      return
      end