pkg/ctrl/ctrl_set_unpack_yz.F


#include "CTRL_CPPOPTIONS.h"

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

c     ==================================================================
c     SUBROUTINE ctrl_set_unpack_yz
c     ==================================================================
c
c     o Unpack the control vector such that land points are filled in.
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 )
      integer nwetglobal(nr,nobcs)
      integer mythid

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

      real*4     cbuff    ( 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_UNPACK_PRECISE
      _RL   weightfldyz( nsx,npx,sny,nsy,npy,nr,nobcs )
#endif

cgg(
      integer igg
      _RL     gg
      character*(80) weightname
cgg)

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

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

      do iobcs=1,nobcs
         call MDSREADFIELD_YZ_GL(
     &        masktype, ctrlprec, 'RL',
     &        Nr, globmskyz(1,1,1,1,1,1,iobcs), iobcs, mythid)
#ifdef CTRL_UNPACK_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
                           weightfldyz(bi,ip,j,bj,jp,k,iobcs) = wbaro
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         endif
#endif
      enddo

      nrec_nl=int(ncvarrecs(ivartype)/snx)
      do irec = 1, nrec_nl
cgg      do iobcs = 1, nobcs
cgg      And now back-calculate what iobcs should be.
         do i=1,snx
            iobcs= mod((irec-1)*snx+i-1,nobcs)+1

            read(cunit) filencvarindex(ivartype)
            if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
     &           then
               print *, 'ctrl_set_unpack_yz:WARNING: wrong ncvarindex ',
     &              filencvarindex(ivartype), ncvarindex(ivartype)
               STOP 'in S/R ctrl_unpack'
            endif
            read(cunit) filej
            read(cunit) filei
            do k = 1, Nr
               cbuffindex = nwetglobal(k,iobcs) 
               if ( cbuffindex .gt. 0 ) then
                  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
                  read(cunit) (cbuff(ii), ii=1,cbuffindex)
               endif
               cbuffindex = 0
               do jp = 1,nPy
                do bj = jtlo,jthi
                 do j = jmin,jmax
                  do ip = 1,nPx
                   do bi = itlo,ithi
                    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
                       globfld3d(ii,bi,ip,jj,bj,jp,kk) = 
     &                      cbuff(cbuffindex)
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                       globfld3d(ii,bi,ip,jj,bj,jp,kk) = 
     &                      globfld3d(ii,bi,ip,jj,bj,jp,kk)/
# ifdef CTRL_UNPACK_PRECISE
     &                      sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
# else
     &                      sqrt(weightfld(k,iobcs))
# endif
#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                    else
                       globfld3d(ii,bi,ip,jj,bj,jp,kk) = 0. _d 0
                    endif
                   enddo
                  enddo
                 enddo
                enddo
               enddo
c
c     -- end of k loop --
            enddo
c     -- end of i loop --
         enddo

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

c     -- end of iobcs loop -- This loop has been removed.
cgg     enddo
c     -- end of irec loop --
      enddo

      do irec = nrec_nl*snx+1,ncvarrecs(ivartype)
         iobcs= mod(irec-1,nobcs)+1

         read(cunit) filencvarindex(ivartype)
         if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
     &        then
            print *, 'ctrl_set_unpack_yz:WARNING: wrong ncvarindex ',
     &           filencvarindex(ivartype), ncvarindex(ivartype)
            STOP 'in S/R ctrl_unpack'
         endif
         read(cunit) filej
         read(cunit) filei
         do k = 1, Nr
            cbuffindex = nwetglobal(k,iobcs) 
            if ( cbuffindex .gt. 0 ) then
               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
               read(cunit) (cbuff(ii), ii=1,cbuffindex)
            endif
            cbuffindex = 0
            do jp = 1,nPy
             do bj = jtlo,jthi
              do j = jmin,jmax
               do ip = 1,nPx
                do bi = itlo,ithi
                  if ( globmskyz(bi,ip,j,bj,jp,k,iobcs) .ne. 0. ) then
                     cbuffindex = cbuffindex + 1
                     globfldyz(bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                     globfldyz(bi,ip,j,bj,jp,k) = 
     &                    globfldyz(bi,ip,j,bj,jp,k)/
# ifdef CTRL_UNPACK_PRECISE
     &                    sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
# else
     &                    sqrt(weightfld(k,iobcs))
# endif
#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
                  else
                     globfldyz(bi,ip,j,bj,jp,k) = 0. _d 0
                  endif
                enddo
               enddo
              enddo
             enddo
            enddo
c
c     -- end of k loop
         enddo
             
         call MDSWRITEFIELD_YZ_GL( fname, ctrlprec, 'RL',
     &                             Nr, globfldyz, irec,
     &                             optimcycle, mythid)

c     -- end of iobcs loop -- This loop has been removed.
cgg     enddo
c     -- end of irec loop --
      enddo

      _END_MASTER( mythid )

      return
      end

1	heimbach	1.2
2			#include "CTRL_CPPOPTIONS.h"
3
4			subroutine ctrl_set_unpack_yz(
5	heimbach	1.3	& cunit, ivartype, fname, masktype, weighttype,
6	heimbach	1.2	& weightfld, nwetglobal, mythid)
7
8			c ==================================================================
9			c SUBROUTINE ctrl_set_unpack_yz
10			c ==================================================================
11			c
12	heimbach	1.3	c o Unpack the control vector such that land points are filled in.
13			c
14			c o Open boundary packing added :
15			c gebbie@mit.edu, 18-Mar-2003
16			c
17			c changed: heimbach@mit.edu 17-Jun-2003
18			c merged Armin's changes to replace write of
19			c nr * globfld2d by 1 * globfld3d
20			c (ad hoc fix to speed up global I/O)
21	heimbach	1.2	c
22			c ==================================================================
23
24			implicit none
25
26			c == global variables ==
27
28			#include "EEPARAMS.h"
29			#include "SIZE.h"
30			#include "PARAMS.h"
31			#include "GRID.h"
32
33			#include "ctrl.h"
34			#include "optim.h"
35
36			c == routine arguments ==
37
38			integer cunit
39			integer ivartype
40			character*( 80) fname
41			character* (9) masktype
42	heimbach	1.3	character*( 80) weighttype
43	heimbach	1.2	_RL weightfld( nr,nobcs )
44			integer nwetglobal(nr,nobcs)
45			integer mythid
46
47			c == local variables ==
48
49			integer bi,bj
50			integer ip,jp
51			integer i,j,k
52	heimbach	1.7	integer ii,jj,kk
53	heimbach	1.2	integer il
54	heimbach	1.3	integer irec,iobcs,nrec_nl
55	heimbach	1.2	integer itlo,ithi
56			integer jtlo,jthi
57			integer jmin,jmax
58			integer imin,imax
59
60			integer cbuffindex
61
62	heimbach	1.7	real4 cbuff ( nsxnpxsnynsy*npy )
63	heimbach	1.2	_RL globfldyz( nsx,npx,sny,nsy,npy,nr )
64	heimbach	1.3	_RL globfld3d( snx,nsx,npx,sny,nsy,npy,nr )
65			_RL globmskyz( nsx,npx,sny,nsy,npy,nr,nobcs )
66			#ifdef CTRL_UNPACK_PRECISE
67			_RL weightfldyz( nsx,npx,sny,nsy,npy,nr,nobcs )
68			#endif
69	heimbach	1.2
70			cgg(
71			integer igg
72			_RL gg
73	heimbach	1.3	character*(80) weightname
74	heimbach	1.2	cgg)
75
76			c == external ==
77
78			integer ilnblnk
79			external ilnblnk
80
81			cc == 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			c Initialise temporary file
93			do k = 1,nr
94			do jp = 1,nPy
95			do bj = jtlo,jthi
96			do j = jmin,jmax
97			do ip = 1,nPx
98			do bi = itlo,ithi
99			globfldyz(bi,ip,j,bj,jp,k) = 0. _d 0
100	heimbach	1.3	do iobcs=1,nobcs
101			globmskyz(bi,ip,j,bj,jp,k,iobcs) = 0. _d 0
102			enddo
103			enddo
104			enddo
105			enddo
106			enddo
107			enddo
108			enddo
109			c Initialise temporary file
110			do k = 1,nr
111			do jp = 1,nPy
112			do bj = jtlo,jthi
113			do j = jmin,jmax
114			do ip = 1,nPx
115			do bi = itlo,ithi
116			do i = imin,imax
117			globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
118			enddo
119	heimbach	1.2	enddo
120			enddo
121			enddo
122			enddo
123			enddo
124			enddo
125
126			c-- Only the master thread will do I/O.
127			_BEGIN_MASTER( mythid )
128
129	heimbach	1.3	do iobcs=1,nobcs
130			call MDSREADFIELD_YZ_GL(
131			& masktype, ctrlprec, 'RL',
132			& Nr, globmskyz(1,1,1,1,1,1,iobcs), iobcs, mythid)
133			#ifdef CTRL_UNPACK_PRECISE
134			il=ilnblnk( weighttype)
135			write(weightname(1:80),'(80a)') ' '
136			write(weightname(1:80),'(a)') weighttype(1:il)
137			call MDSREADFIELD_YZ_GL(
138			& weightname, ctrlprec, 'RL',
139			& Nr, weightfldyz(1,1,1,1,1,1,iobcs), iobcs, mythid)
140			CGG One special exception: barotropic velocity should be nondimensionalized
141			cgg differently. Probably introduce new variable.
142			if (iobcs .eq. 3 .or. iobcs .eq. 4) then
143			k = 1
144			do jp = 1,nPy
145			do bj = jtlo,jthi
146			do j = jmin,jmax
147			do ip = 1,nPx
148			do bi = itlo,ithi
149			weightfldyz(bi,ip,j,bj,jp,k,iobcs) = wbaro
150			enddo
151			enddo
152			enddo
153			enddo
154			enddo
155			endif
156			#endif
157			enddo
158
159			nrec_nl=int(ncvarrecs(ivartype)/snx)
160			do irec = 1, nrec_nl
161	heimbach	1.2	cgg do iobcs = 1, nobcs
162			cgg And now back-calculate what iobcs should be.
163	heimbach	1.3	do i=1,snx
164			iobcs= mod((irec-1)*snx+i-1,nobcs)+1
165
166			read(cunit) filencvarindex(ivartype)
167			if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
168			& then
169			print *, 'ctrl_set_unpack_yz:WARNING: wrong ncvarindex ',
170			& filencvarindex(ivartype), ncvarindex(ivartype)
171			STOP 'in S/R ctrl_unpack'
172			endif
173			read(cunit) filej
174			read(cunit) filei
175			do k = 1, Nr
176			cbuffindex = nwetglobal(k,iobcs)
177			if ( cbuffindex .gt. 0 ) then
178			read(cunit) filencbuffindex
179			if (filencbuffindex .NE. cbuffindex) then
180			print *, 'WARNING: wrong cbuffindex ',
181			& filencbuffindex, cbuffindex
182			STOP 'in S/R ctrl_unpack'
183			endif
184			read(cunit) filek
185			if (filek .NE. k) then
186			print *, 'WARNING: wrong k ',
187			& filek, k
188			STOP 'in S/R ctrl_unpack'
189			endif
190			read(cunit) (cbuff(ii), ii=1,cbuffindex)
191			endif
192			cbuffindex = 0
193			do jp = 1,nPy
194			do bj = jtlo,jthi
195			do j = jmin,jmax
196			do ip = 1,nPx
197			do bi = itlo,ithi
198	heimbach	1.7	ii=mod((i-1)nrsny+(k-1)*sny+j-1,snx)+1
199			jj=mod(((i-1)nrsny+(k-1)*sny+j-1)/snx,sny)+1
200			kk=int((i-1)nrsny+(k-1)sny+j-1)/(snxsny)+1
201	heimbach	1.3	if ( globmskyz(bi,ip,j,bj,jp,k,iobcs) .ne. 0. ) then
202			cbuffindex = cbuffindex + 1
203	heimbach	1.7	globfld3d(ii,bi,ip,jj,bj,jp,kk) =
204	heimbach	1.4	& cbuff(cbuffindex)
205	heimbach	1.3	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
206	heimbach	1.7	globfld3d(ii,bi,ip,jj,bj,jp,kk) =
207			& globfld3d(ii,bi,ip,jj,bj,jp,kk)/
208	heimbach	1.3	# ifdef CTRL_UNPACK_PRECISE
209			& sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
210			# else
211			& sqrt(weightfld(k,iobcs))
212			# endif
213			#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
214			else
215	heimbach	1.7	globfld3d(ii,bi,ip,jj,bj,jp,kk) = 0. _d 0
216	heimbach	1.3	endif
217			enddo
218			enddo
219			enddo
220			enddo
221			enddo
222			c
223			c -- end of k loop --
224			enddo
225			c -- end of i loop --
226			enddo
227	heimbach	1.2
228	heimbach	1.3	call MDSWRITEFIELD_3d_GL( fname, ctrlprec, 'RL',
229			& Nr, globfld3d, irec,
230			& optimcycle, mythid)
231
232			c -- end of iobcs loop -- This loop has been removed.
233			cgg enddo
234			c -- end of irec loop --
235			enddo
236
237			do irec = nrec_nl*snx+1,ncvarrecs(ivartype)
238			iobcs= mod(irec-1,nobcs)+1
239	heimbach	1.2
240			read(cunit) filencvarindex(ivartype)
241			if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
242			& then
243			print *, 'ctrl_set_unpack_yz:WARNING: wrong ncvarindex ',
244			& filencvarindex(ivartype), ncvarindex(ivartype)
245			STOP 'in S/R ctrl_unpack'
246			endif
247			read(cunit) filej
248			read(cunit) filei
249			do k = 1, Nr
250	heimbach	1.3	cbuffindex = nwetglobal(k,iobcs)
251	heimbach	1.2	if ( cbuffindex .gt. 0 ) then
252			read(cunit) filencbuffindex
253			if (filencbuffindex .NE. cbuffindex) then
254			print *, 'WARNING: wrong cbuffindex ',
255			& filencbuffindex, cbuffindex
256			STOP 'in S/R ctrl_unpack'
257			endif
258			read(cunit) filek
259			if (filek .NE. k) then
260			print *, 'WARNING: wrong k ',
261			& filek, k
262			STOP 'in S/R ctrl_unpack'
263			endif
264			read(cunit) (cbuff(ii), ii=1,cbuffindex)
265			endif
266			cbuffindex = 0
267			do jp = 1,nPy
268			do bj = jtlo,jthi
269			do j = jmin,jmax
270			do ip = 1,nPx
271			do bi = itlo,ithi
272	heimbach	1.3	if ( globmskyz(bi,ip,j,bj,jp,k,iobcs) .ne. 0. ) then
273	heimbach	1.2	cbuffindex = cbuffindex + 1
274			globfldyz(bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
275			#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
276			globfldyz(bi,ip,j,bj,jp,k) =
277			& globfldyz(bi,ip,j,bj,jp,k)/
278	heimbach	1.3	# ifdef CTRL_UNPACK_PRECISE
279			& sqrt(weightfldyz(bi,ip,j,bj,jp,k,iobcs))
280			# else
281	heimbach	1.2	& sqrt(weightfld(k,iobcs))
282	heimbach	1.3	# endif
283			#endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
284	heimbach	1.2	else
285			globfldyz(bi,ip,j,bj,jp,k) = 0. _d 0
286			endif
287			enddo
288			enddo
289			enddo
290			enddo
291			enddo
292			c
293	heimbach	1.3	c -- end of k loop
294	heimbach	1.2	enddo
295
296			call MDSWRITEFIELD_YZ_GL( fname, ctrlprec, 'RL',
297			& Nr, globfldyz, irec,
298			& optimcycle, mythid)
299
300			c -- end of iobcs loop -- This loop has been removed.
301			cgg enddo
302			c -- end of irec loop --
303			enddo
304
305			_END_MASTER( mythid )
306
307			return
308			end
309