/[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.28 - (hide annotations) (download)
Tue Jun 16 20:20:26 2015 UTC (8 years, 11 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65n, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, HEAD
Changes since 1.27: +15 -7 lines
- ctrl_map_ini_genarr.F: rename xx_etan0 as xx_etan, consistent with
  the other initial condition controls; add numbers of characters.
- ctrl_set_pack_xy.F, etc.: replace tests for .NOT.ctrlSmoothCorrel2D
  with tests for doPackOld = (.NOT.ctrlSmoothCorrel2D).AND.(.NOT.ctrlUseGen)
  etc. to disentangle packing options from ALLOW_SMOOTH_CORREL2D/3D CPP options.

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

  ViewVC Help
Powered by ViewVC 1.1.22