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

Contents of /MITgcm/pkg/ctrl/ctrl_set_unpack_xyz.F

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


Revision 1.32 - (show 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.31: +14 -6 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 C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_set_unpack_xyz.F,v 1.31 2015/03/23 21:07:37 gforget Exp $
2 C $Name: $
3
4 #include "CTRL_OPTIONS.h"
5
6 subroutine ctrl_set_unpack_xyz(
7 & lxxadxx, cunit, ivartype, fname, masktype, weighttype,
8 & weightfld, nwetglobal, mythid)
9
10 c ==================================================================
11 c SUBROUTINE ctrl_set_unpack_xyz
12 c ==================================================================
13 c
14 c o Unpack the control vector such that land points are filled in.
15 c
16 c o Use a more precise nondimensionalization that depends on (x,y)
17 c Added weighttype to the argument list so that I can geographically
18 c vary the nondimensionalization.
19 c gebbie@mit.edu, 18-Mar-2003
20 c
21 c ==================================================================
22
23 implicit none
24
25 c == global variables ==
26
27 #include "EEPARAMS.h"
28 #include "SIZE.h"
29 #include "PARAMS.h"
30 #include "GRID.h"
31
32 #include "ctrl.h"
33 #include "optim.h"
34
35 c == routine arguments ==
36
37 logical lxxadxx
38 integer cunit
39 integer ivartype
40 character*( 80) fname
41 character*( 9) masktype
42 character*( 80) weighttype
43 _RL weightfld( nr,nsx,nsy )
44 integer nwetglobal(nr)
45 integer mythid
46
47 #ifndef EXCLUDE_CTRL_PACK
48 # ifndef ALLOW_PACKUNPACK_METHOD2
49 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 #ifdef CTRL_UNPACK_PRECISE
66 integer il
67 character*(80) weightname
68 _RL weightfld3d( snx,nsx,npx,sny,nsy,npy,nr )
69 #endif
70 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
74 character*(128) cfile
75
76 _RL delZnorm
77 integer reclen, irectrue
78 integer cunit2, cunit3
79 character*(80) cfile2, cfile3
80
81 LOGICAL doPackOld
82
83 c == external ==
84
85 integer ilnblnk
86 external ilnblnk
87
88 cc == end of interface ==
89
90 jtlo = 1
91 jthi = nsy
92 itlo = 1
93 ithi = nsx
94 jmin = 1
95 jmax = sny
96 imin = 1
97 imax = snx
98
99 #ifdef CTRL_DELZNORM
100 delZnorm = 0.
101 do k = 1, Nr
102 delZnorm = delZnorm + delR(k)/FLOAT(Nr)
103 enddo
104 #endif
105
106 doPackOld = (.NOT.ctrlSmoothCorrel3D).AND.(.NOT.ctrlUseGen)
107
108 c Initialise temporary file
109 do k = 1,nr
110 do jp = 1,nPy
111 do bj = jtlo,jthi
112 do j = jmin,jmax
113 do ip = 1,nPx
114 do bi = itlo,ithi
115 do i = imin,imax
116 globfld3d (i,bi,ip,j,bj,jp,k) = 0. _d 0
117 globmsk (i,bi,ip,j,bj,jp,k) = 0. _d 0
118 globfldtmp2(i,bi,ip,j,bj,jp) = 0. _d 0
119 globfldtmp3(i,bi,ip,j,bj,jp) = 0. _d 0
120 enddo
121 enddo
122 enddo
123 enddo
124 enddo
125 enddo
126 enddo
127
128 c-- Only the master thread will do I/O.
129 _BEGIN_MASTER( mythid )
130
131 #ifdef CTRL_DELZNORM
132 do k = 1, nr
133 print *, 'ph-delznorm ', k, delZnorm, delR(k)
134 print *, 'ph-weight ', weightfld(k,1,1)
135 enddo
136 #endif
137
138 if ( doPackDiag ) then
139 write(cfile2(1:80),'(80a)') ' '
140 write(cfile3(1:80),'(80a)') ' '
141 if ( lxxadxx ) then
142 write(cfile2(1:80),'(a,I3.3,a,I4.4,a)')
143 & 'diag_unpack_nondim_ctrl_',
144 & ivartype, '_', optimcycle, '.bin'
145 write(cfile3(1:80),'(a,I3.3,a,I4.4,a)')
146 & 'diag_unpack_dimens_ctrl_',
147 & ivartype, '_', optimcycle, '.bin'
148 else
149 write(cfile2(1:80),'(a,I3.3,a,I4.4,a)')
150 & 'diag_unpack_nondim_grad_',
151 & ivartype, '_', optimcycle, '.bin'
152 write(cfile3(1:80),'(a,I3.3,a,I4.4,a)')
153 & 'diag_unpack_dimens_grad_',
154 & ivartype, '_', optimcycle, '.bin'
155 endif
156
157 reclen = FLOAT(snx*nsx*npx*sny*nsy*npy*4)
158 call mdsfindunit( cunit2, mythid )
159 open( cunit2, file=cfile2, status='unknown',
160 & access='direct', recl=reclen )
161 call mdsfindunit( cunit3, mythid )
162 open( cunit3, file=cfile3, status='unknown',
163 & access='direct', recl=reclen )
164 endif
165
166 #ifdef CTRL_UNPACK_PRECISE
167 if (weighttype.NE.' ') then
168 il=ilnblnk( weighttype)
169 write(weightname(1:80),'(80a)') ' '
170 write(weightname(1:80),'(a)') weighttype(1:il)
171 call MDSREADFIELD_3D_GL(
172 & weightname, ctrlprec, 'RL',
173 & Nr, weightfld3d, 1, mythid)
174 else
175 do k = 1,nr
176 do jp = 1,nPy
177 do bj = jtlo,jthi
178 do j = jmin,jmax
179 do ip = 1,nPx
180 do bi = itlo,ithi
181 do i = imin,imax
182 weightfld3d(i,bi,ip,j,bj,jp,k) = 1. _d 0
183 enddo
184 enddo
185 enddo
186 enddo
187 enddo
188 enddo
189 enddo
190 endif
191 #endif
192
193 call MDSREADFIELD_3D_GL(
194 & masktype, ctrlprec, 'RL',
195 & Nr, globmsk, 1, mythid)
196
197 do irec = 1, ncvarrecs(ivartype)
198 #ifndef ALLOW_ADMTLM
199 read(cunit) filencvarindex(ivartype)
200 if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
201 & then
202 print *, 'ctrl_set_unpack_xyz:WARNING: wrong ncvarindex ',
203 & filencvarindex(ivartype), ncvarindex(ivartype)
204 STOP 'in S/R ctrl_set_unpack_xyz'
205 endif
206 read(cunit) filej
207 read(cunit) filei
208 #endif /* ALLOW_ADMTLM */
209 do k = 1, Nr
210 irectrue = (irec-1)*nr + k
211 if ( doZscaleUnpack ) then
212 delZnorm = (delR(1)/delR(k))**delZexp
213 else
214 delZnorm = 1. _d 0
215 endif
216 cbuffindex = nwetglobal(k)
217 if ( cbuffindex .gt. 0 ) then
218 #ifndef ALLOW_ADMTLM
219 read(cunit) filencbuffindex
220 if (filencbuffindex .NE. cbuffindex) then
221 print *, 'WARNING: wrong cbuffindex ',
222 & filencbuffindex, cbuffindex
223 STOP 'in S/R ctrl_set_unpack_xyz'
224 endif
225 read(cunit) filek
226 if (filek .NE. k) then
227 print *, 'WARNING: wrong k ',
228 & filek, k
229 STOP 'in S/R ctrl_set_unpack_xyz'
230 endif
231 cph#endif /* ALLOW_ADMTLM */
232 read(cunit) (cbuff(ii), ii=1,cbuffindex)
233 #endif /* ALLOW_ADMTLM */
234 endif
235 c
236 cbuffindex = 0
237 do jp = 1,nPy
238 do bj = jtlo,jthi
239 do j = jmin,jmax
240 do ip = 1,nPx
241 do bi = itlo,ithi
242 do i = imin,imax
243 if ( globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
244 cbuffindex = cbuffindex + 1
245 globfld3d(i,bi,ip,j,bj,jp,k) = cbuff(cbuffindex)
246 cph(
247 globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
248 cph)
249 #ifdef ALLOW_ADMTLM
250 nveccount = nveccount + 1
251 globfld3d(i,bi,ip,j,bj,jp,k) =
252 & phtmpadmtlm(nveccount)
253 cph(
254 globfldtmp2(i,bi,ip,j,bj,jp) =
255 & phtmpadmtlm(nveccount)
256 cph)
257 #endif
258 IF ( doPackOld ) THEN
259 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
260 if ( lxxadxx ) then
261 globfld3d(i,bi,ip,j,bj,jp,k) = delZnorm
262 & * globfld3d(i,bi,ip,j,bj,jp,k)
263 # ifdef CTRL_UNPACK_PRECISE
264 & / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
265 # else
266 & / sqrt(weightfld(k,bi,bj))
267 # endif
268 else
269 globfld3d(i,bi,ip,j,bj,jp,k) = 1/delZnorm
270 & * globfld3d(i,bi,ip,j,bj,jp,k)
271 # ifdef CTRL_UNPACK_PRECISE
272 & * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
273 # else
274 & * sqrt(weightfld(k,bi,bj))
275 # endif
276 endif
277 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
278 ENDIF !IF ( doPackOld ) THEN
279 else
280 globfld3d(i,bi,ip,j,bj,jp,k) = 0. _d 0
281 endif
282 cph(
283 globfldtmp3(i,bi,ip,j,bj,jp) =
284 & globfld3d(i,bi,ip,j,bj,jp,k)
285 cph)
286 enddo
287 enddo
288 enddo
289 enddo
290 enddo
291 enddo
292 c
293 if ( doPackDiag ) then
294 write(cunit2,rec=irectrue) globfldtmp2
295 write(cunit3,rec=irectrue) globfldtmp3
296 endif
297 c
298 enddo
299
300 call MDSWRITEFIELD_3D_GL( fname, ctrlprec, 'RL',
301 & Nr, globfld3d,
302 & irec, optimcycle, mythid)
303
304 enddo
305
306 if ( doPackDiag ) then
307 close ( cunit2 )
308 close ( cunit3 )
309 endif
310
311 _END_MASTER( mythid )
312
313 # else
314 c == local variables ==
315
316 integer bi,bj
317 integer ip,jp
318 integer i,j,k
319 integer ii
320 integer il
321 integer irec
322 integer itlo,ithi
323 integer jtlo,jthi
324
325 integer cbuffindex
326
327 _RL msk3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
328 real*8 msk2d_buf(sNx,sNy,nSx,nSy)
329 real*8 msk2d_buf_glo(Nx,Ny)
330
331 _RL fld3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
332 real*8 fld2d_buf(sNx,sNy,nSx,nSy)
333 real*8 fld2d_buf_glo(Nx,Ny)
334
335 _RL fld3dDim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
336 _RL fld3dNodim(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
337
338 #ifdef CTRL_UNPACK_PRECISE
339 _RL wei3d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
340 #endif
341
342 real*4 cbuff ( snx*nsx*npx*sny*nsy*npy )
343
344 character*(80) weightname
345 _RL delZnorm
346 character*(80) cfile2, cfile3
347 _RL dummy
348
349 LOGICAL doPackOld
350
351 c == external ==
352
353 integer ilnblnk
354 external ilnblnk
355
356 c == end of interface ==
357
358 c-- part 1: preliminary reads and definitions
359
360 doPackOld = (.NOT.ctrlSmoothCorrel3D).AND.(.NOT.ctrlUseGen)
361
362 #ifdef CTRL_UNPACK_PRECISE
363 #ifdef ALLOW_AUTODIFF
364 call active_read_xyz(weighttype, wei3d, 1,
365 & .FALSE., .FALSE., 0 , mythid, dummy)
366 #else
367 CALL READ_REC_XYZ_RL( weighttype, wei3d, 1, 1, myThid )
368 #endif
369 #endif
370
371 #ifdef ALLOW_AUTODIFF
372 call active_read_xyz(masktype, msk3d, 1,
373 & .FALSE., .FALSE., 0 , mythid, dummy)
374 #else
375 CALL READ_REC_XYZ_RL( masktype, msk3d, 1, 1, myThid )
376 #endif
377
378 if ( doPackDiag ) then
379 write(cfile2(1:80),'(80a)') ' '
380 write(cfile3(1:80),'(80a)') ' '
381 il = ilnblnk( fname )
382 if ( lxxadxx ) then
383 write(cfile2(1:80),'(2a)') fname(1:il),'.unpack_ctrl_adim'
384 write(cfile3(1:80),'(2a)') fname(1:il),'.unpack_ctrl_dim'
385 else
386 write(cfile2(1:80),'(2a)') fname(1:il),'.unpack_grad_adim'
387 write(cfile3(1:80),'(2a)') fname(1:il),'.unpack_grad_dim'
388 endif
389 endif
390
391 c-- part 2: loop over records
392
393 do irec = 1, ncvarrecs(ivartype)
394
395 c-- 2.1: array <- buffer <- global buffer <- global file
396
397 #ifndef ALLOW_ADMTLM
398 _BEGIN_MASTER( mythid )
399 IF ( myProcId .eq. 0 ) THEN
400 read(cunit) filencvarindex(ivartype)
401 if (filencvarindex(ivartype) .NE. ncvarindex(ivartype))
402 & then
403 print *, 'ctrl_set_unpack_xyz:WARNING: wrong ncvarindex ',
404 & filencvarindex(ivartype), ncvarindex(ivartype)
405 STOP 'in S/R ctrl_set_unpack_xyz'
406 endif
407 read(cunit) filej
408 read(cunit) filei
409 ENDIF
410 _END_MASTER( mythid )
411 _BARRIER
412 #endif /* ALLOW_ADMTLM */
413
414 do k = 1, nr
415
416 CALL MDS_PASS_R8toRL( msk2d_buf, msk3d,
417 & 0, 0, 1, k, Nr, 0, 0, .FALSE., myThid )
418 CALL BAR2( myThid )
419 CALL GATHER_2D_R8( msk2d_buf_glo, msk2d_buf,
420 & Nx,Ny,.FALSE.,.TRUE.,myThid)
421 CALL BAR2( myThid )
422
423 _BEGIN_MASTER( mythid )
424 cbuffindex = nwetglobal(k)
425 IF ( myProcId .eq. 0 ) THEN
426
427 #ifndef ALLOW_ADMTLM
428 if ( cbuffindex .gt. 0) then
429 read(cunit) filencbuffindex
430 read(cunit) filek
431 if (filencbuffindex .NE. cbuffindex) then
432 print *, 'WARNING: wrong cbuffindex ',
433 & filencbuffindex, cbuffindex
434 STOP 'in S/R ctrl_set_unpack_xyz'
435 endif
436 if (filek .NE. k) then
437 print *, 'WARNING: wrong k ',
438 & filek, k
439 STOP 'in S/R ctrl_set_unpack_xyz'
440 endif
441 read(cunit) (cbuff(ii), ii=1,cbuffindex)
442 endif
443 #endif
444
445 cbuffindex = 0
446 DO j=1,Ny
447 DO i=1,Nx
448 if (msk2d_buf_glo(i,j) .ne. 0. ) then
449 cbuffindex = cbuffindex + 1
450 fld2d_buf_glo(i,j) = cbuff(cbuffindex)
451 endif
452 ENDDO
453 ENDDO
454
455 ENDIF
456 _END_MASTER( mythid )
457 _BARRIER
458
459 CALL BAR2( myThid )
460 CALL SCATTER_2D_R8( fld2d_buf_glo, fld2d_buf,
461 & Nx,Ny,.FALSE.,.TRUE.,myThid)
462 CALL BAR2( myThid )
463 CALL MDS_PASS_R8toRL( fld2d_buf, fld3dNodim,
464 & 0, 0, 1, k, Nr, 0, 0, .TRUE., myThid )
465
466 enddo !do k = 1, nr
467
468 c-- 2.2: normalize field if needed
469 DO bj = myByLo(myThid), myByHi(myThid)
470 DO bi = myBxLo(myThid), myBxHi(myThid)
471 DO k=1,Nr
472 if ( doZscalePack ) then
473 delZnorm = (delR(1)/delR(k))**delZexp
474 else
475 delZnorm = 1. _d 0
476 endif
477 DO j=1,sNy
478 DO i=1,sNx
479 if (msk3d(i,j,k,bi,bj).EQ.0. _d 0) then
480 fld3dDim(i,j,k,bi,bj)=0. _d 0
481 fld3dNodim(i,j,k,bi,bj)=0. _d 0
482 else
483 #ifdef ALLOW_ADMTLM
484 nveccount = nveccount + 1
485 fld3dNodim(i,j,k,bi,bj)=phtmpadmtlm(nveccount)
486 #endif
487 IF ( .NOT.doPackOld ) THEN
488 fld3dDim(i,j,k,bi,bj)=fld3dNodim(i,j,k,bi,bj)
489 ELSE !IF ( .NOT.doPackOld ) THEN
490 # ifndef ALLOW_NONDIMENSIONAL_CONTROL_IO
491 fld3dDim(i,j,k,bi,bj) = fld3dNodim(i,j,k,bi,bj)
492 # else
493 if (lxxadxx) then
494 fld3dDim(i,j,k,bi,bj) =
495 & fld3dNodim(i,j,k,bi,bj) * delZnorm
496 # ifdef CTRL_UNPACK_PRECISE
497 & / sqrt(wei3d(i,j,k,bi,bj))
498 # else
499 & / sqrt(weightfld(k,bi,bj))
500 # endif
501 else
502 fld3dDim(i,j,k,bi,bj) =
503 & fld3dNodim(i,j,k,bi,bj) / delZnorm
504 # ifdef CTRL_UNPACK_PRECISE
505 & * sqrt(wei3d(i,j,k,bi,bj))
506 # else
507 & * sqrt(weightfld(k,bi,bj))
508 # endif
509 endif
510 # endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
511 ENDIF !IF ( .NOT.doPackOld ) THEN
512 endif
513 ENDDO
514 ENDDO
515 ENDDO
516 ENDDO
517 ENDDO
518
519 c-- 2.3:
520 if ( doPackDiag ) then
521 c error: twice the same one
522 call WRITE_REC_3D_RL( cfile2, ctrlprec,
523 & Nr, fld3dNodim, irec, 0, mythid)
524 call WRITE_REC_3D_RL( cfile3, ctrlprec,
525 & Nr, fld3dDim, irec, 0, mythid)
526 endif
527
528 c-- 2.4:
529 call WRITE_REC_3D_RL( fname, ctrlprec,
530 & Nr, fld3dDim, irec, 0, mythid)
531
532
533 enddo !do irec = 1, ncvarrecs(ivartype)
534
535 # endif /* ALLOW_PACKUNPACK_METHOD2 */
536 # endif /* EXCLUDE_CTRL_PACK */
537
538 return
539 end
540

  ViewVC Help
Powered by ViewVC 1.1.22