[MITgcm-devel] Inverted echo sounders
Matthew Mazloff
mmazloff at ucsd.edu
Tue Jan 18 11:41:16 EST 2011
Hi Patrick,
Sure -- I have a simple box model (64x64x8) this works for -- I'll
make sure this works with the latest code then check it in as a new
subdirectory under SOSE -- unless you have a better place.
-Matt
On Jan 18, 2011, at 4:47 AM, Patrick Heimbach wrote:
> 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
>
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list