[MITgcm-devel] Inverted echo sounders
Matthew Mazloff
mmazloff at ucsd.edu
Wed Jan 19 16:22:21 EST 2011
sorry -- OBCS.h has important changes too
-Matt
On Jan 19, 2011, at 1:15 PM, Matthew Mazloff wrote:
> Hi Martin,
>
> The important pieces of code to check in are
> ctrl_getobcse.F ctrl_getobcsn.F ctrl_getobcss.F ctrl_getobcsw.F
> ctrl_init_obcs_variables.F
>
> These all run properly with my setup. Model runs will be unchanged
> unless ALLOW_OBCS_CONTROL_MODES is defined.
>
> It would be nice to also add the simple additions in
> ecco_cost_final.F too
>
> Thanks!
> Matt
>
>
>
>
> On Jan 19, 2011, at 12:37 PM, Matthew Mazloff wrote:
>
>> 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
>>
>>
>> _______________________________________________
>> 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