pkg/ctrl/ctrl_set_unpack_xy.F


#include "CTRL_CPPOPTIONS.h"

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

c     ==================================================================
c     SUBROUTINE ctrl_set_unpack_xy
c     ==================================================================
c
c     o Unpack the control vector such that the land points are filled
c       in.
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 ==

      logical lxxadxx
      integer cunit
      integer ivartype
      character*( 80)   fname, fnameGlobal
      character*(  9)   masktype
      character*( 80) weighttype
      integer nwetglobal(nr)
      integer mythid

c     == local variables ==

      integer bi,bj
      integer ip,jp
      integer i,j,k
      integer ii
      integer il
      integer irec,nrec_nl
      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 )
      _RL     globfld2d( snx,nsx,npx,sny,nsy,npy )
      real*4  cbuff    ( snx*nsx*npx*sny*nsy*npy )

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

      integer reclen,irectrue
      integer cunit2, cunit3
      character*(80) cfile2, cfile3
      real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
      real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )

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

      nbuffGlobal = nbuffGlobal + 1

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_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 ALLOW_NONDIMENSIONAL_CONTROL_IO
      il=ilnblnk( weighttype)
      write(weightname(1:80),'(80a)') ' '
      write(weightname(1:80),'(a)') weighttype(1:il)
      call MDSREADFIELD_2D_GL( 
     &     weightname, ctrlprec, 'RL',
     &     1, globfld2d, 1, mythid)
