/[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.22 - (hide annotations) (download)
Tue Jul 19 12:44:36 2011 UTC (12 years, 10 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint63g, checkpoint63p, checkpoint63q, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c
Changes since 1.21: +3 -4 lines
remove unused variables and some obsolete code

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

  ViewVC Help
Powered by ViewVC 1.1.22