BoxAdj/code_ad/ctrl_init_obcs_variables.F

C
C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_init_obcs_variables.F,v 1.6 2004/12/04 17:12:00 heimbach Exp $
C $Name: checkpoint57a_pre $

#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
CMM(
      logical  exst
c      _RL  tmpmodes
      integer nz,nk 
CMM)

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