pkg/ctrl/ctrl_set_pack_xyz.F


#include "CTRL_CPPOPTIONS.h"

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

c     ==================================================================
c     SUBROUTINE ctrl_set_pack_xyz
c     ==================================================================
c
c     o Compress the control vector such that only ocean points are
c       written to file.
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 ==

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

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_PACK_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*(80) weightname

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

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

#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 )

      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 = 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_PACK_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)

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

         write(cunit) ncvarindex(ivartype)
         write(cunit) 1
         write(cunit) 1
         do k = 1, nr
         irectrue = (irec-1)*nr + k
            if ( doZscalePack ) then
cph               delZnorm = SQRT(delR(1)/delR(k))
               delZnorm = delR(1)/delR(k)
            else
               delZnorm = 1. _d 0
            endif
            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
cph(
                     globfldtmp3(i,bi,ip,j,bj,jp) =
     &                    globfld3d(i,bi,ip,j,bj,jp,k)
cph)
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                     if (lxxadxx) then
                        cbuff(cbuffindex) = delZnorm
     &                       * globfld3d(i,bi,ip,j,bj,jp,k)
# ifdef CTRL_PACK_PRECISE
     &                       * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
# else
     &                       * sqrt(weightfld(k,bi,bj))
# endif
                     else
                        cbuff(cbuffindex) = delZnorm
     &                       * globfld3d(i,bi,ip,j,bj,jp,k)
# ifdef CTRL_PACK_PRECISE
     &                       / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
# else
     &                       / sqrt(weightfld(k,bi,bj))
# endif
                     endif
cph(
                     globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
cph)
#else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
                     cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                  endif
                 enddo
                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
         enddo
c
c     -- end of irec loop --
      enddo

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

      return
      end

