pkg/ctrl/ctrl_set_pack_yz.F

C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_set_pack_yz.F,v 1.12 2008/01/23 21:35:36 heimbach Exp $
C $Name:  $

#include "CTRL_CPPOPTIONS.h"

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

c     ==================================================================
c     SUBROUTINE ctrl_set_pack_yz
c     ==================================================================
c
c     o Compress the control vector such that only ocean points are
c       written to file.
c
c     o Open boundary packing added :
c          gebbie@mit.edu, 18-Mar-2003
c
c     changed: heimbach@mit.edu 17-Jun-2003
c              merged Armin's changes to replace write of
c              nr * globfld2d by 1 * globfld3d
c              (ad hoc fix to speed up global I/O)
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 ==

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

#ifndef EXCLUDE_CTRL_PACK
c     == local variables ==

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

      integer cbuffindex
cgg(
      integer igg
      integer reclen, irectrue
      integer cunit2, cunit3
      _RL     gg
      character*(80) weightname
      character*(80) cfile2, cfile3
cgg)
      real*4  cbuff      ( nsx*npx*sny*nsy*npy )
      real*4  globfldtmp2( nsx,npx,sny,nsy,npy )
      real*4  globfldtmp3( nsx,npx,sny,nsy,npy )
      _RL     globfldyz  ( nsx,npx,sny,nsy,npy,nr )
      _RL     globfld3d  ( snx,nsx,npx,sny,nsy,npy,nr )
      _RL     globmskyz  ( nsx,npx,sny,nsy,npy,nr,nobcs )
#ifdef CTRL_PACK_PRECISE
      _RL     weightfldyz( nsx,npx,sny,nsy,npy,nr,nobcs )
#endif

c     == external ==

      integer  ilnblnk
      external ilnblnk

c     == end of interface ==

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

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
                        globfldyz  (bi,ip,j,bj,jp,k) = 0. _d 0
                        globfldtmp2(bi,ip,j,bj,jp)   = 0.
                        globfldtmp3(bi,ip,j,bj,jp)   = 0.
                        do iobcs=1,nobcs
                           globmskyz(bi,ip,j,bj,jp,k,iobcs) = 0. _d 0
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo
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
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo

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

      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_pack_nonout_ctrl_',
     &           ivartype, '_', optimcycle, '.bin'
            write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
     &           'diag_pack_dimout_ctrl_',
     &           ivartype, '_', optimcycle, '.bin'
         else
            write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
     &           'diag_pack_nonout_grad_',
     &           ivartype, '_', optimcycle, '.bin'
            write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
     &           'diag_pack_dimout_grad_',
     &           ivartype, '_', optimcycle, '.bin'
         endif

         reclen = 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

      do iobcs=1,nobcs
         call MDSREADFIELD_YZ_GL(
     &        masktype, ctrlprec, 'RL',
     &        Nr, globmskyz(1,1,1,1,1,1,iobcs), iobcs,mythid)
#ifdef CTRL_PACK_PRECISE
         il=ilnblnk( weighttype)
         write(weightname(1:80),'(80a)') ' '
         write(weightname(1:80),'(a)') weighttype(1:il)
         call MDSREADFIELD_YZ_GL(
     &       weightname, ctrlprec, 'RL',
     &       Nr, weightfldyz(1,1,1,1,1,1,iobcs), iobcs, mythid)
CGG   One special exception: barotropic velocity should be nondimensionalized
cgg   differently. Probably introduce new variable.
         if (iobcs .eq. 3 .or. iobcs .eq. 4) then
            k = 1
            do jp = 1,nPy
               do bj = jtlo,jthi
                  do j = jmin,jmax
                     do ip = 1,nPx
                        do bi = itlo,ithi
cph                           weightfldyz(bi,ip,j,bj,jp,k,iobcs) = wbaro
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         endif
#endif
      enddo

      if ( useSingleCPUio ) then
C     MDSREADFIELD_YZ_GL does not know about useSingleCPUio, so the faster
C     method that works for .not.useSingleCPUio cannot be used
       nrec_nl = 0
      else
       nrec_nl = int(ncvarrecs(ivartype)/Nx)
      endif
      do irec = 1, nrec_nl
cgg       do iobcs = 1, nobcs
cgg    Need to solve for what iobcs would have been.

         call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
     &        nr, globfld3D, irec, mythid)

         do i=1,snx
            iobcs= mod((irec-1)*snx+i-1,nobcs)+1

