[MITgcm-devel] Inverted echo sounders
Patrick Heimbach
heimbach at MIT.EDU
Tue Jan 18 07:47:41 EST 2011
Hi Matt,
thanks for the appended code.
The most promising way to get this checked in cleanly would be to
1. initially put the relevant code and code modifs in MITgcm_contrib/
2. provide an (as simple as possible) working reference setup as a test config.
Do you think you could do these two parts?
It would then be much simpler to merge it into the main branch
(works often, but not always...).
Cheers
-Patrick
On Jan 17, 2011, at 12:34 PM, Matthew Mazloff wrote:
> 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
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
---
Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
MIT | EAPS 54-1518 | 77 Massachusetts Ave | Cambridge MA 02139 USA
FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach
More information about the MITgcm-devel
mailing list