1	heimbach	1.2
2			#include "CTRL_CPPOPTIONS.h"
3
4			subroutine ctrl_set_pack_xyz(
5	heimbach	1.4	& cunit, ivartype, fname, masktype, weighttype,
6	heimbach	1.2	& weightfld, lxxadxx, mythid)
7
8			c ==================================================================
9			c SUBROUTINE ctrl_set_pack_xyz
10			c ==================================================================
11			c
12			c o Compress the control vector such that only ocean points are
13			c written to file.
14			c
15	heimbach	1.4	c o Use a more precise nondimensionalization that depends on (x,y)
16			c Added weighttype to the argument list so that I can geographically
17			c vary the nondimensionalization.
18			c gebbie@mit.edu, 18-Mar-2003
19			c
20	heimbach	1.2	c ==================================================================
21
22			implicit none
23
24			c == global variables ==
25
26			#include "EEPARAMS.h"
27			#include "SIZE.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30
31			#include "ctrl.h"
32			#include "optim.h"
33
34			c == routine arguments ==
35
36			integer cunit
37			integer ivartype
38			character*( 80) fname
39	heimbach	1.11	character*( 9) masktype
40	heimbach	1.4	character*( 80) weighttype
41	heimbach	1.2	_RL weightfld( nr,nsx,nsy )
42			logical lxxadxx
43			integer mythid
44
45			c == local variables ==
46
47			integer bi,bj
48			integer ip,jp
49			integer i,j,k
50			integer ii
51			integer il
52			integer irec
53			integer itlo,ithi
54			integer jtlo,jthi
55			integer jmin,jmax
56			integer imin,imax
57
58			integer cbuffindex
59
60			_RL globmsk ( snx,nsx,npx,sny,nsy,npy,nr )
61			_RL globfld3d( snx,nsx,npx,sny,nsy,npy,nr )
62	heimbach	1.4	#ifdef CTRL_PACK_PRECISE
63			_RL weightfld3d( snx,nsx,npx,sny,nsy,npy,nr )
64			#endif
65	heimbach	1.9	real4 cbuff ( snxnsxnpxsnynsynpy )
66			real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
67			real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )
68	heimbach	1.4
69			character*(80) weightname
70	heimbach	1.2
71	heimbach	1.9	_RL delZnorm
72			integer reclen, irectrue
73			integer cunit2, cunit3
74			character*(80) cfile2, cfile3
75
76	heimbach	1.2	c == external ==
77
78			integer ilnblnk
79			external ilnblnk
80
81			c == end of interface ==
82
83			jtlo = 1
84			jthi = nsy
85			itlo = 1
86			ithi = nsx
87			jmin = 1
88			jmax = sny
89			imin = 1
90			imax = snx
91
92	heimbach	1.9	#ifdef CTRL_DELZNORM
93			delZnorm = 0.
94			do k = 1, Nr
95			delZnorm = delZnorm + delR(k)/FLOAT(Nr)
96			enddo
97			#endif
98
99	heimbach	1.2	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			do i = imin,imax
107	heimbach	1.9	globfld3d (i,bi,ip,j,bj,jp,k) = 0. _d 0
108			globmsk (i,bi,ip,j,bj,jp,k) = 0. _d 0
109			globfldtmp2(i,bi,ip,j,bj,jp) = 0.
110			globfldtmp3(i,bi,ip,j,bj,jp) = 0.
111	heimbach	1.2	enddo
112			enddo
113			enddo
114			enddo
115			enddo
116			enddo
117			enddo
118
119			c-- Only the master thread will do I/O.
120			_BEGIN_MASTER( mythid )
121
122	heimbach	1.9	if ( doPackDiag ) then
123			write(cfile2(1:80),'(80a)') ' '
124			write(cfile3(1:80),'(80a)') ' '
125			if ( lxxadxx ) then
126			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
127			& 'diag_pack_nonout_ctrl_',
128			& ivartype, '_', optimcycle, '.bin'
129			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
130			& 'diag_pack_dimout_ctrl_',
131			& ivartype, '_', optimcycle, '.bin'
132			else
133			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
134			& 'diag_pack_nonout_grad_',
135			& ivartype, '_', optimcycle, '.bin'
136			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
137			& 'diag_pack_dimout_grad_',
138			& ivartype, '_', optimcycle, '.bin'
139			endif
140
141			reclen = FLOAT(snxnsxnpxsnynsynpy4)
142			call mdsfindunit( cunit2, mythid )
143			open( cunit2, file=cfile2, status='unknown',
144			& access='direct', recl=reclen )
145			call mdsfindunit( cunit3, mythid )
146			open( cunit3, file=cfile3, status='unknown',
147			& access='direct', recl=reclen )
148			endif
149
150	heimbach	1.4	#ifdef CTRL_PACK_PRECISE
151			il=ilnblnk( weighttype)
152			write(weightname(1:80),'(80a)') ' '
153			write(weightname(1:80),'(a)') weighttype(1:il)
154
155			call MDSREADFIELD_3D_GL(
156			& weightname, ctrlprec, 'RL',
157			& Nr, weightfld3d, 1, mythid)
158			#endif
159
160	heimbach	1.2	call MDSREADFIELD_3D_GL(
161			& masktype, ctrlprec, 'RL',
162			& Nr, globmsk, 1, mythid)
163
164			do irec = 1, ncvarrecs(ivartype)
165
166			call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
167			& Nr, globfld3d, irec, mythid)
168
169			write(cunit) ncvarindex(ivartype)
170			write(cunit) 1
171			write(cunit) 1
172			do k = 1, nr
173	heimbach	1.9	irectrue = (irec-1)*nr + k
174	heimbach	1.10	if ( doZscalePack ) then
175	heimbach	1.12	cph delZnorm = SQRT(delR(1)/delR(k))
176			delZnorm = delR(1)/delR(k)
177	heimbach	1.10	else
178			delZnorm = 1. _d 0
179			endif
180	heimbach	1.2	cbuffindex = 0
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			do i = imin,imax
187	heimbach	1.9	if (globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
188	heimbach	1.2	cbuffindex = cbuffindex + 1
189	heimbach	1.9	cph(
190			globfldtmp3(i,bi,ip,j,bj,jp) =
191			& globfld3d(i,bi,ip,j,bj,jp,k)
192			cph)
193	heimbach	1.2	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
194			if (lxxadxx) then
195	heimbach	1.10	cbuff(cbuffindex) = delZnorm
196			& * globfld3d(i,bi,ip,j,bj,jp,k)
197	heimbach	1.4	# ifdef CTRL_PACK_PRECISE
198	heimbach	1.9	& * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
199	heimbach	1.4	# else
200	heimbach	1.9	& * sqrt(weightfld(k,bi,bj))
201	heimbach	1.4	# endif
202	heimbach	1.2	else
203	heimbach	1.10	cbuff(cbuffindex) = delZnorm
204			& * globfld3d(i,bi,ip,j,bj,jp,k)
205	heimbach	1.4	# ifdef CTRL_PACK_PRECISE
206	heimbach	1.9	& / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
207	heimbach	1.4	# else
208	heimbach	1.9	& / sqrt(weightfld(k,bi,bj))
209	heimbach	1.4	# endif
210	heimbach	1.2	endif
211	heimbach	1.9	cph(
212			globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
213			cph)
214	heimbach	1.4	#else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
215	heimbach	1.2	cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
216	heimbach	1.4	#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
217	heimbach	1.2	endif
218			enddo
219			enddo
220			enddo
221			enddo
222			enddo
223			enddo
224			c --> check cbuffindex.
225			if ( cbuffindex .gt. 0) then
226			write(cunit) cbuffindex
227			write(cunit) k
228			write(cunit) (cbuff(ii), ii=1,cbuffindex)
229			endif
230	heimbach	1.9	c
231			if ( doPackDiag ) then
232			write(cunit2,rec=irectrue) globfldtmp2
233			write(cunit3,rec=irectrue) globfldtmp3
234			endif
235			c
236	heimbach	1.2	enddo
237			c
238			c -- end of irec loop --
239			enddo
240
241	heimbach	1.9	if ( doPackDiag ) then
242			close ( cunit2 )
243			close ( cunit3 )
244			endif
245
246	heimbach	1.2	_END_MASTER( mythid )
247
248			return
249			end
250