#endif
      
      call MDSREADFIELD_3D_GL( 
     &     masktype, ctrlprec, 'RL',
     &     Nr, globmsk, 1, mythid)

      nrec_nl=int(ncvarrecs(ivartype)/Nr)
      do irec = 1, nrec_nl
         print *, 'ph-pack nrec_nl = ', irec, nrec_nl, ivartype,
     &          ncvarrecs(ivartype)
         do k = 1,Nr
         irectrue = (irec-1)*nr + k
         read(cunit) filencvarindex(ivartype)
         if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
     &        then
            print *, 'ctrl_set_unpack_xy:WARNING: wrong ncvarindex ',
     &           filencvarindex(ivartype), ncvarindex(ivartype)
            STOP 'in S/R ctrl_unpack'
         endif
         read(cunit) filej
         read(cunit) filei
            cbuffindex = nwetglobal(1)
            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. 1) then
                  print *, 'WARNING: wrong k ',
     &                 filek, 1
                  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
                 do i = imin,imax
                  if ( globmsk(i,bi,ip,j,bj,jp,1) .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_NONDIMENSIONAL_CONTROL_IO
                     if ( lxxadxx ) then
                        globfld3d(i,bi,ip,j,bj,jp,k) =
     &                       globfld3d(i,bi,ip,j,bj,jp,k)/
     &                       sqrt(globfld2d(i,bi,ip,j,bj,jp))
                     else
                        globfld3d(i,bi,ip,j,bj,jp,k) =
     &                       globfld3d(i,bi,ip,j,bj,jp,k)*
     &                       sqrt(globfld2d(i,bi,ip,j,bj,jp))
                     endif
#endif
                  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
cph(
            if ( doPackDiag ) then
               write(cunit2,rec=irectrue) globfldtmp2
               write(cunit3,rec=irectrue) globfldtmp3
            endif
cph)
         enddo

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

      enddo

      do irec = nrec_nl*Nr+1,ncvarrecs(ivartype)
         print  *, 'ph-pack nrec_nl+irec ', irec, nrec_nl, ivartype, 
     &             ncvarrecs(ivartype) 
#ifndef ALLOW_ADMTLM
         read(cunit) filencvarindex(ivartype)
         if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
     &        then
            print *, 'ctrl_set_unpack_xy:WARNING: wrong ncvarindex ',
     &           filencvarindex(ivartype), ncvarindex(ivartype)
            STOP 'in S/R ctrl_unpack'
         endif
         read(cunit) filej
         read(cunit) filei
#endif /* ndef ALLOW_ADMTLM */
         do k = 1,1
            irectrue = irec
            cbuffindex = nwetglobal(k)
#ifndef ALLOW_ADMTLM
            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
#else ALLOW_ADMTLM
            write(fnameGlobal(1:80),'(a)') ' '
            write(fnameGlobal,'(a,i4.4)') 
     &           'admtlm_vector.it', optimcycle
            call mdsreadvector( fnameGlobal, 64, 'RL',
     &           admtlmrec, cbuffGlobal, 1, 1, nbuffGlobal, mythid )
            do ii = 1, cbuffindex
               cbuff(ii) = cbuffGlobal(ii)
            enddo
#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
                     globfld3d(i,bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
cph(
                     globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
cph)
#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
                     if ( lxxadxx ) then
                        globfld3d(i,bi,ip,j,bj,jp,k) = 
     &                       globfld3d(i,bi,ip,j,bj,jp,k)/
     &                       sqrt(globfld2d(i,bi,ip,j,bj,jp))
                     else
                        globfld3d(i,bi,ip,j,bj,jp,k) = 
     &                       globfld3d(i,bi,ip,j,bj,jp,k)*
     &                       sqrt(globfld2d(i,bi,ip,j,bj,jp))
                     endif
#endif
                  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
cph(
            if ( doPackDiag ) then
               write(cunit2,rec=irectrue) globfldtmp2
               write(cunit3,rec=irectrue) globfldtmp3
            endif
cph)
         enddo
             
         call MDSWRITEFIELD_2D_GL( fname, ctrlprec, 'RL',
     &                             1, globfld3d(1,1,1,1,1,1,1),
     &                             irec,  optimcycle, mythid)

      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_unpack_xy(
5	heimbach	1.11	& lxxadxx, cunit, ivartype, fname, masktype, weighttype,
6	heimbach	1.2	& nwetglobal, mythid)
7
8			c ==================================================================
9			c SUBROUTINE ctrl_set_unpack_xy
10			c ==================================================================
11			c
12			c o Unpack the control vector such that the land points are filled
13			c in.
14			c
15	heimbach	1.3	c changed: heimbach@mit.edu 17-Jun-2003
16			c merged Armin's changes to replace write of
17			c nr * globfld2d by 1 * globfld3d
18			c (ad hoc fix to speed up global I/O)
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	heimbach	1.11	logical lxxadxx
37	heimbach	1.2	integer cunit
38			integer ivartype
39	heimbach	1.9	character*( 80) fname, fnameGlobal
40	heimbach	1.14	character*( 9) masktype
41	heimbach	1.2	character*( 80) weighttype
42			integer nwetglobal(nr)
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	heimbach	1.3	integer irec,nrec_nl
53	heimbach	1.2	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			_RL globfld2d( snx,nsx,npx,sny,nsy,npy )
63	heimbach	1.11	real4 cbuff ( snxnsxnpxsnynsynpy )
64	heimbach	1.2
65			character*(128) cfile
66			character*( 80) weightname
67
68	heimbach	1.11	integer reclen,irectrue
69			integer cunit2, cunit3
70			character*(80) cfile2, cfile3
71	heimbach	1.12	real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
72			real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )
73	heimbach	1.11
74	heimbach	1.2	c == external ==
75
76			integer ilnblnk
77			external ilnblnk
78
79			c == end of interface ==
80
81			jtlo = 1
82			jthi = nsy
83			itlo = 1
84			ithi = nsx
85			jmin = 1
86			jmax = sny
87			imin = 1
88			imax = snx
89
90	heimbach	1.9	nbuffGlobal = nbuffGlobal + 1
91
92	heimbach	1.2	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			do i = imin,imax
100	heimbach	1.12	globfld3d (i,bi,ip,j,bj,jp,k) = 0. _d 0
101			globmsk (i,bi,ip,j,bj,jp,k) = 0. _d 0
102			globfldtmp2(i,bi,ip,j,bj,jp) = 0.
103			globfldtmp3(i,bi,ip,j,bj,jp) = 0.
104	heimbach	1.2	enddo
105			enddo
106			enddo
107			enddo
108			enddo
109			enddo
110			enddo
111
112			c-- Only the master thread will do I/O.
113			_BEGIN_MASTER( mythid )
114
115	heimbach	1.11	if ( doPackDiag ) then
116			write(cfile2(1:80),'(80a)') ' '
117			write(cfile3(1:80),'(80a)') ' '
118			if ( lxxadxx ) then
119			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
120			& 'diag_unpack_nondim_ctrl_',
121			& ivartype, '_', optimcycle, '.bin'
122			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
123			& 'diag_unpack_dimens_ctrl_',
124			& ivartype, '_', optimcycle, '.bin'
125			else
126			write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
127			& 'diag_unpack_nondim_grad_',
128			& ivartype, '_', optimcycle, '.bin'
129			write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
130			& 'diag_unpack_dimens_grad_',
131			& ivartype, '_', optimcycle, '.bin'
132			endif
133
134	heimbach	1.12	reclen = FLOAT(snxnsxnpxsnynsynpy4)
135	heimbach	1.11	call mdsfindunit( cunit2, mythid )
136	heimbach	1.13	open( cunit2, file=cfile2, status='unknown',
137	heimbach	1.11	& access='direct', recl=reclen )
138			call mdsfindunit( cunit3, mythid )
139	heimbach	1.13	open( cunit3, file=cfile3, status='unknown',
140	heimbach	1.11	& access='direct', recl=reclen )
141			endif
142
143	heimbach	1.2	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
144			il=ilnblnk( weighttype)
145			write(weightname(1:80),'(80a)') ' '
146			write(weightname(1:80),'(a)') weighttype(1:il)
147			call MDSREADFIELD_2D_GL(
148			& weightname, ctrlprec, 'RL',
149			& 1, globfld2d, 1, mythid)
150			#endif
151
152			call MDSREADFIELD_3D_GL(
153			& masktype, ctrlprec, 'RL',
154			& Nr, globmsk, 1, mythid)
155
156	heimbach	1.3	nrec_nl=int(ncvarrecs(ivartype)/Nr)
157			do irec = 1, nrec_nl
158	heimbach	1.7	print *, 'ph-pack nrec_nl = ', irec, nrec_nl, ivartype,
159			& ncvarrecs(ivartype)
160	heimbach	1.3	do k = 1,Nr
161	heimbach	1.11	irectrue = (irec-1)*nr + k
162	heimbach	1.3	read(cunit) filencvarindex(ivartype)
163			if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
164			& then
165			print *, 'ctrl_set_unpack_xy:WARNING: wrong ncvarindex ',
166			& filencvarindex(ivartype), ncvarindex(ivartype)
167			STOP 'in S/R ctrl_unpack'
168			endif
169			read(cunit) filej
170			read(cunit) filei
171			cbuffindex = nwetglobal(1)
172			if ( cbuffindex .gt. 0 ) then
173			read(cunit) filencbuffindex
174			if (filencbuffindex .NE. cbuffindex) then
175			print *, 'WARNING: wrong cbuffindex ',
176			& filencbuffindex, cbuffindex
177			STOP 'in S/R ctrl_unpack'
178			endif
179			read(cunit) filek
180			if (filek .NE. 1) then
181			print *, 'WARNING: wrong k ',
182			& filek, 1
183			STOP 'in S/R ctrl_unpack'
184			endif
185			read(cunit) (cbuff(ii), ii=1,cbuffindex)
186			endif
187			cbuffindex = 0
188			do jp = 1,nPy
189			do bj = jtlo,jthi
190			do j = jmin,jmax
191			do ip = 1,nPx
192			do bi = itlo,ithi
193			do i = imin,imax
194			if ( globmsk(i,bi,ip,j,bj,jp,1) .ne. 0. ) then
195			cbuffindex = cbuffindex + 1
196			globfld3d(i,bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
197	heimbach	1.12	cph(
198			globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
199			cph)
200	heimbach	1.3	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
201	heimbach	1.11	if ( lxxadxx ) then
202			globfld3d(i,bi,ip,j,bj,jp,k) =
203			& globfld3d(i,bi,ip,j,bj,jp,k)/
204			& sqrt(globfld2d(i,bi,ip,j,bj,jp))
205			else
206			globfld3d(i,bi,ip,j,bj,jp,k) =
207			& globfld3d(i,bi,ip,j,bj,jp,k)*
208			& sqrt(globfld2d(i,bi,ip,j,bj,jp))
209			endif
210	heimbach	1.3	#endif
211			else
212			globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
213	heimbach	1.12	endif
214			cph(
215			globfldtmp3(i,bi,ip,j,bj,jp) =
216			& globfld3d(i,bi,ip,j,bj,jp,k)
217			cph)
218			enddo
219	heimbach	1.3	enddo
220			enddo
221			enddo
222			enddo
223			enddo
224	heimbach	1.12	cph(
225			if ( doPackDiag ) then
226			write(cunit2,rec=irectrue) globfldtmp2
227			write(cunit3,rec=irectrue) globfldtmp3
228			endif
229			cph)
230	heimbach	1.3	enddo
231
232			call MDSWRITEFIELD_3D_GL( fname, ctrlprec, 'RL',
233			& NR, globfld3d,
234			& irec, optimcycle, mythid)
235
236			enddo
237
238			do irec = nrec_nl*Nr+1,ncvarrecs(ivartype)
239	heimbach	1.7	print *, 'ph-pack nrec_nl+irec ', irec, nrec_nl, ivartype,
240			& ncvarrecs(ivartype)
241	heimbach	1.10	#ifndef ALLOW_ADMTLM
242	heimbach	1.2	read(cunit) filencvarindex(ivartype)
243			if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
244			& then
245			print *, 'ctrl_set_unpack_xy:WARNING: wrong ncvarindex ',
246			& filencvarindex(ivartype), ncvarindex(ivartype)
247			STOP 'in S/R ctrl_unpack'
248			endif
249			read(cunit) filej
250			read(cunit) filei
251	heimbach	1.10	#endif /* ndef ALLOW_ADMTLM */
252	heimbach	1.2	do k = 1,1
253	heimbach	1.11	irectrue = irec
254	heimbach	1.2	cbuffindex = nwetglobal(k)
255	heimbach	1.10	#ifndef ALLOW_ADMTLM
256	heimbach	1.2	if ( cbuffindex .gt. 0 ) then
257			read(cunit) filencbuffindex
258			if (filencbuffindex .NE. cbuffindex) then
259			print *, 'WARNING: wrong cbuffindex ',
260			& filencbuffindex, cbuffindex
261			STOP 'in S/R ctrl_unpack'
262			endif
263			read(cunit) filek
264			if (filek .NE. k) then
265			print *, 'WARNING: wrong k ',
266			& filek, k
267			STOP 'in S/R ctrl_unpack'
268			endif
269			read(cunit) (cbuff(ii), ii=1,cbuffindex)
270			endif
271	heimbach	1.10	#else ALLOW_ADMTLM
272	heimbach	1.9	write(fnameGlobal(1:80),'(a)') ' '
273			write(fnameGlobal,'(a,i4.4)')
274			& 'admtlm_vector.it', optimcycle
275			call mdsreadvector( fnameGlobal, 64, 'RL',
276	heimbach	1.10	& admtlmrec, cbuffGlobal, 1, 1, nbuffGlobal, mythid )
277			do ii = 1, cbuffindex
278			cbuff(ii) = cbuffGlobal(ii)
279			enddo
280	heimbach	1.9	#endif
281	heimbach	1.2	cbuffindex = 0
282			do jp = 1,nPy
283			do bj = jtlo,jthi
284			do j = jmin,jmax
285			do ip = 1,nPx
286			do bi = itlo,ithi
287			do i = imin,imax
288			if ( globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
289			cbuffindex = cbuffindex + 1
290			globfld3d(i,bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
291	heimbach	1.12	cph(
292			globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
293			cph)
294	heimbach	1.2	#ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
295	heimbach	1.11	if ( lxxadxx ) then
296			globfld3d(i,bi,ip,j,bj,jp,k) =
297			& globfld3d(i,bi,ip,j,bj,jp,k)/
298			& sqrt(globfld2d(i,bi,ip,j,bj,jp))
299			else
300			globfld3d(i,bi,ip,j,bj,jp,k) =
301			& globfld3d(i,bi,ip,j,bj,jp,k)*
302			& sqrt(globfld2d(i,bi,ip,j,bj,jp))
303			endif
304	heimbach	1.2	#endif
305			else
306			globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
307			endif
308	heimbach	1.12	cph(
309			globfldtmp3(i,bi,ip,j,bj,jp) =
310			& globfld3d(i,bi,ip,j,bj,jp,k)
311			cph)
312	heimbach	1.2	enddo
313			enddo
314			enddo
315			enddo
316			enddo
317			enddo
318	heimbach	1.12	cph(
319			if ( doPackDiag ) then
320			write(cunit2,rec=irectrue) globfldtmp2
321			write(cunit3,rec=irectrue) globfldtmp3
322			endif
323			cph)
324	heimbach	1.2	enddo
325
326			call MDSWRITEFIELD_2D_GL( fname, ctrlprec, 'RL',
327			& 1, globfld3d(1,1,1,1,1,1,1),
328			& irec, optimcycle, mythid)
329
330			enddo
331
332	heimbach	1.11	if ( doPackDiag ) then
333			close ( cunit2 )
334			close ( cunit3 )
335			endif
336
337	heimbach	1.2	_END_MASTER( mythid )
338
339			return
340			end
341