BoxAdj/code_ad/ctrl_init_obcs_variables.F

C
C $Header: /u/gcmpack/MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_init_obcs_variables.F,v 1.1 2011/01/18 19:33:08 mmazloff Exp $
C $Name:  $

#include "CTRL_CPPOPTIONS.h"

      subroutine ctrl_init_obcs_variables( mythid )

c     ==================================================================
c     SUBROUTINE ctrl_init_obcs_variables
c     ==================================================================
c
c     o Set parts of the vector of control variables and initialize the
c       rest to zero.
c
c     started: heimbach@mit.edu 25-Mar-2002
c
c     ==================================================================
c     SUBROUTINE ctrl_init_obcs_variables
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"

#include "ctrl.h"
#ifdef ALLOW_OBCS_CONTROL
# include "OBCS.h"
#endif

c     == routine arguments ==

      integer mythid

c     == local variables ==

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer ntmp
      integer ivarindex
      integer iobcs

c     == end of interface ==

      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)
      jmin = 1-oly
      jmax = sny+oly
      imin = 1-olx
      imax = snx+olx

#ifdef ALLOW_OBCSN_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do i = imin,imax
                xx_obcsn0(i,k,bi,bj,iobcs) = 0. _d 0
                xx_obcsn1(i,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCSS_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do i = imin,imax
                xx_obcss0(i,k,bi,bj,iobcs) = 0. _d 0
                xx_obcss1(i,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCSW_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do j = jmin,jmax
                xx_obcsw0(j,k,bi,bj,iobcs) = 0. _d 0
                xx_obcsw1(j,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCSE_CONTROL
      do iobcs = 1, nobcs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do k = 1,nr
              do j = jmin,jmax
                xx_obcse0(j,k,bi,bj,iobcs) = 0. _d 0
                xx_obcse1(j,k,bi,bj,iobcs) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
      enddo
#endif

#ifdef ALLOW_OBCS_CONTROL
        shiftvel(1) = 0. d0
        shiftvel(2) = 0. d0
#endif

#ifdef ALLOW_OBCS_CONTROL
cih   Get matrices for reconstruction from barotropic-barclinic modes
CMM have not put in a good way to use/not use modes -- better to have
c    runtime option (e.g. defining filename=baro_invmodes.bin) 
c    for now hardcode with ECCO_CPPOPTION
#ifndef ALLOW_OBCS_CONTROL_MODES
CMM  if this is not defined then don't use modes
c    let z and k space be equal
        do k = 1,nr
          do j = 1,nr
            do i = 1,nr
              modesv(i,j,k) = 1. _d 0
            enddo
          enddo
        enddo
#else
CMM -- i want to use mdsreadvector but cannot as it assumes some
c       grid info and will increment record accordingly
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 j = 1,Nr
           read(117,rec=nz) ((modesv(k,i,j), k=1,Nr), i=1,Nr)
        end do
#endif
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


1	mmazloff	1.1	C
2	mmazloff	1.2	C $Header: /u/gcmpack/MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_init_obcs_variables.F,v 1.1 2011/01/18 19:33:08 mmazloff Exp $
3			C $Name: $
4	mmazloff	1.1
5			#include "CTRL_CPPOPTIONS.h"
6
7			subroutine ctrl_init_obcs_variables( mythid )
8
9			c ==================================================================
10			c SUBROUTINE ctrl_init_obcs_variables
11			c ==================================================================
12			c
13			c o Set parts of the vector of control variables and initialize the
14			c rest to zero.
15			c
16			c started: heimbach@mit.edu 25-Mar-2002
17			c
18			c ==================================================================
19			c SUBROUTINE ctrl_init_obcs_variables
20			c ==================================================================
21
22			implicit none
23
24			c == global variables ==
25
26			#include "EEPARAMS.h"
27			#include "SIZE.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30
31			#include "ctrl.h"
32			#ifdef ALLOW_OBCS_CONTROL
33			# include "OBCS.h"
34			#endif
35
36			c == routine arguments ==
37
38			integer mythid
39
40			c == local variables ==
41
42			integer bi,bj
43			integer i,j,k
44			integer itlo,ithi
45			integer jtlo,jthi
46			integer jmin,jmax
47			integer imin,imax
48			integer ntmp
49			integer ivarindex
50			integer iobcs
51
52			c == end of interface ==
53
54			jtlo = mybylo(mythid)
55			jthi = mybyhi(mythid)
56			itlo = mybxlo(mythid)
57			ithi = mybxhi(mythid)
58			jmin = 1-oly
59			jmax = sny+oly
60			imin = 1-olx
61			imax = snx+olx
62
63			#ifdef ALLOW_OBCSN_CONTROL
64			do iobcs = 1, nobcs
65			do bj = jtlo,jthi
66			do bi = itlo,ithi
67			do k = 1,nr
68			do i = imin,imax
69			xx_obcsn0(i,k,bi,bj,iobcs) = 0. _d 0
70			xx_obcsn1(i,k,bi,bj,iobcs) = 0. _d 0
71			enddo
72			enddo
73			enddo
74			enddo
75			enddo
76			#endif
77
78			#ifdef ALLOW_OBCSS_CONTROL
79			do iobcs = 1, nobcs
80			do bj = jtlo,jthi
81			do bi = itlo,ithi
82			do k = 1,nr
83			do i = imin,imax
84			xx_obcss0(i,k,bi,bj,iobcs) = 0. _d 0
85			xx_obcss1(i,k,bi,bj,iobcs) = 0. _d 0
86			enddo
87			enddo
88			enddo
89			enddo
90			enddo
91			#endif
92
93			#ifdef ALLOW_OBCSW_CONTROL
94			do iobcs = 1, nobcs
95			do bj = jtlo,jthi
96			do bi = itlo,ithi
97			do k = 1,nr
98			do j = jmin,jmax
99			xx_obcsw0(j,k,bi,bj,iobcs) = 0. _d 0
100			xx_obcsw1(j,k,bi,bj,iobcs) = 0. _d 0
101			enddo
102			enddo
103			enddo
104			enddo
105			enddo
106			#endif
107
108			#ifdef ALLOW_OBCSE_CONTROL
109			do iobcs = 1, nobcs
110			do bj = jtlo,jthi
111			do bi = itlo,ithi
112			do k = 1,nr
113			do j = jmin,jmax
114			xx_obcse0(j,k,bi,bj,iobcs) = 0. _d 0
115			xx_obcse1(j,k,bi,bj,iobcs) = 0. _d 0
116			enddo
117			enddo
118			enddo
119			enddo
120			enddo
121			#endif
122
123			#ifdef ALLOW_OBCS_CONTROL
124			shiftvel(1) = 0. d0
125			shiftvel(2) = 0. d0
126			#endif
127
128			#ifdef ALLOW_OBCS_CONTROL
129			cih Get matrices for reconstruction from barotropic-barclinic modes
130			CMM have not put in a good way to use/not use modes -- better to have
131			c runtime option (e.g. defining filename=baro_invmodes.bin)
132			c for now hardcode with ECCO_CPPOPTION
133			#ifndef ALLOW_OBCS_CONTROL_MODES
134			CMM if this is not defined then don't use modes
135			c let z and k space be equal
136			do k = 1,nr
137			do j = 1,nr
138			do i = 1,nr
139			modesv(i,j,k) = 1. _d 0
140			enddo
141			enddo
142			enddo
143			#else
144			CMM -- i want to use mdsreadvector but cannot as it assumes some
145			c grid info and will increment record accordingly
146			c so for now use this -- it works -- and fix in future....
147			open(117, file='baro_invmodes.bin', access='direct',
148			& recl=2WORDLENGTHNr*Nr, status='old')
149	mmazloff	1.2	do j = 1,Nr
150			read(117,rec=nz) ((modesv(k,i,j), k=1,Nr), i=1,Nr)
151	mmazloff	1.1	end do
152			#endif
153			CMM AND HERE IS README PART OF CODE:
154			CMM modesv should be orthonormal modes -- ideally solution of
155			c planetary vertical mode equation using model mean dRho/dz
156			c modesv is nr X nr X nr
157			c dim one is z-space
158			c dim two is mode space
159			c dim three is the total depth for which this set of modes applies
160			c so for example modesv(:,2,nr) will be the second mode
161			c in z-space for the full model depth
162			c It is important to note that the modes must be orthogonal
163			c when weighted by dz!
164			c i.e. if f_i(z) = mode i, sum_j(f_i(z_j)f_j(z_j)dz_j = delta_ij
165			c first mode should also be constant in depth...barotropic
166			c here is matlab code to ensure modes are orthonormal
167			c ________________________________________________________
168			C for k = 2:NZ;
169			C iz = 1:k;% iz = vector of depth levels
170			C NZ-k % print out a progress number
171			C % have to weight by dz! Make a vector of fractional dz's
172			C zwt = DRF(iz)/sum(DRF(iz),1);
173			C %ensure first mode is barotropic (constant in depth)
174			C avm1=sum(modesv(iz,1,k).*zwt,1);
175			C modesv(iz,1,k)=avm1;
176			C %make all modes orthogonal weighted by delta z
177			C % use gram-schmidt leaving first one the same
178			C for mds = 1:k-1
179			C R = sum((modesv(iz,mds,k).^2).*zwt,1);
180			C R2 = (modesv(iz,mds,k).zwt)'modesv(iz,mds+1:k,k)/R;
181			C modesv(iz,mds+1:k,k) = modesv(iz,mds+1:k,k) - ...
182			C modesv(iz,mds,k)*R2;
183			C end
184			C %All now orthognal, now normalize
185			C for mds = 1:k
186			C R = sqrt(sum((modesv(iz,mds,k).^2).*zwt,1));
187			C if R < 1e-8
188			C fprintf('Small R!! for mds = %2i and k = %2i\n',mds,k);
189			C R = inf;
190			C end
191			C modesv(iz,mds,k) = modesv(iz,mds,k)./R;
192			C end
193			C end;% end of loop over depth level k
194			C fid = fopen('baro_invmodes.bin','w','b');
195			C fwrite(fid,modesv,'single');%or double depending on control prec
196			C fclose(fid)
197			C if 1 %plot first 5 modes for deepest case
198			C figure
199			C clf;plot(modesv(:,[1:5],NZ),RC(:));
200			C title('output modes at deepest point')
201			C end
202			C if 1 % test orthogonality
203			C % do whole matrix (need to include dz!)
204			C k = NZ;
205			C cmm = (squeeze(modesv(iz,iz,k)).repmat(zwt,[1 k]))' ...
206			C squeeze(modesv(iz,iz,k));
207			C figure;imagesc(cmm);colorbar
208			C title([num2str(k) ' mode orthogonality min, max diag: ' ...
209			C num2str(min(diag(cmm))) ', ' ...
210			C num2str(max(diag(cmm)))])
211			C end
212			C save baro_invmodes modesv
213			c ________________________________________________________
214
215			#endif
216
217
218			return
219			end
220
221
222