<div dir="ltr">Thanks for the response, Martin. Using myBxLo(myThid),myBxHi(myThid) and myByLo(myThid),myByHi(myThid) instead of 1,nSx and 1,nSy doesn't make a difference. The code I sent you happened to use the latter because I suspected that it would make a difference too, but it doesn't. Ngrid definitely has the correct value. Also, I figured that I should show you an example input to the penta-diagonal solver subroutine.<div><br></div><div>Here is the forcing term (which is variable B in the subroutine code I sent you):</div><div><img src="cid:ii_ktu95lu71" alt="image.png" width="499" height="365"><br></div><div>I've n-ple-checked that this is correct. Here is the diagonal entry of the penta-diagonal matrix (which is variable D in the subroutine code I sent you):<div><img src="cid:ii_ktu9vmjo2" alt="image.png" width="499" height="377"><br></div><div>You can see that this is also correct, as it should correlate highly with the seafloor depths.</div><div><br></div><div>You mentioned doing a run without blank tiles. I didn't know this was an option. It looked to me like I just needed to comment out blankList in data.exch2. I tried this and I still see the "hole" over Asia. What else do I need to do to try this? I suspect that the "hole" isn't the main problem, though. My main problem is my penta-diagonal matrix solver subroutine. I can't just invert the matrix using some LU decomposition subroutine (too much memory). So that leaves me stuck because I don't expect anyone to know how to correctly program the very specific type of matrix solver I need in Fortran...</div><div>Cheers,</div><div>David</div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Sep 21, 2021 at 5:17 AM Martin Losch <<a href="mailto:Martin.Losch@awi.de" target="_blank">Martin.Losch@awi.de</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,204);padding-left:1ex">Hi David,<br>
<br>
have you tried to do this in a run without blank tiles? That way the “hole” over Asia has some reasonable coordinates (and values=0 instead of NaN).<br>
Also you may want to use myBxLo(myThid),myBxHi(myThid) instead of 1,nSx and similar for y?<br>
<br>
I am assuming that Ngrid has the correct value. Do you need the overlap here? You I,J loop ranges include the overlaps, but are these required for the solution (as they are just duplicates from the neighboring tiles). If you do need overlaps, then the code implies that you only need one and not OLx/OLy?<br>
<br>
Martin<br>
<br>
> On 14. Sep 2021, at 18:01, David Trossman <<a href="mailto:david.s.trossman@gmail.com" target="_blank">david.s.trossman@gmail.com</a>> wrote:<br>
> <br>
> Hi all,<br>
> I've been banging my head against the wall with this problem for so long I may have too much brain damage to come up with a solution.<br>
> <br>
> The problem is this: I want to solve some equations to solve for the magnetic field associated with the ocean circulation, and these equations involve the inverse of a penta-diagonal matrix. I translated what would be the following in Matlab (I know, I know...):<br>
> <br>
> difFluxw = spdiags([-Cw(:),[Cw(2:end);nan]], [0, -1], Ngrid,Ngrid);<br>
> difFluxe = spdiags([-Ce(:),[nan;Ce(1:end-1)]], [0, 1], Ngrid,Ngrid);<br>
> difFluxs = spdiags([-Cs(:),[Cs(1+Nlon:end);nan*ones(Nlon,1)]],[0, -Nlon], Ngrid,Ngrid);<br>
> difFluxn = spdiags([-Cn(:),[nan*ones(Nlon,1);Cn(1:end-Nlon)]],[0, Nlon], Ngrid,Ngrid);<br>
> advFluxw = spdiags([-aCw(:),-[aCw(2:end);nan]], [0, -1], Ngrid,Ngrid);<br>
> advFluxe = spdiags([ aCe(:), [nan;aCe(1:end-1)]], [0, 1], Ngrid,Ngrid);<br>
> advFluxs = spdiags([-aCs(:),-[aCs(1+Nlon:end);nan*ones(Nlon,1)]],[0, -Nlon], Ngrid,Ngrid);<br>
> advFluxn = spdiags([ aCn(:), [nan*ones(Nlon,1);aCn(1:end-Nlon)]],[0, Nlon], Ngrid,Ngrid);<br>
> Fluxw = difFluxw + advFluxw; <br>
> Fluxe = difFluxe + advFluxe; <br>
> Fluxs = difFluxs + advFluxs; <br>
> Fluxn = difFluxn + advFluxn; <br>
> bjunk=find(lon==min(lon(:)) & lat>min(lat(:)) & lat<max(lat(:))); <br>
> Fluxw(bjunk,:) = 0;<br>
> Fluxw = Fluxw + sparse(bjunk,bjunk, -Cw(bjunk), Ngrid,Ngrid);<br>
> Fluxw = Fluxw + sparse(bjunk,bjunk-1+Nlon, Cw(bjunk), Ngrid,Ngrid);<br>
> % flux through east side of cell:<br>
> bjunk=find(lon==max(lon(:)) & lat > min(lat(:)) & lat<max(lat(:)));<br>
> Fluxe(bjunk,:) = 0;<br>
> Fluxe = Fluxe + sparse(bjunk,bjunk, -Ce(bjunk), Ngrid,Ngrid);<br>
> Fluxe = Fluxe + sparse(bjunk,bjunk+1-Nlon, Ce(bjunk),Ngrid,Ngrid);<br>
> integHcoef = spdiags(Hcoef(:).*vol(:), 0, Ngrid,Ngrid);<br>
> metrik.Lfluxes = Fluxw + Fluxe + Fluxs + Fluxn ;<br>
> metrik.L = metrik.Lfluxes + integHcoef;<br>
> <br>
> into Fortran. Nevermind what the variables mean; all that's important is that the variable metrik.L is a pentadiagonal matrix that I save in Fortran as two off-diagonal arrays that are Ngrid-Nlon in length, two off-diagonal arrays that are Ngrid-1 in length, and a diagonal array of length Ngrid. I've written a pent-diagonal matrix solver for this type of situation as follows:<br>
> <br>
> SUBROUTINE PENTA(myThid,Ngrid,Nlon,E0,A0,D0,C0,F0,B0,X0)<br>
> <br>
> c RESULTS: matrix has 5 bands, EADCF, with D being the main diagonal,<br>
> c E and A are the lower diagonals, and C and F are the upper diagonals.<br>
> <br>
> c E is defined for rows i = n1:Ngrid, but is defined as E(1) to E(Ngrid-Nlon)<br>
> c A is defined for rows i = 2:Ngrid, but is defined as A(1) to A(Ngrid-1)<br>
> c D is defined for rows i = 1:Ngrid<br>
> c C is defined for rows i = 1:Ngrid-1, but the last element isn't used<br>
> c F is defined for rows i = 1:(Ngrid-Nlon), but the last Nlon-1 elements aren't used<br>
> <br>
> c B is the right-hand side<br>
> c X is the solution vector<br>
> <br>
> IMPLICIT NONE<br>
> #include "SIZE.h"<br>
> #include "EEPARAMS.h"<br>
> #include "GRID.h"<br>
> #include "SURFACE.h"<br>
> <br>
> INTEGER, intent(in) :: myThid<br>
> INTEGER, intent(in) :: Ngrid<br>
> INTEGER, intent(in) :: Nlon<br>
> _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),<br>
> & intent(in) :: E0<br>
> _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),<br>
> & intent(in) :: A0<br>
> _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),<br>
> & intent(in) :: D0<br>
> _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),<br>
> & intent(in) :: C0<br>
> _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),<br>
> & intent(in) :: F0<br>
> _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),<br>
> & intent(out) :: B0<br>
> _RL, dimension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nSx,1:nSy),<br>
> & intent(out) :: X0<br>
> <br>
> INTEGER bi,bj<br>
> INTEGER I,J,IJ<br>
> INTEGER Nn<br>
> _RL XMULT<br>
> _RL, dimension(1:Ngrid*nSx*nSy) :: E<br>
> _RL, dimension(1:Ngrid*nSx*nSy) :: A<br>
> _RL, dimension(1:Ngrid*nSx*nSy) :: D<br>
> _RL, dimension(1:Ngrid*nSx*nSy) :: C<br>
> _RL, dimension(1:Ngrid*nSx*nSy) :: F<br>
> _RL, dimension(1:Ngrid*nSx*nSy) :: B<br>
> _RL, dimension(1:Ngrid*nSx*nSy) :: X<br>
> <br>
> Nn=Ngrid*nSx*nSy<br>
> IJ=0<br>
> DO bj=1,nSy<br>
> DO bi=1,nSx<br>
> DO I = 1-OLx,sNx+OLx<br>
> DO J = 1-OLy,sNy+OLy<br>
> IJ=IJ+1<br>
> E(IJ)=E0(I,J,bi,bj)<br>
> A(IJ)=A0(I,J,bi,bj)<br>
> D(IJ)=D0(I,J,bi,bj)<br>
> C(IJ)=C0(I,J,bi,bj)<br>
> F(IJ)=F0(I,J,bi,bj)<br>
> ENDDO<br>
> ENDDO<br>
> ENDDO<br>
> ENDDO<br>
> DO IJ = 2,Nn-1<br>
> XMULT = A(IJ-1)/D(IJ-1)<br>
> D(IJ) = D(IJ) - XMULT*C(IJ-1)<br>
> C(IJ) = C(IJ) - XMULT*F(IJ-1)<br>
> B(IJ) = B(IJ) - XMULT*B(IJ-1)<br>
> XMULT = E(IJ-1)/D(IJ-1)<br>
> A(IJ) = A(IJ) - XMULT*C(IJ-1)<br>
> D(IJ+1) = D(IJ+1) - XMULT*F(IJ-1)<br>
> B(IJ+1) = B(IJ+1) - XMULT*B(IJ-1)<br>
> ENDDO<br>
> XMULT = A(Nn-1)/D(Nn-1)<br>
> D(Nn) = D(Nn) - XMULT*C(Nn-1)<br>
> X(Nn) = (B(Nn) - XMULT*B(Nn-1))/D(Nn)<br>
> DO IJ = Nn-1,Nn-Nlon+1,-1<br>
> X(IJ) = (B(IJ) - C(IJ)*X(IJ+1))/D(IJ)<br>
> ENDDO<br>
> DO IJ = Nn-Nlon,1,-1<br>
> X(IJ) = (B(IJ) - F(IJ)*X(IJ+n1) -<br>
> & C(IJ)*X(IJ+1))/D(IJ)<br>
> ENDDO<br>
> IJ=0<br>
> DO bj=1,nSy<br>
> DO bi=1,nSx<br>
> DO I = 1-OLx,sNx+OLx<br>
> DO J = 1-OLy,sNy+OLy<br>
> IJ=IJ+1<br>
> B0(I,J,bi,bj)=B(IJ)<br>
> X0(I,J,bi,bj)=X(IJ)<br>
> ENDDO<br>
> ENDDO<br>
> ENDDO<br>
> ENDDO<br>
> RETURN<br>
> END<br>
> <br>
> X0 is what I want, but the above code isn't correct. I get something like the following for the magnetic field:<br>
> <image.png><br>
> It seems like I'm not considering an ordering for the tiling correctly, or the pent-diagonal matrix solver isn't quite correct in other ways (although, I'm adapting it from an existing code with two off-diagonal arrays of length Ngrid-2, two off-diagonal arrays of length Ngrid-1, and a diagonal array of length Ngrid: <a href="http://www.math.uakron.edu/~kreider/anpde/penta.f" rel="noreferrer" target="_blank">http://www.math.uakron.edu/~kreider/anpde/penta.f</a>). It's also possible that I need to interpolate over the "hole" (of NaNs) over Asia in the MITgcm. I recognize that you don't know whether my inputs to the subroutine are correct, but just assume they are.<br>
> <br>
> One possible approach to figuring out what's wrong is to interpolate from the LLC grid to a regularly spaced lat-lon grid, but I'm not sure if this capability exists inline as the model runs. No, I don't want to do this offline using averaged output because I'm trying to make this calculation adjoint-able (plus, a function of the average of each field is not the same as the average of a function of each field). I figured that someone on this list would be able to point out what I'm doing incorrectly in the above Fortran code or knows how to do the interpolation to a regular lat-lon grid, but maybe I'm asking for too much... In any case, thanks for any help.<br>
> -David<br>
> _______________________________________________<br>
> MITgcm-support mailing list<br>
> <a href="mailto:MITgcm-support@mitgcm.org" target="_blank">MITgcm-support@mitgcm.org</a><br>
> <a href="http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support" rel="noreferrer" target="_blank">http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support</a><br>
<br>
_______________________________________________<br>
MITgcm-support mailing list<br>
<a href="mailto:MITgcm-support@mitgcm.org" target="_blank">MITgcm-support@mitgcm.org</a><br>
<a href="http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support" rel="noreferrer" target="_blank">http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support</a><br>
</blockquote></div>