/[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.25 - (hide annotations) (download)
Fri Feb 1 19:25:32 2013 UTC (11 years, 3 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64l, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f
Changes since 1.24: +41 -24 lines
More additions to enable simple rescaling of generic control variables

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

  ViewVC Help
Powered by ViewVC 1.1.22