/[MITgcm]/MITgcm/pkg/ctrl/ctrl_set_pack_xyz.F
ViewVC logotype

Annotation of /MITgcm/pkg/ctrl/ctrl_set_pack_xyz.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.21 - (hide annotations) (download)
Thu Dec 23 02:43:43 2010 UTC (13 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62w, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint63
Changes since 1.20: +9 -9 lines
change arg. list of S/R MDSIO_PASS_R4/8toRL/S.

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

  ViewVC Help
Powered by ViewVC 1.1.22