pkg/ctrl/ctrl_set_unpack_xyz.F

C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_set_unpack_xyz.F,v 1.21 2010/09/25 23:04:07 gforget Exp $
C $Name:  $

#include "CTRL_CPPOPTIONS.h"

      subroutine ctrl_set_unpack_xyz(
     &     lxxadxx, cunit, ivartype, fname, masktype, weighttype,
     &     weightfld, nwetglobal, mythid)

c     ==================================================================
c     SUBROUTINE ctrl_set_unpack_xyz
c     ==================================================================
c
c     o Unpack the control vector such that land points are filled in.
c
c     o Use a more precise nondimensionalization that depends on (x,y)
c       Added weighttype to the argument list so that I can geographically
c       vary the nondimensionalization.
c       gebbie@mit.edu, 18-Mar-2003
c
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"

#include "ctrl.h"
#include "optim.h"

c     == routine arguments ==

      logical lxxadxx
      integer cunit
      integer ivartype
      character*( 80)   fname
      character*(  9) masktype
      character*( 80) weighttype
      _RL     weightfld( nr,nsx,nsy )
      integer nwetglobal(nr)
      integer mythid

#ifndef EXCLUDE_CTRL_PACK
# ifndef ALLOW_PACKUNPACK_METHOD2
c     == local variables ==

      integer bi,bj
      integer ip,jp
      integer i,j,k
      integer ii
      integer il
      integer irec
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax

      integer cbuffindex

      _RL     globmsk  ( snx,nsx,npx,sny,nsy,npy,nr )
      _RL     globfld3d( snx,nsx,npx,sny,nsy,npy,nr )
#ifdef CTRL_UNPACK_PRECISE
      _RL   weightfld3d( snx,nsx,npx,sny,nsy,npy,nr )
#endif
      real*4 cbuff      ( snx*nsx*npx*sny*nsy*npy )
      real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
      real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )

      character*(128)   cfile
      character*(80) weightname

      _RL delZnorm
      integer reclen, irectrue
      integer cunit2, cunit3
      character*(80) cfile2, cfile3

c     == external ==

      integer  ilnblnk
      external ilnblnk

cc     == end of interface ==

      jtlo = 1
      jthi = nsy
      itlo = 1
      ithi = nsx
      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

#ifdef CTRL_DELZNORM
      delZnorm = 0.
      do k = 1, Nr
         delZnorm = delZnorm + delR(k)/FLOAT(Nr)
      enddo
#endif

c     Initialise temporary file
      do k = 1,nr
         do jp = 1,nPy
            do bj = jtlo,jthi
               do j = jmin,jmax
                  do ip = 1,nPx
                     do bi = itlo,ithi
                        do i = imin,imax
                           globfld3d  (i,bi,ip,j,bj,jp,k) = 0. _d 0
                           globmsk    (i,bi,ip,j,bj,jp,k) = 0. _d 0
                           globfldtmp2(i,bi,ip,j,bj,jp)   = 0.
                           globfldtmp3(i,bi,ip,j,bj,jp)   = 0.
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo

c--   Only the master thread will do I/O.
      _BEGIN_MASTER( mythid )

#ifdef CTRL_DELZNORM
      do k = 1, nr
         print *, 'ph-delznorm ', k, delZnorm, delR(k)
         print *, 'ph-weight   ', weightfld(k,1,1)
      enddo
#endif

      if ( doPackDiag ) then
         write(cfile2(1:80),'(80a)') ' '
         write(cfile3(1:80),'(80a)') ' '
         if ( lxxadxx ) then
            write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
     &           'diag_unpack_nondim_ctrl_',
     &           ivartype, '_', optimcycle, '.bin'
            write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
     &           'diag_unpack_dimens_ctrl_',
     &           ivartype, '_', optimcycle, '.bin'
         else
            write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
     &           'diag_unpack_nondim_grad_',
     &           ivartype, '_', optimcycle, '.bin'
            write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
     &           'diag_unpack_dimens_grad_',
     &           ivartype, '_', optimcycle, '.bin'
         endif

         reclen = FLOAT(snx*nsx*npx*sny*nsy*npy*4)
         call mdsfindunit( cunit2, mythid )
         open( cunit2, file=cfile2, status='unknown',
     &        access='direct', recl=reclen )
         call mdsfindunit( cunit3, mythid )
         open( cunit3, file=cfile3, status='unknown',
     &        access='direct', recl=reclen )
      endif

