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