[MITgcm-devel] Inverted echo sounders
Matthew Mazloff
mmazloff at ucsd.edu
Wed Jan 19 15:37:51 EST 2011
Hi Martin,
Thanks for checking that in!
So I'm merging your changes and also ensuring backward compatibility
with ALLOW_CTRL_OBCS_BALANCE (even though I'm not sure the status of
that piece of code).
I'll check in my updates after they've been checked...
-Matt
ps> you put an #if instead of an #ifdef on line 94 of ctrl_get_gen_rec.F
On Jan 19, 2011, at 5:59 AM, Martin Losch wrote:
> Hi Matt,
> as you may have seen I have checked in the cost_obcs*.F files (I
> modified some comments to be a little more in line with general
> standards, I hope) from your directory. These were exactly the same
> files that I have been using for a while now, and never got around
> to checking in.
>
> This moring, I also checked in changes to ctrl_getobcs*.F that
> should be compatible with what you have in MITgcm_contrib/SOSE/
> BoxAdj/code_ad/ (basically got rid of exf_swapfields), but you need
> to merge them.
>
> Martin
>
> On Jan 18, 2011, at 8:50 PM, Matthew Mazloff wrote:
>
>> Hi Patrick,
>>
>> This code:
>> /u/gcmpack/MITgcm_contrib/SOSE/BoxAdj
>> is compatible with the latest code.
>>
>> You can compile it and run it with the files in 'input'. (gendata
>> is included too)
>>
>> The only cost term right now is cost_theta. To actually see the
>> true benefit of this code one must also enable ssh cost
>> constraints, which then make the barotropic sensitivity swamp the
>> baroclinic (cost_theta) sensitivity without this code. But I
>> wanted to keep the setup simple for now...let me know if you want
>> me to add this though
>>
>> Anyway, it works. The only place the code could use the help of an
>> expert is for subroutine
>> ctrl_init_obcs_variables.F
>> where the reading in of the modes is unrefined...perhaps you can
>> figure out a cleaner way to do it
>>
>> Let me know what you think -- it would be great to have this in the
>> main branch
>>
>> Thanks
>> 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
>>
>>
>> _______________________________________________
>> 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
More information about the MITgcm-devel
mailing list