/[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.19 by jmc, Tue Oct 9 00:00:01 2007 UTC revision 1.26 by gforget, Tue Aug 6 20:57:00 2013 UTC
# Line 1  Line 1 
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 45  c     == routine arguments == Line 45  c     == routine arguments ==
45        integer mythid        integer mythid
46    
47  #ifndef EXCLUDE_CTRL_PACK  #ifndef EXCLUDE_CTRL_PACK
48    # ifndef ALLOW_PACKUNPACK_METHOD2
49  c     == local variables ==  c     == local variables ==
50    
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 63  c     == local variables == Line 63  c     == local variables ==
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 )        real*4 cbuff      ( snx*nsx*npx*sny*nsy*npy )
71        real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )        real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
72        real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )        real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )
73    
       character*(80) weightname  
   
74        _RL delZnorm        _RL delZnorm
75        integer reclen, irectrue        integer reclen, irectrue
76        integer cunit2, cunit3        integer cunit2, cunit3
# Line 101  c     == end of interface == Line 101  c     == end of interface ==
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                             globfldtmp2(i,bi,ip,j,bj,jp)   = 0.               globfldtmp2(i,bi,ip,j,bj,jp)   = 0. _d 0
113                             globfldtmp3(i,bi,ip,j,bj,jp)   = 0.               globfldtmp3(i,bi,ip,j,bj,jp)   = 0. _d 0
                         enddo  
                      enddo  
                   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.
