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

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
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....
c       should fix rel too - not universal!
        open(117, file='baro_invmodes.bin', access='direct',
     &     recl=2*WORDLENGTH*Nr*Nr, status='old')
        do j = 1,Nr
           read(117,rec=j) ((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
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 =(transpose(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     end
C     if 1  % test orthogonality
C       % do whole matrix (need to include dz!)
C       k = NZ;
C       cmm=(transpose(squeeze(modesv(iz,iz,k)).*repmat(zwt,[1 k])))* ...
C             squeeze(modesv(iz,iz,k));
C       figure;imagesc(cmm);colorbar
C     end
C     save baro_invmodes modesv
c    ________________________________________________________

#endif


      return
      end


1	mmazloff	1.1	C
2	mmazloff	1.4	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	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			cih Get matrices for reconstruction from barotropic-barclinic modes
125			CMM have not put in a good way to use/not use modes -- better to have
126			c runtime option (e.g. defining filename=baro_invmodes.bin)
127			c for now hardcode with ECCO_CPPOPTION
128			#ifndef ALLOW_OBCS_CONTROL_MODES
129			CMM if this is not defined then don't use modes
130			c let z and k space be equal
131			do k = 1,nr
132			do j = 1,nr
133			do i = 1,nr
134			modesv(i,j,k) = 1. _d 0
135			enddo
136			enddo
137			enddo
138			#else
139			CMM -- i want to use mdsreadvector but cannot as it assumes some
140			c grid info and will increment record accordingly
141			c so for now use this -- it works -- and fix in future....
142	mmazloff	1.4	c should fix rel too - not universal!
143	mmazloff	1.1	open(117, file='baro_invmodes.bin', access='direct',
144			& recl=2WORDLENGTHNr*Nr, status='old')
145	mmazloff	1.2	do j = 1,Nr
146	mmazloff	1.3	read(117,rec=j) ((modesv(k,i,j), k=1,Nr), i=1,Nr)
147	mmazloff	1.1	end do
148			#endif
149			CMM AND HERE IS README PART OF CODE:
150			CMM modesv should be orthonormal modes -- ideally solution of
151			c planetary vertical mode equation using model mean dRho/dz
152			c modesv is nr X nr X nr
153			c dim one is z-space
154			c dim two is mode space
155			c dim three is the total depth for which this set of modes applies
156			c so for example modesv(:,2,nr) will be the second mode
157			c in z-space for the full model depth
158			c It is important to note that the modes must be orthogonal
159			c when weighted by dz!
160			c i.e. if f_i(z) = mode i, sum_j(f_i(z_j)f_j(z_j)dz_j = delta_ij
161			c first mode should also be constant in depth...barotropic
162			c here is matlab code to ensure modes are orthonormal
163			c ________________________________________________________
164			C for k = 2:NZ;
165			C iz = 1:k;% iz = vector of depth levels
166			C NZ-k % print out a progress number
167	mmazloff	1.4	C % have to weight by dz! Make a vector of fractional dz
168	mmazloff	1.1	C zwt = DRF(iz)/sum(DRF(iz),1);
169			C %ensure first mode is barotropic (constant in depth)
170			C avm1=sum(modesv(iz,1,k).*zwt,1);
171			C modesv(iz,1,k)=avm1;
172			C %make all modes orthogonal weighted by delta z
173			C % use gram-schmidt leaving first one the same
174			C for mds = 1:k-1
175			C R = sum((modesv(iz,mds,k).^2).*zwt,1);
176	mmazloff	1.4	C R2 =(transpose(modesv(iz,mds,k).zwt))modesv(iz,mds+1:k,k)/R;
177	mmazloff	1.1	C modesv(iz,mds+1:k,k) = modesv(iz,mds+1:k,k) - ...
178			C modesv(iz,mds,k)*R2;
179			C end
180			C %All now orthognal, now normalize
181			C for mds = 1:k
182			C R = sqrt(sum((modesv(iz,mds,k).^2).*zwt,1));
183			C if R < 1e-8
184	mmazloff	1.5	C fprintf('Small R for mds = %2i and k = %2i\n',mds,k);
185	mmazloff	1.1	C R = inf;
186			C end
187			C modesv(iz,mds,k) = modesv(iz,mds,k)./R;
188			C end
189			C end;% end of loop over depth level k
190			C fid = fopen('baro_invmodes.bin','w','b');
191			C fwrite(fid,modesv,'single');%or double depending on control prec
192			C fclose(fid)
193			C if 1 %plot first 5 modes for deepest case
194			C figure
195			C clf;plot(modesv(:,[1:5],NZ),RC(:));
196			C end
197			C if 1 % test orthogonality
198			C % do whole matrix (need to include dz!)
199			C k = NZ;
200	mmazloff	1.4	C cmm=(transpose(squeeze(modesv(iz,iz,k)).repmat(zwt,[1 k]))) ...
201	mmazloff	1.1	C squeeze(modesv(iz,iz,k));
202			C figure;imagesc(cmm);colorbar
203			C end
204			C save baro_invmodes modesv
205			c ________________________________________________________
206
207			#endif
208
209
210			return
211			end
212
213
214