/[MITgcm]/MITgcm/pkg/ctrl/ctrl_set_pack_xyz.F
ViewVC logotype

Diff of /MITgcm/pkg/ctrl/ctrl_set_pack_xyz.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by heimbach, Sat Jul 13 02:47:32 2002 UTC revision 1.19 by jmc, Tue Oct 9 00:00:01 2007 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "CTRL_CPPOPTIONS.h"  #include "CTRL_CPPOPTIONS.h"
5    
   
6        subroutine ctrl_set_pack_xyz(        subroutine ctrl_set_pack_xyz(
7       &     cunit, ivartype, fname, masktype,       &     cunit, ivartype, fname, masktype, weighttype,
8       &     weightfld, lxxadxx, mythid)       &     weightfld, lxxadxx, mythid)
9    
10  c     ==================================================================  c     ==================================================================
# Line 13  c Line 14  c
14  c     o Compress the control vector such that only ocean points are  c     o Compress the control vector such that only ocean points are
15  c       written to file.  c       written to file.
16  c  c
17    c     o Use a more precise nondimensionalization that depends on (x,y)
18    c       Added weighttype to the argument list so that I can geographically
19    c       vary the nondimensionalization.
20    c       gebbie@mit.edu, 18-Mar-2003
21    c
22  c     ==================================================================  c     ==================================================================
23    
24        implicit none        implicit none
# Line 25  c     == global variables == Line 31  c     == global variables ==
31  #include "GRID.h"  #include "GRID.h"
32    
33  #include "ctrl.h"  #include "ctrl.h"
 #include "cost.h"  
   
 #ifdef ALLOW_ECCO_OPTIMIZATION  
34  #include "optim.h"  #include "optim.h"
 #endif  
35    
36  c     == routine arguments ==  c     == routine arguments ==
37    
38        integer cunit        integer cunit
39        integer ivartype        integer ivartype
40        character*( 80) fname        character*( 80) fname
41        character*  (5) masktype        character*(  9) masktype
42          character*( 80) weighttype
43        _RL     weightfld( nr,nsx,nsy )        _RL     weightfld( nr,nsx,nsy )
44        logical lxxadxx        logical lxxadxx
45        integer mythid        integer mythid
46    
47    #ifndef EXCLUDE_CTRL_PACK
48  c     == local variables ==  c     == local variables ==
49    
 #ifndef ALLOW_ECCO_OPTIMIZATION  
       integer optimcycle  
 #endif  
   
50        integer bi,bj        integer bi,bj
51        integer ip,jp        integer ip,jp
52        integer i,j,k        integer i,j,k
# Line 60  c     == local variables == Line 60  c     == local variables ==
60    
61        integer cbuffindex        integer cbuffindex
62    
       _RL     cbuff    ( snx*nsx*npx*sny*nsy*npy )  
63        _RL     globmsk  ( snx,nsx,npx,sny,nsy,npy,nr )        _RL     globmsk  ( snx,nsx,npx,sny,nsy,npy,nr )
64        _RL     globfld3d( snx,nsx,npx,sny,nsy,npy,nr )        _RL     globfld3d( snx,nsx,npx,sny,nsy,npy,nr )
65    #ifdef CTRL_PACK_PRECISE
66          _RL   weightfld3d( snx,nsx,npx,sny,nsy,npy,nr )
67    #endif
68          real*4 cbuff      ( snx*nsx*npx*sny*nsy*npy )
69          real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
70          real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )
71    
72          character*(80) weightname
73    
74          _RL delZnorm
75          integer reclen, irectrue
76          integer cunit2, cunit3
77          character*(80) cfile2, cfile3
78    
79  c     == external ==  c     == external ==
80    
# Line 71  c     == external == Line 83  c     == external ==
83    
84  c     == end of interface ==  c     == end of interface ==
85    
 #ifndef ALLOW_ECCO_OPTIMIZATION  
       optimcycle = 0  
 #endif  
   
86        jtlo = 1        jtlo = 1
87        jthi = nsy        jthi = nsy
88        itlo = 1        itlo = 1
# Line 84  c     == end of interface == Line 92  c     == end of interface ==
92        imin = 1        imin = 1
93        imax = snx        imax = snx
94    
95    #ifdef CTRL_DELZNORM
96          delZnorm = 0.
97          do k = 1, Nr
98             delZnorm = delZnorm + delR(k)/FLOAT(Nr)
99          enddo
100    #endif
101    
102  c     Initialise temporary file  c     Initialise temporary file
103        do k = 1,nr        do k = 1,nr
104           do jp = 1,nPy           do jp = 1,nPy
# Line 92  c     Initialise temporary file Line 107  c     Initialise temporary file
107                    do ip = 1,nPx                    do ip = 1,nPx
108                       do bi = itlo,ithi                       do bi = itlo,ithi
109                          do i = imin,imax                          do i = imin,imax
110                             globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0                             globfld3d  (i,bi,ip,j,bj,jp,k) = 0. _d 0
111                             globmsk  (i,bi,ip,j,bj,jp,k) = 0. _d 0                             globmsk    (i,bi,ip,j,bj,jp,k) = 0. _d 0
112                               globfldtmp2(i,bi,ip,j,bj,jp)   = 0.
113                               globfldtmp3(i,bi,ip,j,bj,jp)   = 0.
114                          enddo                          enddo
115                       enddo                       enddo
116                    enddo                    enddo
# Line 105  c     Initialise temporary file Line 122  c     Initialise temporary file
122  c--   Only the master thread will do I/O.  c--   Only the master thread will do I/O.
123        _BEGIN_MASTER( mythid )        _BEGIN_MASTER( mythid )
124    
125        call MDSREADFIELD_3D_GL(        if ( doPackDiag ) then
126             write(cfile2(1:80),'(80a)') ' '
127             write(cfile3(1:80),'(80a)') ' '
128             if ( lxxadxx ) then
129                write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
130         &           'diag_pack_nonout_ctrl_',
131         &           ivartype, '_', optimcycle, '.bin'
132                write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
133         &           'diag_pack_dimout_ctrl_',
134         &           ivartype, '_', optimcycle, '.bin'
135             else
136                write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
137         &           'diag_pack_nonout_grad_',
138         &           ivartype, '_', optimcycle, '.bin'
139                write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
140         &           'diag_pack_dimout_grad_',
141         &           ivartype, '_', optimcycle, '.bin'
142             endif
143    
144             reclen = FLOAT(snx*nsx*npx*sny*nsy*npy*4)
145             call mdsfindunit( cunit2, mythid )
146             open( cunit2, file=cfile2, status='unknown',
147         &        access='direct', recl=reclen )
148             call mdsfindunit( cunit3, mythid )
149             open( cunit3, file=cfile3, status='unknown',
150         &        access='direct', recl=reclen )
151          endif
152    
153    #ifdef CTRL_PACK_PRECISE
154          il=ilnblnk( weighttype)
155          write(weightname(1:80),'(80a)') ' '
156          write(weightname(1:80),'(a)') weighttype(1:il)
157    
158          call MDSREADFIELD_3D_GL(
159         &     weightname, ctrlprec, 'RL',
160         &     Nr, weightfld3d, 1, mythid)
161    #endif
162    
163          call MDSREADFIELD_3D_GL(
164       &     masktype, ctrlprec, 'RL',       &     masktype, ctrlprec, 'RL',
165       &     Nr, globmsk, 1, mythid)       &     Nr, globmsk, 1, mythid)
166    
# Line 114  c--   Only the master thread will do I/O Line 169  c--   Only the master thread will do I/O
169           call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',           call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
170       &        Nr, globfld3d, irec, mythid)       &        Nr, globfld3d, irec, mythid)
171    
172    #ifndef ALLOW_ADMTLM
173           write(cunit) ncvarindex(ivartype)           write(cunit) ncvarindex(ivartype)
174           write(cunit) 1           write(cunit) 1
175           write(cunit) 1           write(cunit) 1
176    #endif
177           do k = 1, nr           do k = 1, nr
178             irectrue = (irec-1)*nr + k
179                if ( doZscalePack ) then
180                   delZnorm = (delR(1)/delR(k))**delZexp
181                else
182                   delZnorm = 1. _d 0
183                endif
184              cbuffindex = 0              cbuffindex = 0
185              do jp = 1,nPy              do jp = 1,nPy
186               do bj = jtlo,jthi               do bj = jtlo,jthi
# Line 125  c--   Only the master thread will do I/O Line 188  c--   Only the master thread will do I/O
188                 do ip = 1,nPx                 do ip = 1,nPx
189                  do bi = itlo,ithi                  do bi = itlo,ithi
190                   do i = imin,imax                   do i = imin,imax
191                    if (globmsk(i,bi,ip,j,bj,jp,k)  .ne. 0. ) then                    if (globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
192                       cbuffindex = cbuffindex + 1                       cbuffindex = cbuffindex + 1
193    cph(
194                         globfldtmp3(i,bi,ip,j,bj,jp) =
195         &                    globfld3d(i,bi,ip,j,bj,jp,k)
196    cph)
197    #ifndef ALLOW_SMOOTH_CORREL3D
198  #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO  #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
199                       if (lxxadxx) then                       if (lxxadxx) then
200                          cbuff(cbuffindex) =                          cbuff(cbuffindex) = 1/delZnorm
201       &                       globfld3d(i,bi,ip,j,bj,jp,k) *       &                       * globfld3d(i,bi,ip,j,bj,jp,k)
202       &                       sqrt(weightfld(k,bi,bj))  # ifdef CTRL_PACK_PRECISE
203         &                       * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
204    # else
205         &                       * sqrt(weightfld(k,bi,bj))
206    # endif
207                       else                       else
208                            cbuff(cbuffindex) = delZnorm
209         &                       * globfld3d(i,bi,ip,j,bj,jp,k)
210    # ifdef CTRL_PACK_PRECISE
211         &                       / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
212    # else
213         &                       / sqrt(weightfld(k,bi,bj))
214    # endif
215                         endif
216  cph(  cph(
217                          print *, 'ph-nondim bef. ', k, j, i,                       globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
      &                       globfld3d(i,bi,ip,j,bj,jp,k),  
      &                       weightfld(k,bi,bj)  
 cph)  
                         cbuff(cbuffindex) =  
      &                       globfld3d(i,bi,ip,j,bj,jp,k) /  
      &                       sqrt(weightfld(k,bi,bj))  
 cph(  
                         write(6,'(A,4I5,F10.2)'), 'ph-nondim aft. ',  
      &                       k, j, i, cbuffindex,  
      &                       cbuff(cbuffindex)  
218  cph)  cph)
219                       endif  #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
220  #else                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
221    #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
222    #else  /* ALLOW_SMOOTH_CORREL3D */
223                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
224    #endif /* ALLOW_SMOOTH_CORREL3D */
225    #ifdef ALLOW_ADMTLM
226                         nveccount = nveccount + 1
227                         phtmpadmtlm(nveccount) = cbuff(cbuffindex)
228  #endif  #endif
229                    endif                    endif
230                   enddo                   enddo
# Line 159  cph) Line 235  cph)
235              enddo              enddo
236  c           --> check cbuffindex.  c           --> check cbuffindex.
237              if ( cbuffindex .gt. 0) then              if ( cbuffindex .gt. 0) then
238    #ifndef ALLOW_ADMTLM
239                 write(cunit) cbuffindex                 write(cunit) cbuffindex
240                 write(cunit) k                 write(cunit) k
241    cph#endif
242                 write(cunit) (cbuff(ii), ii=1,cbuffindex)                 write(cunit) (cbuff(ii), ii=1,cbuffindex)
243    #endif
244                endif
245    c
246                if ( doPackDiag ) then
247                   write(cunit2,rec=irectrue) globfldtmp2
248                   write(cunit3,rec=irectrue) globfldtmp3
249              endif              endif
250    c
251           enddo           enddo
252  c  c
253  c     -- end of irec loop --  c     -- end of irec loop --
254        enddo        enddo
255    
256          if ( doPackDiag ) then
257             close ( cunit2 )
258             close ( cunit3 )
259          endif
260    
261        _END_MASTER( mythid )        _END_MASTER( mythid )
262    
263    #endif
264    
265        return        return
266        end        end
267    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22