# Line 126  c--   Only the master thread will do I/O Line 126  c--   Only the master thread will do I/O
126           write(cfile2(1:80),'(80a)') ' '           write(cfile2(1:80),'(80a)') ' '
127           write(cfile3(1:80),'(80a)') ' '           write(cfile3(1:80),'(80a)') ' '
128           if ( lxxadxx ) then           if ( lxxadxx ) then
129              write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')              write(cfile2(1:80),'(a,I3.3,a,I4.4,a)')
130       &           'diag_pack_nonout_ctrl_',       &           'diag_pack_nonout_ctrl_',
131       &           ivartype, '_', optimcycle, '.bin'       &           ivartype, '_', optimcycle, '.bin'
132              write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')              write(cfile3(1:80),'(a,I3.3,a,I4.4,a)')
133       &           'diag_pack_dimout_ctrl_',       &           'diag_pack_dimout_ctrl_',
134       &           ivartype, '_', optimcycle, '.bin'       &           ivartype, '_', optimcycle, '.bin'
135           else           else
136              write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')              write(cfile2(1:80),'(a,I3.3,a,I4.4,a)')
137       &           'diag_pack_nonout_grad_',       &           'diag_pack_nonout_grad_',
138       &           ivartype, '_', optimcycle, '.bin'       &           ivartype, '_', optimcycle, '.bin'
139              write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')              write(cfile3(1:80),'(a,I3.3,a,I4.4,a)')
140       &           'diag_pack_dimout_grad_',       &           'diag_pack_dimout_grad_',
141       &           ivartype, '_', optimcycle, '.bin'       &           ivartype, '_', optimcycle, '.bin'
142           endif           endif
# Line 151  c--   Only the master thread will do I/O Line 151  c--   Only the master thread will do I/O
151        endif        endif
152    
153  #ifdef CTRL_PACK_PRECISE  #ifdef CTRL_PACK_PRECISE
154        il=ilnblnk( weighttype)        if (weighttype.NE.' ') then
155        write(weightname(1:80),'(80a)') ' '         il=ilnblnk( weighttype)
156        write(weightname(1:80),'(a)') weighttype(1:il)         write(weightname(1:80),'(80a)') ' '
157           write(weightname(1:80),'(a)') weighttype(1:il)
158        call MDSREADFIELD_3D_GL(         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(
# Line 194  cph( Line 211  cph(
211                       globfldtmp3(i,bi,ip,j,bj,jp) =                       globfldtmp3(i,bi,ip,j,bj,jp) =
212       &                    globfld3d(i,bi,ip,j,bj,jp,k)       &                    globfld3d(i,bi,ip,j,bj,jp,k)
213  cph)  cph)
214  #ifndef ALLOW_SMOOTH_CORREL3D                    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) = 1/delZnorm                          cbuff(cbuffindex) = 1/delZnorm
# Line 219  cph) Line 236  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  /* ALLOW_SMOOTH_CORREL3D */                       ELSE !IF ( .NOT.ctrlSmoothCorrel3D ) THEN
240                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)                       cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
241  #endif /* ALLOW_SMOOTH_CORREL3D */                       ENDIF !IF ( .NOT.ctrlSmoothCorrel3D ) THEN
242  #ifdef ALLOW_ADMTLM  #ifdef ALLOW_ADMTLM
243                       nveccount = nveccount + 1                       nveccount = nveccount + 1
244                       phtmpadmtlm(nveccount) = cbuff(cbuffindex)                       phtmpadmtlm(nveccount) = cbuff(cbuffindex)
# Line 260  c     -- end of irec loop -- Line 277  c     -- end of irec loop --
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          call active_read_xyz(weighttype, wei3d, 1,
327         &    .FALSE., .FALSE., 0 , mythid, dummy)
328    #endif
329    
330          call active_read_xyz(masktype, msk3d, 1,
331         &    .FALSE., .FALSE., 0 , mythid, dummy)
332    
333          if ( doPackDiag ) then
334             write(cfile2(1:80),'(80a)') ' '
335             write(cfile3(1:80),'(80a)') ' '
336             il = ilnblnk( fname )
337             if ( lxxadxx ) then
338                write(cfile2(1:80),'(2a)') fname(1:il),'.pack_ctrl_adim'
339                write(cfile3(1:80),'(2a)') fname(1:il),'.pack_ctrl_dim'
340             else
341                write(cfile2(1:80),'(2a)') fname(1:il),'.pack_grad_adim'
342                write(cfile3(1:80),'(2a)') fname(1:il),'.pack_grad_dim'
343             endif
344          endif
345    
346    c-- part 2: loop over records
347    
348          do irec = 1, ncvarrecs(ivartype)
349    
350    c-- 2.1:
351          call READ_REC_3D_RL( fname, ctrlprec,
352         &        Nr, fld3dDim, irec, 0, mythid)
353    
354    c-- 2.2: normalize field if needed
355          DO bj = myByLo(myThid), myByHi(myThid)
356           DO bi = myBxLo(myThid), myBxHi(myThid)
357            DO k=1,Nr
358             if ( doZscalePack ) then
359                delZnorm = (delR(1)/delR(k))**delZexp
360             else
361                delZnorm = 1. _d 0
362             endif
363             DO j=1,sNy
364              DO i=1,sNx
365               if (msk3d(i,j,k,bi,bj).EQ.0. _d 0) then
366                fld3dDim(i,j,k,bi,bj)=0. _d 0
367                fld3dNodim(i,j,k,bi,bj)=0. _d 0
368               else
369               IF ( ctrlSmoothCorrel3D ) THEN
370                fld3dNodim(i,j,k,bi,bj)=fld3dDim(i,j,k,bi,bj)
371               ELSE !IF ( ctrlSmoothCorrel3D ) THEN
372    # ifndef ALLOW_NONDIMENSIONAL_CONTROL_IO
373                fld3dNodim(i,j,k,bi,bj) = fld3dDim(i,j,k,bi,bj)
374    # else
375                if (lxxadxx) then
376                   fld3dNodim(i,j,k,bi,bj) =
377         &              fld3dDim(i,j,k,bi,bj) / delZnorm
378    #  ifdef CTRL_PACK_PRECISE
379         &              * sqrt(wei3d(i,j,k,bi,bj))
380    #  else
381         &              * sqrt(weightfld(k,bi,bj))
382    #  endif
383                else
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                endif
392    # endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
393               ENDIF !IF ( ctrlSmoothCorrel3D ) THEN
394               endif
395              ENDDO
396             ENDDO
397            ENDDO
398           ENDDO
399          ENDDO
400    
401    c-- 2.3:
402          if ( doPackDiag ) then
403    c error: twice the same one
404          call WRITE_REC_3D_RL( cfile2, ctrlprec,
405         &        Nr, fld3dNodim, irec, 0, mythid)
406          call WRITE_REC_3D_RL( cfile3, ctrlprec,
407         &        Nr, fld3dDim, irec, 0, mythid)
408          endif
409    
410    c-- 2.4: array -> buffer -> global buffer -> global file
411    
412    #ifndef ALLOW_ADMTLM
413          _BEGIN_MASTER( mythid )
414          IF ( myProcId .eq. 0 ) THEN
415             write(cunit) ncvarindex(ivartype)
416             write(cunit) 1
417             write(cunit) 1
418          ENDIF
419          _END_MASTER( mythid )
420          _BARRIER
421    #endif
422    
423          do k = 1, nr
424    
425            CALL MDS_PASS_R8toRL( fld2d_buf, fld3dNodim,
426         &                 0, 0, 1, k, Nr, 0, 0, .FALSE., myThid )
427            CALL BAR2( myThid )
428            CALL GATHER_2D_R8( fld2d_buf_glo, fld2d_buf,
429         &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
430            CALL BAR2( myThid )
431    
432            CALL MDS_PASS_R8toRL( msk2d_buf, msk3d,
433         &                 0, 0, 1, k, Nr, 0, 0, .FALSE., myThid )
434            CALL BAR2( myThid )
435            CALL GATHER_2D_R8( msk2d_buf_glo, msk2d_buf,
436         &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
437            CALL BAR2( myThid )
438    
439            _BEGIN_MASTER( mythid )
440            cbuffindex = 0
441            IF ( myProcId .eq. 0 ) THEN
442    
443            DO j=1,Ny
444              DO i=1,Nx
445                if (msk2d_buf_glo(i,j) .ne. 0. ) then
446                   cbuffindex = cbuffindex + 1
447                   cbuff(cbuffindex) = fld2d_buf_glo(i,j)
448    #ifdef ALLOW_ADMTLM
449                   nveccount = nveccount + 1
450                   phtmpadmtlm(nveccount) = cbuff(cbuffindex)
451    #endif
452                endif
453              ENDDO
454            ENDDO
455    
456    #ifndef ALLOW_ADMTLM
457            if ( cbuffindex .gt. 0) then
458              write(cunit) cbuffindex
459              write(cunit) k
460              write(cunit) (cbuff(ii), ii=1,cbuffindex)
461            endif
462  #endif  #endif
463    
464            ENDIF
465            _END_MASTER( mythid )
466            _BARRIER
467    
468          enddo
469          enddo
470    
471    # endif /* ALLOW_PACKUNPACK_METHOD2 */
472    # endif /* EXCLUDE_CTRL_PACK */
473    
474        return        return
475        end        end
   
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22