#ifdef CTRL_UNPACK_PRECISE
      il=ilnblnk( weighttype)
      write(weightname(1:80),'(80a)') ' '
      write(weightname(1:80),'(a)') weighttype(1:il)

      call MDSREADFIELD_3D_GL(
     &     weightname, ctrlprec, 'RL',
     &     Nr, weightfld3d, 1, mythid)
#endif

      call MDSREADFIELD_3D_GL(
     &     masktype, ctrlprec, 'RL',
     &     Nr, globmsk, 1, mythid)

      do irec = 1, ncvarrecs(ivartype)
#ifndef ALLOW_ADMTLM
         read(cunit) filencvarindex(ivartype)
         if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
     &        then
            print *, 'ctrl_set_unpack_xyz:WARNING: wrong ncvarindex ',
     &           filencvarindex(ivartype), ncvarindex(ivartype)
            STOP 'in S/R ctrl_unpack'
         endif
         read(cunit) filej
         read(cunit) filei
#endif /* ALLOW_ADMTLM */
         do k = 1, Nr
         irectrue = (irec-1)*nr + k
            if ( doZscaleUnpack ) then
               delZnorm = (delR(1)/delR(k))**delZexp
            else
               delZnorm = 1. _d 0
            endif
            cbuffindex = nwetglobal(k)
            if ( cbuffindex .gt. 0 ) then
#ifndef ALLOW_ADMTLM
               read(cunit) filencbuffindex
               if (filencbuffindex .NE. cbuffindex) then
                  print *, 'WARNING: wrong cbuffindex ',
     &                 filencbuffindex, cbuffindex
                  STOP 'in S/R ctrl_unpack'
               endif
               read(cunit) filek
               if (filek .NE. k) then
                  print *, 'WARNING: wrong k ',
     &                 filek, k
                  STOP 'in S/R ctrl_unpack'
               endif
cph#endif /* ALLOW_ADMTLM */
               read(cunit) (cbuff(ii), ii=1,cbuffindex)
#endif /* ALLOW_ADMTLM */
            endif
