[MITgcm-devel] Inverted echo sounders

[Matthew Mazloff]

Hi Martin,

That, but also the ability to control the obcs using a modal  
decomposition.  This is important when optimizing to both SSH and in  
situ obs as barotropic sensitivities from SSH cost overwhelms  
baroclinic sensitivities.  Thus the only way to optimize both is to  
separate the components and pre-condition accordingly.

The code changes involve
--  adding the mode matrix to ctrl.h
--  reading in the mode matrix in ctrl_init_obcs_variables.F
--  changing ctrl_getobcs[n,s,w,e].F to do the decomposition

Its fairly simple...it would be great if you wanted to check it in.

The only issue is one likely would want to fix the reading in of the  
modes to be more consistent with other MITgcm I/O practices.  I will  
attach this piece of the code below.  Also, in this part of the code I  
put the only "Read Me" aspect...so a lot of (hopefully useful) comments

Let me know what you think,
Matt




This is starting at line ~131 of ctrl_init_obcs_variables.F




#ifdef ALLOW_OBCS_CONTROL
cih   Get matrices for reconstruction from barotropic-barclinic modes
CMM( Likely want to define filename=baro_invmodes.bin at runtime
c    using data.ctrl...hardcoded with no checks for now...)
c    if baro_invmodes.bin is given then will use modes
c      otherwise z and k space equal -- no modes
c      really need more checks on this though!!!!

          do k = 1,nr
            do j = 1,nr
              do i = 1,nr
                modesv(i,j,k) = 1. _d 0
              enddo
            enddo
          enddo
          inquire( file='baro_invmodes.bin', exist=exst )
          IF ( exst ) THEN

CMM -- i want to use mdsreadvector but cannot as it assumes some
c       grid info and will increment record based on processor..
c       so for now use this -- it works -- and fix in future....

         open(117, file='baro_invmodes.bin', access='direct',
      &     recl=2*WORDLENGTH*Nr*Nr, status='old')
         do nz = 1,Nr
            read(117,rec=nz) ((modesv(k,nk,nz), k=1,Nr), nk=1,Nr)
         end do

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's
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 = (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       title('output modes at deepest point')
C     end
C     if 1  % test orthogonality
C       % do whole matrix (need to include dz!)
C       k = NZ;
C       cmm = (squeeze(modesv(iz,iz,k)).*repmat(zwt,[1 k]))'* ...
C             squeeze(modesv(iz,iz,k));
C       figure;imagesc(cmm);colorbar
C       title([num2str(k) ' mode orthogonality min, max diag: ' ...
C              num2str(min(diag(cmm))) ', ' ...
C              num2str(max(diag(cmm)))])
C     end
C     save baro_invmodes modesv
c    ________________________________________________________

#endif


       return
       end



~





On Jan 17, 2011, at 12:29 AM, Martin Losch wrote:

> Matt,
>
> I have obcs code to check in: regularization terms for the ECCO cost  
> function, which I got from you (ecco_cost_obcs[n,s,w,e].F). Is that  
> the change that you would like to see in the repository?
>
> Martin
>
> On Jan 15, 2011, at 12:27 AM, Matthew Mazloff wrote:
>
>> ps> would like to get obcs mode code checked in too...
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel