/[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.1.4.1 by heimbach, Fri Jul 12 15:43:54 2002 UTC revision 1.28 by gforget, Tue Jun 16 20:20:26 2015 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "CTRL_CPPOPTIONS.h"  #include "CTRL_OPTIONS.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    # ifndef ALLOW_PACKUNPACK_METHOD2
49  c     == local variables ==  c     == local variables ==
50    
 #ifndef ALLOW_ECCO_OPTIMIZATION  
       integer optimcycle  
 #endif  
   
51        integer bi,bj        integer bi,bj
52        integer ip,jp        integer ip,jp
53        integer i,j,k        integer i,j,k
54        integer ii        integer ii
       integer il  
55        integer irec        integer irec
56        integer itlo,ithi        integer itlo,ithi
57        integer jtlo,jthi        integer jtlo,jthi
# 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          integer il
67          character*(80) weightname
68          _RL   weightfld3d( snx,nsx,npx,sny,nsy,npy,nr )
69    #endif
70          real*4 cbuff      ( snx*nsx*npx*sny*nsy*npy )
71          real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
72          real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )
73    
74          _RL delZnorm
75          integer reclen, irectrue
76          integer cunit2, cunit3
77          character*(80) cfile2, cfile3
78    
79          LOGICAL doPackOld
80    
81  c     == external ==  c     == external ==
82    
# Line 71  c     == external == Line 85  c     == external ==
85    
86  c     == end of interface ==  c     == end of interface ==
87    
 #ifndef ALLOW_ECCO_OPTIMIZATION  
       optimcycle = 0  
 #endif  
   
88        jtlo = 1        jtlo = 1
89        jthi = nsy        jthi = nsy
90        itlo = 1        itlo = 1
# Line 84  c     == end of interface == Line 94  c     == end of interface ==
94        imin = 1        imin = 1
95        imax = snx        imax = snx
96    
97          doPackOld = (.NOT.ctrlSmoothCorrel3D).AND.(.NOT.ctrlUseGen)
98    
99    #ifdef CTRL_DELZNORM
100          delZnorm = 0.
101          do k = 1, Nr
102             delZnorm = delZnorm + delR(k)/FLOAT(Nr)
103          enddo
104    #endif
105    
106  c     Initialise temporary file  c     Initialise temporary file
107        do k = 1,nr        do k = 1,nr
108           do jp = 1,nPy         do jp = 1,nPy
109              do bj = jtlo,jthi          do bj = jtlo,jthi
110                 do j = jmin,jmax           do j = jmin,jmax
111                    do ip = 1,nPx            do ip = 1,nPx
112                       do bi = itlo,ithi             do bi = itlo,ithi
113                          do i = imin,imax              do i = imin,imax
114                             globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0               globfld3d  (i,bi,ip,j,bj,jp,k) = 0. _d 0
115                             globmsk  (i,bi,ip,j,bj,jp,k) = 0. _d 0               globmsk    (i,bi,ip,j,bj,jp,k) = 0. _d 0
116                          enddo               globfldtmp2(i,bi,ip,j,bj,jp)   = 0. _d 0
117                       enddo               globfldtmp3(i,bi,ip,j,bj,jp)   = 0. _d 0
                   enddo  
                enddo  
118              enddo              enddo
119               enddo
120              enddo
121           enddo           enddo
122            enddo
123           enddo
124        enddo        enddo
125    
126  c--   Only the master thread will do I/O.  c--   Only the master thread will do I/O.
127        _BEGIN_MASTER( mythid )        _BEGIN_MASTER( mythid )
128    
129        call MDSREADFIELD_3D_GL(        if ( doPackDiag ) then
130             write(cfile2(1:80),'(80a)') ' '
131             write(cfile3(1:80),'(80a)') ' '
132             if ( lxxadxx ) then
133                write(cfile2(1:80),'(a,I3.3,a,I4.4,a)')
134         &           'diag_pack_nonout_ctrl_',
135         &           ivartype, '_', optimcycle, '.bin'
136                write(cfile3(1:80),'(a,I3.3,a,I4.4,a)')
137         &           'diag_pack_dimout_ctrl_',
138         &           ivartype, '_', optimcycle, '.bin'
139             else
140                write(cfile2(1:80),'(a,I3.3,a,I4.4,a)')
141         &           'diag_pack_nonout_grad_',
142         &           ivartype, '_', optimcycle, '.bin'
143                write(cfile3(1:80),'(a,I3.3,a,I4.4,a)')
144         &           'diag_pack_dimout_grad_',
145         &           ivartype, '_', optimcycle, '.bin'
146             endif
147    
148             reclen = FLOAT(snx*nsx*npx*sny*nsy*npy*4)
149             call mdsfindunit( cunit2, mythid )
150             open( cunit2, file=cfile2, status='unknown',
151         &        access='direct', recl=reclen )
152             call mdsfindunit( cunit3, mythid )
153             open( cunit3, file=cfile3, status='unknown',
154         &        access='direct', recl=reclen )
155          endif
156    
157    #ifdef CTRL_PACK_PRECISE
158          if (weighttype.NE.' ') then
159           il=ilnblnk( weighttype)
160           write(weightname(1:80),'(80a)') ' '
161           write(weightname(1:80),'(a)') weighttype(1:il)
162           call MDSREADFIELD_3D_GL(
163         &     weightname, ctrlprec, 'RL',
164         &     Nr, weightfld3d, 1, mythid)
165          else
166           do k = 1,nr
167            do jp = 1,nPy
168             do bj = jtlo,jthi
169              do j = jmin,jmax
170               do ip = 1,nPx
171                do bi = itlo,ithi
172                 do i = imin,imax
173                  weightfld3d(i,bi,ip,j,bj,jp,k) = 1. _d 0
174                 enddo
175                enddo
176               enddo
177              enddo
178             enddo
179            enddo
180           enddo
181          endif
182    #endif
183    
184          call MDSREADFIELD_3D_GL(
185       &     masktype, ctrlprec, 'RL',       &     masktype, ctrlprec, 'RL',
186       &     Nr, globmsk, 1, mythid)       &     Nr, globmsk, 1, mythid)
187    
# Line 114  c--   Only the master thread will do I/O Line 190  c--   Only the master thread will do I/O
190           call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',           call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
191       &        Nr, globfld3d, irec, mythid)       &        Nr, globfld3d, irec, mythid)
192    
193    #ifndef ALLOW_ADMTLM
194           write(cunit) ncvarindex(ivartype)           write(cunit) ncvarindex(ivartype)
195           write(cunit) 1           write(cunit) 1
196           write(cunit) 1           write(cunit) 1
197    #endif
198           do k = 1, nr           do k = 1, nr
199             irectrue = (irec-1)*nr + k
200                if ( doZscalePack ) then
201                   delZnorm = (delR(1)/delR(k))**delZexp
202                else
203                   delZnorm = 1. _d 0
204                endif
205              cbuffindex = 0              cbuffindex = 0
206              do jp = 1,nPy              do jp = 1,nPy
207               do bj = jtlo,jthi               do bj = jtlo,jthi
# Line 125  c--   Only the master thread will do I/O Line 209  c--   Only the master thread will do I/O
209                 do ip = 1,nPx                 do ip = 1,nPx
210                  do bi = itlo,ithi                  do bi = itlo,ithi
211                   do i = imin,imax                   do i = imin,imax
212                    if (globmsk(i,bi,ip,j,bj,jp,k)  .ne. 0. ) then                    if (globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
213                       cbuffindex = cbuffindex + 1                       cbuffindex = cbuffindex + 1
214    cph(
215                         globfldtmp3(i,bi,ip,j,bj,jp) =
216         &                    globfld3d(i,bi,ip,j,bj,jp,k)
217    cph)
218                      IF ( doPackOld ) THEN
219  #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO  #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
220                       if (lxxadxx) then                       if (lxxadxx) then
221                          cbuff(cbuffindex) =                          cbuff(cbuffindex) = 1/delZnorm
222       &                       globfld3d(i,bi,ip,j,bj,jp,k) *       &                       * globfld3d(i,bi,ip,j,bj,jp,k)
223       &                       sqrt(weightfld(k,bi,bj))  # ifdef CTRL_PACK_PRECISE
224         &                       * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
225    # else
226         &                       * sqrt(weightfld(k,bi,bj))
227    # endif
228                       else                       else
229                            cbuff(cbuffindex) = delZnorm
230         &                       * globfld3d(i,bi,ip,j,bj,jp,k)
231    # ifdef CTRL_PACK_PRECISE
232         &                       / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
233    # else
234         &                       / sqrt(weightfld(k,bi,bj))
235    # endif
236                         endif
237  cph(  cph(
238                          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)  
239  cph)  cph)
240                       endif  #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
241  #else                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
242    #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
243                         ELSE !IF ( doPackOld ) THEN
244                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
245                         ENDIF !IF ( doPackOld ) THEN
246    #ifdef ALLOW_ADMTLM
247                         nveccount = nveccount + 1
248                         phtmpadmtlm(nveccount) = cbuff(cbuffindex)
249  #endif  #endif
250                    endif                    endif
251                   enddo                   enddo
# Line 159  cph) Line 256  cph)
256              enddo              enddo
257  c           --> check cbuffindex.  c           --> check cbuffindex.
258              if ( cbuffindex .gt. 0) then              if ( cbuffindex .gt. 0) then
259    #ifndef ALLOW_ADMTLM
260                 write(cunit) cbuffindex                 write(cunit) cbuffindex
261                 write(cunit) k                 write(cunit) k
262    cph#endif
263                 write(cunit) (cbuff(ii), ii=1,cbuffindex)                 write(cunit) (cbuff(ii), ii=1,cbuffindex)
264    #endif
265                endif
266    c
267                if ( doPackDiag ) then
268                   write(cunit2,rec=irectrue) globfldtmp2
269                   write(cunit3,rec=irectrue) globfldtmp3
270              endif              endif
271    c
272           enddo           enddo
273  c  c
274  c     -- end of irec loop --  c     -- end of irec loop --
275        enddo        enddo
276    
277          if ( doPackDiag ) then
278             close ( cunit2 )
279             close ( cunit3 )
280          endif
281    
282        _END_MASTER( mythid )        _END_MASTER( mythid )
283    
284    # else
285    c     == local variables ==
286    
287          integer bi,bj
288          integer ip,jp
289          integer i,j,k
290          integer ii
291          integer il
292          integer irec
293          integer itlo,ithi
294          integer jtlo,jthi
295    
296          integer cbuffindex
297    
298          _RL msk3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
299          real*8 msk2d_buf(sNx,sNy,nSx,nSy)
300          real*8 msk2d_buf_glo(Nx,Ny)
301    
302          _RL fld3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
303          real*8 fld2d_buf(sNx,sNy,nSx,nSy)
304          real*8 fld2d_buf_glo(Nx,Ny)
305    
306          _RL fld3dDim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
307          _RL fld3dNodim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
308    
309    #ifdef CTRL_PACK_PRECISE
310          _RL wei3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
311    #endif
312    
313          real*4 cbuff      ( snx*nsx*npx*sny*nsy*npy )
314    
315          character*(80) weightname
316          _RL delZnorm
317          character*(80) cfile2, cfile3
318          _RL dummy
319    
320          LOGICAL doPackOld
321    
322    c     == external ==
323    
324          integer  ilnblnk
325          external ilnblnk
326    
327    c     == end of interface ==
328    
329    c-- part 1: preliminary reads and definitions
330    
331          doPackOld = (.NOT.ctrlSmoothCorrel3D).AND.(.NOT.ctrlUseGen)
332    
333    #ifdef CTRL_PACK_PRECISE
334    #ifdef ALLOW_AUTODIFF
335          call active_read_xyz(weighttype, wei3d, 1,
336         &    .FALSE., .FALSE., 0 , mythid, dummy)
337    #else
338          CALL READ_REC_XYZ_RL( weighttype, wei3d, 1, 1, myThid )
339    #endif
340    #endif
341    
342    #ifdef ALLOW_AUTODIFF
343          call active_read_xyz(masktype, msk3d, 1,
344         &    .FALSE., .FALSE., 0 , mythid, dummy)
345    #else
346          CALL READ_REC_XYZ_RL( masktype, msk3d, 1, 1, myThid )
347    #endif
348    
349          if ( doPackDiag ) then
350             write(cfile2(1:80),'(80a)') ' '
351             write(cfile3(1:80),'(80a)') ' '
352             il = ilnblnk( fname )
353             if ( lxxadxx ) then
354                write(cfile2(1:80),'(2a)') fname(1:il),'.pack_ctrl_adim'
355                write(cfile3(1:80),'(2a)') fname(1:il),'.pack_ctrl_dim'
356             else
357                write(cfile2(1:80),'(2a)') fname(1:il),'.pack_grad_adim'
358                write(cfile3(1:80),'(2a)') fname(1:il),'.pack_grad_dim'
359             endif
360          endif
361    
362    c-- part 2: loop over records
363    
364          do irec = 1, ncvarrecs(ivartype)
365    
366    c-- 2.1:
367          call READ_REC_3D_RL( fname, ctrlprec,
368         &        Nr, fld3dDim, irec, 0, mythid)
369    
370    c-- 2.2: normalize field if needed
371          DO bj = myByLo(myThid), myByHi(myThid)
372           DO bi = myBxLo(myThid), myBxHi(myThid)
373            DO k=1,Nr
374             if ( doZscalePack ) then
375                delZnorm = (delR(1)/delR(k))**delZexp
376             else
377                delZnorm = 1. _d 0
378             endif
379             DO j=1,sNy
380              DO i=1,sNx
381               if (msk3d(i,j,k,bi,bj).EQ.0. _d 0) then
382                fld3dDim(i,j,k,bi,bj)=0. _d 0
383                fld3dNodim(i,j,k,bi,bj)=0. _d 0
384               else
385               IF ( .NOT.doPackOld ) THEN
386                fld3dNodim(i,j,k,bi,bj)=fld3dDim(i,j,k,bi,bj)
387               ELSE !IF ( .NOT.doPackOld ) THEN
388    # ifndef ALLOW_NONDIMENSIONAL_CONTROL_IO
389                fld3dNodim(i,j,k,bi,bj) = fld3dDim(i,j,k,bi,bj)
390    # else
391                if (lxxadxx) then
392                   fld3dNodim(i,j,k,bi,bj) =
393         &              fld3dDim(i,j,k,bi,bj) / delZnorm
394    #  ifdef CTRL_PACK_PRECISE
395         &              * sqrt(wei3d(i,j,k,bi,bj))
396    #  else
397         &              * sqrt(weightfld(k,bi,bj))
398    #  endif
399                else
400                   fld3dNodim(i,j,k,bi,bj) =
401         &              fld3dDim(i,j,k,bi,bj) * delZnorm
402    #  ifdef CTRL_PACK_PRECISE
403         &              / sqrt(wei3d(i,j,k,bi,bj))
404    #  else
405         &              / sqrt(weightfld(k,bi,bj))
406    #  endif
407                endif
408    # endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
409               ENDIF !IF ( .NOT.doPackOld ) THEN
410               endif
411              ENDDO
412             ENDDO
413            ENDDO
414           ENDDO
415          ENDDO
416    
417    c-- 2.3:
418          if ( doPackDiag ) then
419    c error: twice the same one
420          call WRITE_REC_3D_RL( cfile2, ctrlprec,
421         &        Nr, fld3dNodim, irec, 0, mythid)
422          call WRITE_REC_3D_RL( cfile3, ctrlprec,
423         &        Nr, fld3dDim, irec, 0, mythid)
424          endif
425    
426    c-- 2.4: array -> buffer -> global buffer -> global file
427    
428    #ifndef ALLOW_ADMTLM
429          _BEGIN_MASTER( mythid )
430          IF ( myProcId .eq. 0 ) THEN
431             write(cunit) ncvarindex(ivartype)
432             write(cunit) 1
433             write(cunit) 1
434          ENDIF
435          _END_MASTER( mythid )
436          _BARRIER
437    #endif
438    
439          do k = 1, nr
440    
441            CALL MDS_PASS_R8toRL( fld2d_buf, fld3dNodim,
442         &                 0, 0, 1, k, Nr, 0, 0, .FALSE., myThid )
443            CALL BAR2( myThid )
444            CALL GATHER_2D_R8( fld2d_buf_glo, fld2d_buf,
445         &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
446            CALL BAR2( myThid )
447    
448            CALL MDS_PASS_R8toRL( msk2d_buf, msk3d,
449         &                 0, 0, 1, k, Nr, 0, 0, .FALSE., myThid )
450            CALL BAR2( myThid )
451            CALL GATHER_2D_R8( msk2d_buf_glo, msk2d_buf,
452         &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
453            CALL BAR2( myThid )
454    
455            _BEGIN_MASTER( mythid )
456            cbuffindex = 0
457            IF ( myProcId .eq. 0 ) THEN
458    
459            DO j=1,Ny
460              DO i=1,Nx
461                if (msk2d_buf_glo(i,j) .ne. 0. ) then
462                   cbuffindex = cbuffindex + 1
463                   cbuff(cbuffindex) = fld2d_buf_glo(i,j)
464    #ifdef ALLOW_ADMTLM
465                   nveccount = nveccount + 1
466                   phtmpadmtlm(nveccount) = cbuff(cbuffindex)
467    #endif
468                endif
469              ENDDO
470            ENDDO
471    
472    #ifndef ALLOW_ADMTLM
473            if ( cbuffindex .gt. 0) then
474              write(cunit) cbuffindex
475              write(cunit) k
476              write(cunit) (cbuff(ii), ii=1,cbuffindex)
477            endif
478    #endif
479    
480            ENDIF
481            _END_MASTER( mythid )
482            _BARRIER
483    
484          enddo
485          enddo
486    
487    # endif /* ALLOW_PACKUNPACK_METHOD2 */
488    # endif /* EXCLUDE_CTRL_PACK */
489    
490        return        return
491        end        end
   
Legend:
Removed from v.1.1.4.1  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22