CGG   One special exception: barotropic velocity should be nondimensionalized
cgg   differently. Probably introduce new variable.
            if (iobcs .eq. 3 .or. iobcs .eq. 4) then
               k = 1
               do jp = 1,nPy
                  do bj = jtlo,jthi
                     do j = jmin,jmax
                        do ip = 1,nPx
                           do bi = itlo,ithi
#ifdef NO_CONTROL_BAROTROPIC_VELOCITY
                              if (.not. lxxadxx) then
cgg    Get rid of any sensitivity to barotropic velocity.
                                 globfld3d(i,bi,ip,j,bj,jp,k) = 0.
                              endif
#endif
                           enddo
                        enddo
                     enddo
                  enddo
               enddo
            endif

            write(cunit) ncvarindex(ivartype)
            write(cunit) 1
            write(cunit) 1
            do k = 1,nr
             irectrue = (irec-1)*nobcs*nr + (iobcs-1)*nr + k
             cbuffindex = 0
             do jp = 1,nPy
              do bj = jtlo,jthi
               do ip = 1,nPx
                do bi = itlo,ithi
                 do j = jmin,jmax
                  ii=mod ( (i-1)*nr*sny+(k-1)*sny+j-1      , snx ) + 1
                  jj=mod( ((i-1)*nr*sny+(k-1)*sny+j-1)/snx , sny ) + 1
                  kk=int((i-1)*nr*sny+(k-1)*sny+j-1)/(snx*sny) + 1
                  if (globmskyz(bi,ip,j,bj,jp,k,iobcs)  .ne. 0. ) then
                     cbuffindex = cbuffindex + 1
cph(
                     globfldtmp3(bi,ip,j,bj,jp) =
     &                    globfld3d(ii,bi,ip,jj,bj,jp,kk)
cph)
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                     if (lxxadxx) then
                        cbuff(cbuffindex) =
     &                       globfld3d(ii,bi,ip,jj,bj,jp,kk) *
#ifdef CTRL_PACK_PRECISE
     &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
#else
     &                       sqrt(weightfld(k,iobcs))
#endif
                     else
                        cbuff(cbuffindex) =
     &                       globfld3d(ii,bi,ip,jj,bj,jp,kk) /
#ifdef CTRL_PACK_PRECISE
     &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
#else
     &                       sqrt(weightfld(k,iobcs))
#endif
                     endif
cph(
                     globfldtmp2(bi,ip,j,bj,jp) = cbuff(cbuffindex)
cph)
#else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
                     cbuff(cbuffindex) = globfld3d(ii,bi,ip,jj,bj,jp,kk)
#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                  endif
                 enddo
                enddo
               enddo
              enddo
             enddo
c           --> check cbuffindex.
             if ( cbuffindex .gt. 0) then
                write(cunit) cbuffindex
                write(cunit) k
                write(cunit) (cbuff(ii), ii=1,cbuffindex)
             endif
c
            if ( doPackDiag ) then
               write(cunit2,rec=irectrue) globfldtmp2
               write(cunit3,rec=irectrue) globfldtmp3
            endif
c
c     -- end of k loop --
            enddo
c     -- end of iobcs loop --
cgg       enddo
c     -- end of i loop --
         enddo
c     -- end of irec loop --
      enddo

      do irec = nrec_nl*nx+1, ncvarrecs(ivartype)
cgg       do iobcs = 1, nobcs
cgg    Need to solve for what iobcs would have been.
         iobcs= mod(irec-1,nobcs)+1

         call MDSREADFIELD_YZ_GL( fname, ctrlprec, 'RL',
     &        nr, globfldyz, irec, mythid)

CGG   One special exception: barotropic velocity should be nondimensionalized
cgg   differently. Probably introduce new variable.
         if (iobcs .eq. 3 .or. iobcs .eq. 4) then
            k = 1
            do jp = 1,nPy
               do bj = jtlo,jthi
                  do j = jmin,jmax
                     do ip = 1,nPx
                        do bi = itlo,ithi
#ifdef NO_CONTROL_BAROTROPIC_VELOCITY
                           if (.not. lxxadxx) then
cgg    Get rid of any sensitivity to barotropic velocity.
                              globfldyz(bi,ip,j,bj,jp,k) = 0.
                           endif
#endif
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         endif

         write(cunit) ncvarindex(ivartype)
         write(cunit) 1
         write(cunit) 1
         do k = 1,nr
            irectrue = (irec-1)*nobcs*nr + (iobcs-1)*nr + k
            cbuffindex = 0
            do jp = 1,nPy
             do bj = jtlo,jthi
              do ip = 1,nPx
               do bi = itlo,ithi
                do j = jmin,jmax
                 if (globmskyz(bi,ip,j,bj,jp,k,iobcs)  .ne. 0. ) then
                     cbuffindex = cbuffindex + 1
cph(
                     globfldtmp3(bi,ip,j,bj,jp) =
     &                    globfldyz(bi,ip,j,bj,jp,k)
cph)
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                     if (lxxadxx) then
                        cbuff(cbuffindex) =
     &                       globfldyz(bi,ip,j,bj,jp,k) *
#ifdef CTRL_PACK_PRECISE
     &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
#else
     &                       sqrt(weightfld(k,iobcs))
#endif
                     else
                        cbuff(cbuffindex) =
     &                       globfldyz(bi,ip,j,bj,jp,k) /
#ifdef CTRL_PACK_PRECISE
     &                       sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
#else
     &                       sqrt(weightfld(k,iobcs))
#endif
                     endif
cph(
                     globfldtmp2(bi,ip,j,bj,jp) = cbuff(cbuffindex)
cph)
#else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
                     cbuff(cbuffindex) = globfldyz(bi,ip,j,bj,jp,k)
#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                 endif
                enddo
               enddo
              enddo
             enddo
            enddo
c           --> check cbuffindex.
            if ( cbuffindex .gt. 0) then
               write(cunit) cbuffindex
               write(cunit) k
               write(cunit) (cbuff(ii), ii=1,cbuffindex)
            endif
c
            if ( doPackDiag ) then
               write(cunit2,rec=irectrue) globfldtmp2
               write(cunit3,rec=irectrue) globfldtmp3
            endif
c
c     -- end of k loop --
         enddo
c     -- end of iobcs loop --
cgg       enddo
c     -- end of irec loop --
      enddo

      _END_MASTER( mythid )

#endif

      return
      end

1	mlosch	1.13	C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_set_pack_yz.F,v 1.12 2008/01/23 21:35:36 heimbach Exp $
2	jmc	1.11	C $Name: $
3	heimbach	1.2
4			#include "CTRL_CPPOPTIONS.h"
5
6			subroutine ctrl_set_pack_yz(
7	heimbach	1.3	& cunit, ivartype, fname, masktype, weighttype,
8	heimbach	1.2	& weightfld, lxxadxx, mythid)
9
10			c ==================================================================
11			c SUBROUTINE ctrl_set_pack_yz
12			c ==================================================================
13			c
14			c o Compress the control vector such that only ocean points are
15			c written to file.
16			c
17	heimbach	1.3	c o Open boundary packing added :
18			c gebbie@mit.edu, 18-Mar-2003
19			c
20			c changed: heimbach@mit.edu 17-Jun-2003
21			c merged Armin's changes to replace write of
22			c nr * globfld2d by 1 * globfld3d
23			c (ad hoc fix to speed up global I/O)
24			c
25	heimbach	1.2	c ==================================================================
26
27			implicit none
28
29			c == global variables ==
30
31			#include "EEPARAMS.h"
32			#include "SIZE.h"
33			#include "PARAMS.h"
34			#include "GRID.h"
35
36			#include "ctrl.h"
37			#include "optim.h"
38
39			c == routine arguments ==
40
41			integer cunit
42			integer ivartype
43			character*( 80) fname
44			character*( 9) masktype
45	heimbach	1.3	character*( 80) weighttype
46	heimbach	1.2	_RL weightfld( nr,nobcs )
47			logical lxxadxx
48			integer mythid
49
50	heimbach	1.10	#ifndef EXCLUDE_CTRL_PACK
51	heimbach	1.2	c == local variables ==
52
53			integer bi,bj
54			integer ip,jp
55			integer i,j,k
56	heimbach	1.7	integer ii,jj,kk
57	heimbach	1.2	integer il
58	heimbach	1.3	integer irec,iobcs,nrec_nl
59	heimbach	1.2	integer itlo,ithi
60			integer jtlo,jthi
61			integer jmin,jmax
62			integer imin,imax
63
64			integer cbuffindex
65			cgg(
66			integer igg
67	heimbach	1.12	integer reclen, irectrue
68			integer cunit2, cunit3
69	heimbach	1.2	_RL gg
70	heimbach	1.3	character*(80) weightname
71	heimbach	1.12	character*(80) cfile2, cfile3
72	heimbach	1.2	cgg)
73	heimbach	1.12	real4 cbuff ( nsxnpxsnynsy*npy )
74			real*4 globfldtmp2( nsx,npx,sny,nsy,npy )
75			real*4 globfldtmp3( nsx,npx,sny,nsy,npy )
76	heimbach	1.2	_RL globfldyz ( nsx,npx,sny,nsy,npy,nr )
77	heimbach	1.3	_RL globfld3d ( snx,nsx,npx,sny,nsy,npy,nr )
78			_RL globmskyz ( nsx,npx,sny,nsy,npy,nr,nobcs )
79			#ifdef CTRL_PACK_PRECISE
80			_RL weightfldyz( nsx,npx,sny,nsy,npy,nr,nobcs )
81			#endif
82	heimbach	1.2
83			c == external ==
84
85			integer ilnblnk
86			external ilnblnk
87
88			c == end of interface ==
89
90			jtlo = 1
91			jthi = nsy
92			itlo = 1
93			ithi = nsx
94			jmin = 1
95			jmax = sny
96			imin = 1
97			imax = snx
98
99			c Initialise temporary file
100			do k = 1,nr
101			do jp = 1,nPy
102			do bj = jtlo,jthi
103			do j = jmin,jmax
104			do ip = 1,nPx
105			do bi = itlo,ithi
106	heimbach	1.12	globfldyz (bi,ip,j,bj,jp,k) = 0. _d 0
107			globfldtmp2(bi,ip,j,bj,jp) = 0.
108			globfldtmp3(bi,ip,j,bj,jp) = 0.
109	heimbach	1.3	do iobcs=1,nobcs
110			globmskyz(bi,ip,j,bj,jp,k,iobcs) = 0. _d 0
111			enddo
112			enddo
113			enddo
114			enddo
115			enddo
116			enddo
117			enddo
118			c Initialise temporary file
119			do k = 1,nr
120			do jp = 1,nPy
121			do bj = jtlo,jthi
122			do j = jmin,jmax
123			do ip = 1,nPx
124			do bi = itlo,ithi
125			do i = imin,imax
126			globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
127			enddo
128	heimbach	1.2	enddo
129			enddo
130			enddo
131			enddo
132			enddo
133			enddo
134
135			c-- Only the master thread will do I/O.
136			_BEGIN_MASTER( mythid )
137
138	heimbach	1.12	if ( doPackDiag ) then
139			write(cfile2(1:80),'(80a)') ' '
140			write(cfile3(1:80),'(80a)') ' '
141			if ( lxxadxx ) then
142			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
143			& 'diag_pack_nonout_ctrl_',
144			& ivartype, '_', optimcycle, '.bin'
145			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
146			& 'diag_pack_dimout_ctrl_',
147			& ivartype, '_', optimcycle, '.bin'
148			else
149			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
150			& 'diag_pack_nonout_grad_',
151			& ivartype, '_', optimcycle, '.bin'
152			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
153			& 'diag_pack_dimout_grad_',
154			& ivartype, '_', optimcycle, '.bin'
155			endif
156
157			reclen = nsxnpxsnynsynpy*4
158			call mdsfindunit( cunit2, mythid )
159			open( cunit2, file=cfile2, status='unknown',
160			& access='direct', recl=reclen )
161			call mdsfindunit( cunit3, mythid )
162			open( cunit3, file=cfile3, status='unknown',
163			& access='direct', recl=reclen )
164			endif
165
166	heimbach	1.3	do iobcs=1,nobcs
167	jmc	1.11	call MDSREADFIELD_YZ_GL(
168	heimbach	1.3	& masktype, ctrlprec, 'RL',
169			& Nr, globmskyz(1,1,1,1,1,1,iobcs), iobcs,mythid)
170			#ifdef CTRL_PACK_PRECISE
171			il=ilnblnk( weighttype)
172			write(weightname(1:80),'(80a)') ' '
173			write(weightname(1:80),'(a)') weighttype(1:il)
174			call MDSREADFIELD_YZ_GL(
175			& weightname, ctrlprec, 'RL',
176			& Nr, weightfldyz(1,1,1,1,1,1,iobcs), iobcs, mythid)
177			CGG One special exception: barotropic velocity should be nondimensionalized
178			cgg differently. Probably introduce new variable.
179			if (iobcs .eq. 3 .or. iobcs .eq. 4) then
180			k = 1
181			do jp = 1,nPy
182			do bj = jtlo,jthi
183			do j = jmin,jmax
184			do ip = 1,nPx
185			do bi = itlo,ithi
186	heimbach	1.9	cph weightfldyz(bi,ip,j,bj,jp,k,iobcs) = wbaro
187	heimbach	1.3	enddo
188			enddo
189			enddo
190			enddo
191			enddo
192			endif
193			#endif
194			enddo
195	heimbach	1.2
196	mlosch	1.13	if ( useSingleCPUio ) then
197			C MDSREADFIELD_YZ_GL does not know about useSingleCPUio, so the faster
198			C method that works for .not.useSingleCPUio cannot be used
199			nrec_nl = 0
200			else
201			nrec_nl = int(ncvarrecs(ivartype)/Nx)
202			endif
203	heimbach	1.3	do irec = 1, nrec_nl
204	heimbach	1.2	cgg do iobcs = 1, nobcs
205			cgg Need to solve for what iobcs would have been.
206
207	heimbach	1.3	call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
208			& nr, globfld3D, irec, mythid)
209
210			do i=1,snx
211			iobcs= mod((irec-1)*snx+i-1,nobcs)+1
212
213			CGG One special exception: barotropic velocity should be nondimensionalized
214			cgg differently. Probably introduce new variable.
215			if (iobcs .eq. 3 .or. iobcs .eq. 4) then
216			k = 1
217			do jp = 1,nPy
218			do bj = jtlo,jthi
219			do j = jmin,jmax
220			do ip = 1,nPx
221			do bi = itlo,ithi
222			#ifdef NO_CONTROL_BAROTROPIC_VELOCITY
223			if (.not. lxxadxx) then
224			cgg Get rid of any sensitivity to barotropic velocity.
225			globfld3d(i,bi,ip,j,bj,jp,k) = 0.
226			endif
227			#endif
228			enddo
229			enddo
230			enddo
231			enddo
232			enddo
233			endif
234
235			write(cunit) ncvarindex(ivartype)
236			write(cunit) 1
237			write(cunit) 1
238			do k = 1,nr
239	heimbach	1.12	irectrue = (irec-1)nobcsnr + (iobcs-1)*nr + k
240	heimbach	1.3	cbuffindex = 0
241			do jp = 1,nPy
242			do bj = jtlo,jthi
243			do ip = 1,nPx
244			do bi = itlo,ithi
245			do j = jmin,jmax
246	heimbach	1.12	ii=mod ( (i-1)nrsny+(k-1)*sny+j-1 , snx ) + 1
247			jj=mod( ((i-1)nrsny+(k-1)*sny+j-1)/snx , sny ) + 1
248			kk=int((i-1)nrsny+(k-1)sny+j-1)/(snxsny) + 1
249	heimbach	1.3	if (globmskyz(bi,ip,j,bj,jp,k,iobcs) .ne. 0. ) then
250			cbuffindex = cbuffindex + 1
251	heimbach	1.12	cph(
252			globfldtmp3(bi,ip,j,bj,jp) =
253			& globfld3d(ii,bi,ip,jj,bj,jp,kk)
254			cph)
255	heimbach	1.3	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
256			if (lxxadxx) then
257	jmc	1.11	cbuff(cbuffindex) =
258	heimbach	1.7	& globfld3d(ii,bi,ip,jj,bj,jp,kk) *
259	heimbach	1.3	#ifdef CTRL_PACK_PRECISE
260			& sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
261			#else
262			& sqrt(weightfld(k,iobcs))
263			#endif
264			else
265	jmc	1.11	cbuff(cbuffindex) =
266	heimbach	1.7	& globfld3d(ii,bi,ip,jj,bj,jp,kk) /
267	heimbach	1.3	#ifdef CTRL_PACK_PRECISE
268			& sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
269			#else
270			& sqrt(weightfld(k,iobcs))
271			#endif
272			endif
273	heimbach	1.12	cph(
274			globfldtmp2(bi,ip,j,bj,jp) = cbuff(cbuffindex)
275			cph)
276	heimbach	1.3	#else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
277	heimbach	1.7	cbuff(cbuffindex) = globfld3d(ii,bi,ip,jj,bj,jp,kk)
278	heimbach	1.3	#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
279			endif
280			enddo
281			enddo
282			enddo
283			enddo
284			enddo
285			c --> check cbuffindex.
286			if ( cbuffindex .gt. 0) then
287			write(cunit) cbuffindex
288			write(cunit) k
289			write(cunit) (cbuff(ii), ii=1,cbuffindex)
290			endif
291			c
292	heimbach	1.12	if ( doPackDiag ) then
293			write(cunit2,rec=irectrue) globfldtmp2
294			write(cunit3,rec=irectrue) globfldtmp3
295			endif
296			c
297	heimbach	1.3	c -- end of k loop --
298			enddo
299			c -- end of iobcs loop --
300			cgg enddo
301			c -- end of i loop --
302			enddo
303			c -- end of irec loop --
304			enddo
305
306	mlosch	1.13	do irec = nrec_nl*nx+1, ncvarrecs(ivartype)
307	heimbach	1.3	cgg do iobcs = 1, nobcs
308			cgg Need to solve for what iobcs would have been.
309			iobcs= mod(irec-1,nobcs)+1
310	heimbach	1.2
311			call MDSREADFIELD_YZ_GL( fname, ctrlprec, 'RL',
312			& nr, globfldyz, irec, mythid)
313
314	heimbach	1.3	CGG One special exception: barotropic velocity should be nondimensionalized
315			cgg differently. Probably introduce new variable.
316			if (iobcs .eq. 3 .or. iobcs .eq. 4) then
317			k = 1
318			do jp = 1,nPy
319			do bj = jtlo,jthi
320			do j = jmin,jmax
321			do ip = 1,nPx
322			do bi = itlo,ithi
323			#ifdef NO_CONTROL_BAROTROPIC_VELOCITY
324			if (.not. lxxadxx) then
325			cgg Get rid of any sensitivity to barotropic velocity.
326			globfldyz(bi,ip,j,bj,jp,k) = 0.
327			endif
328			#endif
329			enddo
330			enddo
331			enddo
332			enddo
333			enddo
334			endif
335
336	heimbach	1.2	write(cunit) ncvarindex(ivartype)
337			write(cunit) 1
338			write(cunit) 1
339			do k = 1,nr
340	heimbach	1.12	irectrue = (irec-1)nobcsnr + (iobcs-1)*nr + k
341	heimbach	1.2	cbuffindex = 0
342			do jp = 1,nPy
343			do bj = jtlo,jthi
344			do ip = 1,nPx
345			do bi = itlo,ithi
346			do j = jmin,jmax
347	heimbach	1.3	if (globmskyz(bi,ip,j,bj,jp,k,iobcs) .ne. 0. ) then
348	heimbach	1.2	cbuffindex = cbuffindex + 1
349	heimbach	1.12	cph(
350			globfldtmp3(bi,ip,j,bj,jp) =
351			& globfldyz(bi,ip,j,bj,jp,k)
352			cph)
353	heimbach	1.2	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
354			if (lxxadxx) then
355	jmc	1.11	cbuff(cbuffindex) =
356	heimbach	1.2	& globfldyz(bi,ip,j,bj,jp,k) *
357	heimbach	1.3	#ifdef CTRL_PACK_PRECISE
358			& sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
359			#else
360	heimbach	1.2	& sqrt(weightfld(k,iobcs))
361	heimbach	1.3	#endif
362	heimbach	1.2	else
363	jmc	1.11	cbuff(cbuffindex) =
364	heimbach	1.2	& globfldyz(bi,ip,j,bj,jp,k) /
365	heimbach	1.3	#ifdef CTRL_PACK_PRECISE
366			& sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
367			#else
368	heimbach	1.2	& sqrt(weightfld(k,iobcs))
369	heimbach	1.3	#endif
370	heimbach	1.2	endif
371	heimbach	1.12	cph(
372			globfldtmp2(bi,ip,j,bj,jp) = cbuff(cbuffindex)
373			cph)
374	heimbach	1.3	#else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
375	heimbach	1.2	cbuff(cbuffindex) = globfldyz(bi,ip,j,bj,jp,k)
376	heimbach	1.3	#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
377	heimbach	1.2	endif
378			enddo
379			enddo
380			enddo
381			enddo
382			enddo
383			c --> check cbuffindex.
384			if ( cbuffindex .gt. 0) then
385			write(cunit) cbuffindex
386			write(cunit) k
387			write(cunit) (cbuff(ii), ii=1,cbuffindex)
388			endif
389	heimbach	1.3	c
390	heimbach	1.12	if ( doPackDiag ) then
391			write(cunit2,rec=irectrue) globfldtmp2
392			write(cunit3,rec=irectrue) globfldtmp3
393			endif
394			c
395	heimbach	1.3	c -- end of k loop --
396	heimbach	1.2	enddo
397			c -- end of iobcs loop --
398			cgg enddo
399			c -- end of irec loop --
400			enddo
401
402			_END_MASTER( mythid )
403
404	heimbach	1.10	#endif
405
406	heimbach	1.2	return
407			end
408