c
            cbuffindex = 0
            do jp = 1,nPy
             do bj = jtlo,jthi
              do j = jmin,jmax
               do ip = 1,nPx
                do bi = itlo,ithi
                 do i = imin,imax
                  if ( globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
                     cbuffindex = cbuffindex + 1
                     globfld3d(i,bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
cph(
                     globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
cph)
#ifdef ALLOW_ADMTLM
                     nveccount = nveccount + 1
                     globfld3d(i,bi,ip,j,bj,jp,k) =
     &                 phtmpadmtlm(nveccount)
cph(
                     globfldtmp2(i,bi,ip,j,bj,jp) =
     &                 phtmpadmtlm(nveccount)
cph)
#endif
#ifndef ALLOW_SMOOTH_CORREL3D
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                     if ( lxxadxx ) then
                        globfld3d(i,bi,ip,j,bj,jp,k) = delZnorm
     &                       * globfld3d(i,bi,ip,j,bj,jp,k)
# ifdef CTRL_UNPACK_PRECISE
     &                       / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
# else
     &                       / sqrt(weightfld(k,bi,bj))
# endif
                     else
                        globfld3d(i,bi,ip,j,bj,jp,k) = 1/delZnorm
     &                       * globfld3d(i,bi,ip,j,bj,jp,k)
# ifdef CTRL_UNPACK_PRECISE
     &                       * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
# else
     &                       * sqrt(weightfld(k,bi,bj))
# endif
                     endif
#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
#endif /* ALLOW_SMOOTH_CORREL3D */
                  else
                     globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
                  endif
cph(
                  globfldtmp3(i,bi,ip,j,bj,jp) =
     &                 globfld3d(i,bi,ip,j,bj,jp,k)
cph)
                 enddo
                enddo
               enddo
              enddo
             enddo
            enddo
c
            if ( doPackDiag ) then
               write(cunit2,rec=irectrue) globfldtmp2
               write(cunit3,rec=irectrue) globfldtmp3
            endif
c
         enddo

         call MDSWRITEFIELD_3D_GL( fname, ctrlprec, 'RL',
     &                             Nr, globfld3d,
     &                             irec,  optimcycle, mythid)

      enddo

      if ( doPackDiag ) then
         close ( cunit2 )
         close ( cunit3 )
      endif

      _END_MASTER( mythid )

# else
c     == local variables ==

      integer bi,bj
      integer ip,jp
      integer i,j,k
      integer ii
      integer il
      integer irec
      integer itlo,ithi
      integer jtlo,jthi

      integer cbuffindex

      _RL msk3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      real*8 msk2d_buf(sNx,sNy,nSx,nSy)
      real*8 msk2d_buf_glo(Nx,Ny)

      _RL fld3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      real*8 fld2d_buf(sNx,sNy,nSx,nSy)
      real*8 fld2d_buf_glo(Nx,Ny)

      _RL fld3dDim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL fld3dNodim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)

#ifdef CTRL_PACK_PRECISE
      _RL wei3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
#endif

      real*4 cbuff      ( snx*nsx*npx*sny*nsy*npy )

      character*(80) weightname
      _RL delZnorm
      character*(80) cfile2, cfile3

c     == external ==

      integer  ilnblnk
      external ilnblnk

c     == end of interface ==

c-- part 1: preliminary reads and definitions

#ifdef CTRL_PACK_PRECISE
      call active_read_xyz(weighttype, wei3d, 1,
     &    .FALSE., .FALSE., 0 , mythid, 1)
#endif

      call active_read_xyz(masktype, msk3d, 1, 
     &    .FALSE., .FALSE., 0 , mythid, 1)

      if ( doPackDiag ) then
         write(cfile2(1:80),'(80a)') ' '
         write(cfile3(1:80),'(80a)') ' '
         il = ilnblnk( fname )
         if ( lxxadxx ) then
            write(cfile2(1:80),'(2a)') fname(1:il),'.unpack_ctrl_adim'
            write(cfile3(1:80),'(2a)') fname(1:il),'.unpack_ctrl_dim'
         else
            write(cfile2(1:80),'(2a)') fname(1:il),'.unpack_grad_adim'
            write(cfile3(1:80),'(2a)') fname(1:il),'.unpack_grad_dim'
         endif
      endif

c-- part 2: loop over records

      do irec = 1, ncvarrecs(ivartype)

c-- 2.1: array <- buffer <- global buffer <- global file

#ifndef ALLOW_ADMTLM
      _BEGIN_MASTER( mythid )
        IF ( myProcId .eq. 0 ) THEN
         read(cunit) filencvarindex(ivartype)
         if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
     &        then
            print *, 'ctrl_set_unpack_xyz:WARNING: wrong ncvarindex ',
     &           filencvarindex(ivartype), ncvarindex(ivartype)
            STOP 'in S/R ctrl_unpack'
         endif
         read(cunit) filej
         read(cunit) filei
        ENDIF
      _END_MASTER( mythid )
      _BARRIER
#endif /* ALLOW_ADMTLM */

      do k = 1, nr

        CALL MDS_PASS_R8toRL( msk2d_buf, msk3d,
     &                       1, k, Nr, 0, 0, .FALSE., myThid )
        CALL BAR2( myThid )
        CALL GATHER_2D_R8( msk2d_buf_glo, msk2d_buf,
     &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
        CALL BAR2( myThid )

        _BEGIN_MASTER( mythid )
         cbuffindex = nwetglobal(k)
        IF ( myProcId .eq. 0 ) THEN

#ifndef ALLOW_ADMTLM
        if ( cbuffindex .gt. 0) then
          read(cunit) filencbuffindex
          read(cunit) filek
          if (filencbuffindex .NE. cbuffindex) then
            print *, 'WARNING: wrong cbuffindex ',
     &                 filencbuffindex, cbuffindex
            STOP 'in S/R ctrl_unpack'
          endif
          if (filek .NE. k) then
            print *, 'WARNING: wrong k ',
     &                 filek, k
            STOP 'in S/R ctrl_unpack'
          endif
          read(cunit) (cbuff(ii), ii=1,cbuffindex)
        endif
#endif

        cbuffindex = 0
        DO j=1,Ny
          DO i=1,Nx
            if (msk2d_buf_glo(i,j) .ne. 0. ) then
               cbuffindex = cbuffindex + 1
               fld2d_buf_glo(i,j) = cbuff(cbuffindex)
            endif
          ENDDO
        ENDDO

        ENDIF
        _END_MASTER( mythid )
        _BARRIER

        CALL BAR2( myThid )
        CALL SCATTER_2D_R8( fld2d_buf_glo, fld2d_buf,
     &                       Nx,Ny,.FALSE.,.TRUE.,myThid)
        CALL BAR2( myThid )
        CALL MDS_PASS_R8toRL( fld2d_buf, fld3dNodim,
     &                       1, k, nr, 0, 0, .TRUE., myThid )

      enddo !do k = 1, nr

c-- 2.2: normalize field if needed
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Nr
         if ( doZscalePack ) then
            delZnorm = (delR(1)/delR(k))**delZexp
         else
            delZnorm = 1. _d 0
         endif
         DO j=1,sNy
          DO i=1,sNx
           if (msk3d(i,j,k,bi,bj).EQ.0. _d 0) then
            fld3dDim(i,j,k,bi,bj)=0. _d 0
            fld3dNodim(i,j,k,bi,bj)=0. _d 0
           else
#ifdef ALLOW_ADMTLM
            nveccount = nveccount + 1
            fld3dNodim(i,j,k,bi,bj)=phtmpadmtlm(nveccount)
#endif           
#ifdef ALLOW_SMOOTH_CORREL3D
            fld3dDim(i,j,k,bi,bj)=fld3dNodim(i,j,k,bi,bj)
#else
# ifndef ALLOW_NONDIMENSIONAL_CONTROL_IO
            fld3dDim(i,j,k,bi,bj) = fld3dNodim(i,j,k,bi,bj)
# else
            if (lxxadxx) then
               fld3dDim(i,j,k,bi,bj) = 
     &              fld3dNodim(i,j,k,bi,bj) * delZnorm 
#  ifdef CTRL_PACK_PRECISE
     &              / sqrt(wei3d(i,j,k,bi,bj))
#  else
     &              / sqrt(weightfld(k,bi,bj))
#  endif
            else
               fld3dDim(i,j,k,bi,bj) =
     &              fld3dNodim(i,j,k,bi,bj) / delZnorm
#  ifdef CTRL_PACK_PRECISE
     &              * sqrt(wei3d(i,j,k,bi,bj))
#  else
     &              * sqrt(weightfld(k,bi,bj))
#  endif
            endif
# endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
#endif /* ALLOW_SMOOTH_CORREL3D */
           endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c-- 2.3:
      if ( doPackDiag ) then
c error: twice the same one
      call WRITE_REC_3D_RL( cfile2, ctrlprec,
     &        Nr, fld3dNodim, irec, 0, mythid)
      call WRITE_REC_3D_RL( cfile3, ctrlprec,
     &        Nr, fld3dDim, irec, 0, mythid)
      endif

c-- 2.4: 
      call WRITE_REC_3D_RL( fname, ctrlprec, 
     &        Nr, fld3dDim, irec, 0, mythid)


      enddo !do irec = 1, ncvarrecs(ivartype)

# endif /* ALLOW_PACKUNPACK_METHOD2 */ 
# endif /* EXCLUDE_CTRL_PACK */

      return
      end

1	gforget	1.22	C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_set_unpack_xyz.F,v 1.21 2010/09/25 23:04:07 gforget Exp $
2	jmc	1.19	C $Name: $
3	heimbach	1.2
4			#include "CTRL_CPPOPTIONS.h"
5
6			subroutine ctrl_set_unpack_xyz(
7	heimbach	1.9	& lxxadxx, cunit, ivartype, fname, masktype, weighttype,
8	heimbach	1.2	& weightfld, nwetglobal, mythid)
9
10			c ==================================================================
11			c SUBROUTINE ctrl_set_unpack_xyz
12			c ==================================================================
13			c
14	heimbach	1.3	c o Unpack the control vector such that land points are filled in.
15			c
16			c o Use a more precise nondimensionalization that depends on (x,y)
17			c Added weighttype to the argument list so that I can geographically
18			c vary the nondimensionalization.
19			c gebbie@mit.edu, 18-Mar-2003
20	heimbach	1.2	c
21			c ==================================================================
22
23			implicit none
24
25			c == global variables ==
26
27			#include "EEPARAMS.h"
28			#include "SIZE.h"
29			#include "PARAMS.h"
30			#include "GRID.h"
31
32			#include "ctrl.h"
33			#include "optim.h"
34
35			c == routine arguments ==
36
37	heimbach	1.9	logical lxxadxx
38	heimbach	1.2	integer cunit
39			integer ivartype
40			character*( 80) fname
41	heimbach	1.11	character*( 9) masktype
42	heimbach	1.3	character*( 80) weighttype
43	heimbach	1.2	_RL weightfld( nr,nsx,nsy )
44			integer nwetglobal(nr)
45			integer mythid
46
47	heimbach	1.17	#ifndef EXCLUDE_CTRL_PACK
48	gforget	1.22	# ifndef ALLOW_PACKUNPACK_METHOD2
49	heimbach	1.2	c == local variables ==
50
51			integer bi,bj
52			integer ip,jp
53			integer i,j,k
54			integer ii
55			integer il
56			integer irec
57			integer itlo,ithi
58			integer jtlo,jthi
59			integer jmin,jmax
60			integer imin,imax
61
62			integer cbuffindex
63
64			_RL globmsk ( snx,nsx,npx,sny,nsy,npy,nr )
65			_RL globfld3d( snx,nsx,npx,sny,nsy,npy,nr )
66	heimbach	1.3	#ifdef CTRL_UNPACK_PRECISE
67			_RL weightfld3d( snx,nsx,npx,sny,nsy,npy,nr )
68			#endif
69	heimbach	1.9	real4 cbuff ( snxnsxnpxsnynsynpy )
70			real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
71			real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )
72	heimbach	1.2
73			character*(128) cfile
74	heimbach	1.3	character*(80) weightname
75	heimbach	1.2
76	heimbach	1.9	_RL delZnorm
77			integer reclen, irectrue
78			integer cunit2, cunit3
79			character*(80) cfile2, cfile3
80
81	heimbach	1.2	c == external ==
82
83			integer ilnblnk
84			external ilnblnk
85
86			cc == end of interface ==
87
88			jtlo = 1
89			jthi = nsy
90			itlo = 1
91			ithi = nsx
92			jmin = 1
93			jmax = sny
94			imin = 1
95			imax = snx
96
97	heimbach	1.9	#ifdef CTRL_DELZNORM
98			delZnorm = 0.
99			do k = 1, Nr
100			delZnorm = delZnorm + delR(k)/FLOAT(Nr)
101			enddo
102			#endif
103
104	heimbach	1.2	c Initialise temporary file
105			do k = 1,nr
106			do jp = 1,nPy
107			do bj = jtlo,jthi
108			do j = jmin,jmax
109			do ip = 1,nPx
110			do bi = itlo,ithi
111			do i = imin,imax
112	heimbach	1.9	globfld3d (i,bi,ip,j,bj,jp,k) = 0. _d 0
113			globmsk (i,bi,ip,j,bj,jp,k) = 0. _d 0
114			globfldtmp2(i,bi,ip,j,bj,jp) = 0.
115			globfldtmp3(i,bi,ip,j,bj,jp) = 0.
116	heimbach	1.2	enddo
117			enddo
118			enddo
119			enddo
120			enddo
121			enddo
122			enddo
123
124			c-- Only the master thread will do I/O.
125			_BEGIN_MASTER( mythid )
126
127	heimbach	1.9	#ifdef CTRL_DELZNORM
128			do k = 1, nr
129			print *, 'ph-delznorm ', k, delZnorm, delR(k)
130			print *, 'ph-weight ', weightfld(k,1,1)
131			enddo
132			#endif
133
134			if ( doPackDiag ) then
135			write(cfile2(1:80),'(80a)') ' '
136			write(cfile3(1:80),'(80a)') ' '
137			if ( lxxadxx ) then
138			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
139	jmc	1.19	& 'diag_unpack_nondim_ctrl_',
140	heimbach	1.9	& ivartype, '_', optimcycle, '.bin'
141			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
142	jmc	1.19	& 'diag_unpack_dimens_ctrl_',
143	heimbach	1.9	& ivartype, '_', optimcycle, '.bin'
144			else
145			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
146	jmc	1.19	& 'diag_unpack_nondim_grad_',
147	heimbach	1.9	& ivartype, '_', optimcycle, '.bin'
148			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
149	jmc	1.19	& 'diag_unpack_dimens_grad_',
150	heimbach	1.9	& ivartype, '_', optimcycle, '.bin'
151			endif
152
153			reclen = FLOAT(snxnsxnpxsnynsynpy4)
154			call mdsfindunit( cunit2, mythid )
155	heimbach	1.10	open( cunit2, file=cfile2, status='unknown',
156	heimbach	1.9	& access='direct', recl=reclen )
157			call mdsfindunit( cunit3, mythid )
158	heimbach	1.10	open( cunit3, file=cfile3, status='unknown',
159	heimbach	1.9	& access='direct', recl=reclen )
160			endif
161
162	heimbach	1.3	#ifdef CTRL_UNPACK_PRECISE
163			il=ilnblnk( weighttype)
164			write(weightname(1:80),'(80a)') ' '
165			write(weightname(1:80),'(a)') weighttype(1:il)
166
167			call MDSREADFIELD_3D_GL(
168			& weightname, ctrlprec, 'RL',
169			& Nr, weightfld3d, 1, mythid)
170			#endif
171
172	jmc	1.19	call MDSREADFIELD_3D_GL(
173	heimbach	1.2	& masktype, ctrlprec, 'RL',
174			& Nr, globmsk, 1, mythid)
175
176			do irec = 1, ncvarrecs(ivartype)
177	heimbach	1.13	#ifndef ALLOW_ADMTLM
178	heimbach	1.2	read(cunit) filencvarindex(ivartype)
179			if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
180			& then
181			print *, 'ctrl_set_unpack_xyz:WARNING: wrong ncvarindex ',
182			& filencvarindex(ivartype), ncvarindex(ivartype)
183			STOP 'in S/R ctrl_unpack'
184			endif
185			read(cunit) filej
186			read(cunit) filei
187	heimbach	1.13	#endif /* ALLOW_ADMTLM */
188	heimbach	1.2	do k = 1, Nr
189	heimbach	1.9	irectrue = (irec-1)*nr + k
190	heimbach	1.10	if ( doZscaleUnpack ) then
191	heimbach	1.16	delZnorm = (delR(1)/delR(k))**delZexp
192	heimbach	1.10	else
193			delZnorm = 1. _d 0
194			endif
195	heimbach	1.2	cbuffindex = nwetglobal(k)
196			if ( cbuffindex .gt. 0 ) then
197	heimbach	1.13	#ifndef ALLOW_ADMTLM
198	heimbach	1.2	read(cunit) filencbuffindex
199			if (filencbuffindex .NE. cbuffindex) then
200			print *, 'WARNING: wrong cbuffindex ',
201			& filencbuffindex, cbuffindex
202			STOP 'in S/R ctrl_unpack'
203			endif
204			read(cunit) filek
205			if (filek .NE. k) then
206			print *, 'WARNING: wrong k ',
207			& filek, k
208			STOP 'in S/R ctrl_unpack'
209			endif
210	heimbach	1.15	cph#endif /* ALLOW_ADMTLM */
211			read(cunit) (cbuff(ii), ii=1,cbuffindex)
212	heimbach	1.13	#endif /* ALLOW_ADMTLM */
213	heimbach	1.2	endif
214	heimbach	1.13	c
215	heimbach	1.2	cbuffindex = 0
216			do jp = 1,nPy
217			do bj = jtlo,jthi
218			do j = jmin,jmax
219			do ip = 1,nPx
220			do bi = itlo,ithi
221			do i = imin,imax
222			if ( globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
223			cbuffindex = cbuffindex + 1
224			globfld3d(i,bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
225	heimbach	1.9	cph(
226			globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
227			cph)
228	heimbach	1.13	#ifdef ALLOW_ADMTLM
229			nveccount = nveccount + 1
230	jmc	1.19	globfld3d(i,bi,ip,j,bj,jp,k) =
231	heimbach	1.14	& phtmpadmtlm(nveccount)
232			cph(
233	jmc	1.19	globfldtmp2(i,bi,ip,j,bj,jp) =
234	heimbach	1.14	& phtmpadmtlm(nveccount)
235			cph)
236	heimbach	1.13	#endif
237	gforget	1.18	#ifndef ALLOW_SMOOTH_CORREL3D
238	heimbach	1.2	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
239	heimbach	1.9	if ( lxxadxx ) then
240	heimbach	1.10	globfld3d(i,bi,ip,j,bj,jp,k) = delZnorm
241			& * globfld3d(i,bi,ip,j,bj,jp,k)
242	heimbach	1.9	# ifdef CTRL_UNPACK_PRECISE
243	heimbach	1.10	& / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
244	heimbach	1.9	# else
245	heimbach	1.10	& / sqrt(weightfld(k,bi,bj))
246	heimbach	1.9	# endif
247			else
248	gforget	1.21	globfld3d(i,bi,ip,j,bj,jp,k) = 1/delZnorm
249	heimbach	1.10	& * globfld3d(i,bi,ip,j,bj,jp,k)
250	heimbach	1.3	# ifdef CTRL_UNPACK_PRECISE
251	heimbach	1.9	& * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
252	heimbach	1.3	# else
253	heimbach	1.9	& * sqrt(weightfld(k,bi,bj))
254	heimbach	1.3	# endif
255	heimbach	1.9	endif
256	heimbach	1.3	#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
257	jmc	1.20	#endif /* ALLOW_SMOOTH_CORREL3D */
258	heimbach	1.2	else
259			globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
260			endif
261	heimbach	1.9	cph(
262			globfldtmp3(i,bi,ip,j,bj,jp) =
263			& globfld3d(i,bi,ip,j,bj,jp,k)
264			cph)
265	heimbach	1.2	enddo
266			enddo
267			enddo
268			enddo
269			enddo
270			enddo
271			c
272	heimbach	1.9	if ( doPackDiag ) then
273			write(cunit2,rec=irectrue) globfldtmp2
274			write(cunit3,rec=irectrue) globfldtmp3
275			endif
276			c
277	heimbach	1.2	enddo
278	jmc	1.19
279	heimbach	1.2	call MDSWRITEFIELD_3D_GL( fname, ctrlprec, 'RL',
280			& Nr, globfld3d,
281			& irec, optimcycle, mythid)
282
283			enddo
284
285	heimbach	1.9	if ( doPackDiag ) then
286			close ( cunit2 )
287			close ( cunit3 )
288			endif
289
290	heimbach	1.2	_END_MASTER( mythid )
291
292	gforget	1.22	# else
293			c == local variables ==
294
295			integer bi,bj
296			integer ip,jp
297			integer i,j,k
298			integer ii
299			integer il
300			integer irec
301			integer itlo,ithi
302			integer jtlo,jthi
303
304			integer cbuffindex
305
306			_RL msk3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
307			real*8 msk2d_buf(sNx,sNy,nSx,nSy)
308			real*8 msk2d_buf_glo(Nx,Ny)
309
310			_RL fld3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
311			real*8 fld2d_buf(sNx,sNy,nSx,nSy)
312			real*8 fld2d_buf_glo(Nx,Ny)
313
314			_RL fld3dDim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
315			_RL fld3dNodim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
316
317			#ifdef CTRL_PACK_PRECISE
318			_RL wei3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
319	heimbach	1.17	#endif
320
321	gforget	1.22	real4 cbuff ( snxnsxnpxsnynsynpy )
322
323			character*(80) weightname
324			_RL delZnorm
325			character*(80) cfile2, cfile3
326
327			c == external ==
328
329			integer ilnblnk
330			external ilnblnk
331
332			c == end of interface ==
333
334			c-- part 1: preliminary reads and definitions
335
336			#ifdef CTRL_PACK_PRECISE
337			call active_read_xyz(weighttype, wei3d, 1,
338			& .FALSE., .FALSE., 0 , mythid, 1)
339			#endif
340
341			call active_read_xyz(masktype, msk3d, 1,
342			& .FALSE., .FALSE., 0 , mythid, 1)
343
344			if ( doPackDiag ) then
345			write(cfile2(1:80),'(80a)') ' '
346			write(cfile3(1:80),'(80a)') ' '
347			il = ilnblnk( fname )
348			if ( lxxadxx ) then
349			write(cfile2(1:80),'(2a)') fname(1:il),'.unpack_ctrl_adim'
350			write(cfile3(1:80),'(2a)') fname(1:il),'.unpack_ctrl_dim'
351			else
352			write(cfile2(1:80),'(2a)') fname(1:il),'.unpack_grad_adim'
353			write(cfile3(1:80),'(2a)') fname(1:il),'.unpack_grad_dim'
354			endif
355			endif
356
357			c-- part 2: loop over records
358
359			do irec = 1, ncvarrecs(ivartype)
360
361			c-- 2.1: array <- buffer <- global buffer <- global file
362
363			#ifndef ALLOW_ADMTLM
364			_BEGIN_MASTER( mythid )
365			IF ( myProcId .eq. 0 ) THEN
366			read(cunit) filencvarindex(ivartype)
367			if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
368			& then
369			print *, 'ctrl_set_unpack_xyz:WARNING: wrong ncvarindex ',
370			& filencvarindex(ivartype), ncvarindex(ivartype)
371			STOP 'in S/R ctrl_unpack'
372			endif
373			read(cunit) filej
374			read(cunit) filei
375			ENDIF
376			_END_MASTER( mythid )
377			_BARRIER
378			#endif /* ALLOW_ADMTLM */
379
380			do k = 1, nr
381
382			CALL MDS_PASS_R8toRL( msk2d_buf, msk3d,
383			& 1, k, Nr, 0, 0, .FALSE., myThid )
384			CALL BAR2( myThid )
385			CALL GATHER_2D_R8( msk2d_buf_glo, msk2d_buf,
386			& Nx,Ny,.FALSE.,.TRUE.,myThid)
387			CALL BAR2( myThid )
388
389			_BEGIN_MASTER( mythid )
390			cbuffindex = nwetglobal(k)
391			IF ( myProcId .eq. 0 ) THEN
392
393			#ifndef ALLOW_ADMTLM
394			if ( cbuffindex .gt. 0) then
395			read(cunit) filencbuffindex
396			read(cunit) filek
397			if (filencbuffindex .NE. cbuffindex) then
398			print *, 'WARNING: wrong cbuffindex ',
399			& filencbuffindex, cbuffindex
400			STOP 'in S/R ctrl_unpack'
401			endif
402			if (filek .NE. k) then
403			print *, 'WARNING: wrong k ',
404			& filek, k
405			STOP 'in S/R ctrl_unpack'
406			endif
407			read(cunit) (cbuff(ii), ii=1,cbuffindex)
408			endif
409			#endif
410
411			cbuffindex = 0
412			DO j=1,Ny
413			DO i=1,Nx
414			if (msk2d_buf_glo(i,j) .ne. 0. ) then
415			cbuffindex = cbuffindex + 1
416			fld2d_buf_glo(i,j) = cbuff(cbuffindex)
417			endif
418			ENDDO
419			ENDDO
420
421			ENDIF
422			_END_MASTER( mythid )
423			_BARRIER
424
425			CALL BAR2( myThid )
426			CALL SCATTER_2D_R8( fld2d_buf_glo, fld2d_buf,
427			& Nx,Ny,.FALSE.,.TRUE.,myThid)
428			CALL BAR2( myThid )
429			CALL MDS_PASS_R8toRL( fld2d_buf, fld3dNodim,
430			& 1, k, nr, 0, 0, .TRUE., myThid )
431
432			enddo !do k = 1, nr
433
434			c-- 2.2: normalize field if needed
435			DO bj = myByLo(myThid), myByHi(myThid)
436			DO bi = myBxLo(myThid), myBxHi(myThid)
437			DO k=1,Nr
438			if ( doZscalePack ) then
439			delZnorm = (delR(1)/delR(k))**delZexp
440			else
441			delZnorm = 1. _d 0
442			endif
443			DO j=1,sNy
444			DO i=1,sNx
445			if (msk3d(i,j,k,bi,bj).EQ.0. _d 0) then
446			fld3dDim(i,j,k,bi,bj)=0. _d 0
447			fld3dNodim(i,j,k,bi,bj)=0. _d 0
448			else
449			#ifdef ALLOW_ADMTLM
450			nveccount = nveccount + 1
451			fld3dNodim(i,j,k,bi,bj)=phtmpadmtlm(nveccount)
452			#endif
453			#ifdef ALLOW_SMOOTH_CORREL3D
454			fld3dDim(i,j,k,bi,bj)=fld3dNodim(i,j,k,bi,bj)
455			#else
456			# ifndef ALLOW_NONDIMENSIONAL_CONTROL_IO
457			fld3dDim(i,j,k,bi,bj) = fld3dNodim(i,j,k,bi,bj)
458			# else
459			if (lxxadxx) then
460			fld3dDim(i,j,k,bi,bj) =
461			& fld3dNodim(i,j,k,bi,bj) * delZnorm
462			# ifdef CTRL_PACK_PRECISE
463			& / sqrt(wei3d(i,j,k,bi,bj))
464			# else
465			& / sqrt(weightfld(k,bi,bj))
466			# endif
467			else
468			fld3dDim(i,j,k,bi,bj) =
469			& fld3dNodim(i,j,k,bi,bj) / delZnorm
470			# ifdef CTRL_PACK_PRECISE
471			& * sqrt(wei3d(i,j,k,bi,bj))
472			# else
473			& * sqrt(weightfld(k,bi,bj))
474			# endif
475			endif
476			# endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
477			#endif /* ALLOW_SMOOTH_CORREL3D */
478			endif
479			ENDDO
480			ENDDO
481			ENDDO
482			ENDDO
483			ENDDO
484
485			c-- 2.3:
486			if ( doPackDiag ) then
487			c error: twice the same one
488			call WRITE_REC_3D_RL( cfile2, ctrlprec,
489			& Nr, fld3dNodim, irec, 0, mythid)
490			call WRITE_REC_3D_RL( cfile3, ctrlprec,
491			& Nr, fld3dDim, irec, 0, mythid)
492			endif
493
494			c-- 2.4:
495			call WRITE_REC_3D_RL( fname, ctrlprec,
496			& Nr, fld3dDim, irec, 0, mythid)
497
498
499			enddo !do irec = 1, ncvarrecs(ivartype)
500
501			# endif /* ALLOW_PACKUNPACK_METHOD2 */
502			# endif /* EXCLUDE_CTRL_PACK */
503
504	heimbach	1.2	return
505			end
506