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

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

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


Revision 1.13 - (show annotations) (download)
Tue Nov 1 04:09:46 2005 UTC (18 years, 7 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58b_post, checkpoint57y_post, checkpoint58, checkpoint58d_post, checkpoint58a_post, checkpoint57z_post, checkpoint57y_pre, checkpoint57w_post, checkpoint57x_post, checkpoint58c_post
Changes since 1.12: +8 -0 lines
Completely restructured the arpack2model interface.
Now (again) only 1-d wetpoint vector is passed to ARPACK.
ctrl_unpack/pack are mimiced by admtlm_dsvd2model/model2dsvd

1
2 #include "CTRL_CPPOPTIONS.h"
3
4 subroutine ctrl_set_pack_xyz(
5 & cunit, ivartype, fname, masktype, weighttype,
6 & weightfld, lxxadxx, mythid)
7
8 c ==================================================================
9 c SUBROUTINE ctrl_set_pack_xyz
10 c ==================================================================
11 c
12 c o Compress the control vector such that only ocean points are
13 c written to file.
14 c
15 c o Use a more precise nondimensionalization that depends on (x,y)
16 c Added weighttype to the argument list so that I can geographically
17 c vary the nondimensionalization.
18 c gebbie@mit.edu, 18-Mar-2003
19 c
20 c ==================================================================
21
22 implicit none
23
24 c == global variables ==
25
26 #include "EEPARAMS.h"
27 #include "SIZE.h"
28 #include "PARAMS.h"
29 #include "GRID.h"
30
31 #include "ctrl.h"
32 #include "optim.h"
33
34 c == routine arguments ==
35
36 integer cunit
37 integer ivartype
38 character*( 80) fname
39 character*( 9) masktype
40 character*( 80) weighttype
41 _RL weightfld( nr,nsx,nsy )
42 logical lxxadxx
43 integer mythid
44
45 c == local variables ==
46
47 integer bi,bj
48 integer ip,jp
49 integer i,j,k
50 integer ii
51 integer il
52 integer irec
53 integer itlo,ithi
54 integer jtlo,jthi
55 integer jmin,jmax
56 integer imin,imax
57
58 integer cbuffindex
59
60 _RL globmsk ( snx,nsx,npx,sny,nsy,npy,nr )
61 _RL globfld3d( snx,nsx,npx,sny,nsy,npy,nr )
62 #ifdef CTRL_PACK_PRECISE
63 _RL weightfld3d( snx,nsx,npx,sny,nsy,npy,nr )
64 #endif
65 real*4 cbuff ( snx*nsx*npx*sny*nsy*npy )
66 real*4 globfldtmp2( snx,nsx,npx,sny,nsy,npy )
67 real*4 globfldtmp3( snx,nsx,npx,sny,nsy,npy )
68
69 character*(80) weightname
70
71 _RL delZnorm
72 integer reclen, irectrue
73 integer cunit2, cunit3
74 character*(80) cfile2, cfile3
75
76 c == external ==
77
78 integer ilnblnk
79 external ilnblnk
80
81 c == end of interface ==
82
83 jtlo = 1
84 jthi = nsy
85 itlo = 1
86 ithi = nsx
87 jmin = 1
88 jmax = sny
89 imin = 1
90 imax = snx
91
92 #ifdef CTRL_DELZNORM
93 delZnorm = 0.
94 do k = 1, Nr
95 delZnorm = delZnorm + delR(k)/FLOAT(Nr)
96 enddo
97 #endif
98
99 c Initialise temporary file
100 do k = 1,nr
101 do jp = 1,nPy
102 do bj = jtlo,jthi
103 do j = jmin,jmax
104 do ip = 1,nPx
105 do bi = itlo,ithi
106 do i = imin,imax
107 globfld3d (i,bi,ip,j,bj,jp,k) = 0. _d 0
108 globmsk (i,bi,ip,j,bj,jp,k) = 0. _d 0
109 globfldtmp2(i,bi,ip,j,bj,jp) = 0.
110 globfldtmp3(i,bi,ip,j,bj,jp) = 0.
111 enddo
112 enddo
113 enddo
114 enddo
115 enddo
116 enddo
117 enddo
118
119 c-- Only the master thread will do I/O.
120 _BEGIN_MASTER( mythid )
121
122 if ( doPackDiag ) then
123 write(cfile2(1:80),'(80a)') ' '
124 write(cfile3(1:80),'(80a)') ' '
125 if ( lxxadxx ) then
126 write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
127 & 'diag_pack_nonout_ctrl_',
128 & ivartype, '_', optimcycle, '.bin'
129 write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
130 & 'diag_pack_dimout_ctrl_',
131 & ivartype, '_', optimcycle, '.bin'
132 else
133 write(cfile2(1:80),'(a,I2.2,a,I4.4,a)')
134 & 'diag_pack_nonout_grad_',
135 & ivartype, '_', optimcycle, '.bin'
136 write(cfile3(1:80),'(a,I2.2,a,I4.4,a)')
137 & 'diag_pack_dimout_grad_',
138 & ivartype, '_', optimcycle, '.bin'
139 endif
140
141 reclen = FLOAT(snx*nsx*npx*sny*nsy*npy*4)
142 call mdsfindunit( cunit2, mythid )
143 open( cunit2, file=cfile2, status='unknown',
144 & access='direct', recl=reclen )
145 call mdsfindunit( cunit3, mythid )
146 open( cunit3, file=cfile3, status='unknown',
147 & access='direct', recl=reclen )
148 endif
149
150 #ifdef CTRL_PACK_PRECISE
151 il=ilnblnk( weighttype)
152 write(weightname(1:80),'(80a)') ' '
153 write(weightname(1:80),'(a)') weighttype(1:il)
154
155 call MDSREADFIELD_3D_GL(
156 & weightname, ctrlprec, 'RL',
157 & Nr, weightfld3d, 1, mythid)
158 #endif
159
160 call MDSREADFIELD_3D_GL(
161 & masktype, ctrlprec, 'RL',
162 & Nr, globmsk, 1, mythid)
163
164 do irec = 1, ncvarrecs(ivartype)
165
166 call MDSREADFIELD_3D_GL( fname, ctrlprec, 'RL',
167 & Nr, globfld3d, irec, mythid)
168
169 #ifndef ALLOW_ADMTLM
170 write(cunit) ncvarindex(ivartype)
171 write(cunit) 1
172 write(cunit) 1
173 #endif
174 do k = 1, nr
175 irectrue = (irec-1)*nr + k
176 if ( doZscalePack ) then
177 cph delZnorm = SQRT(delR(1)/delR(k))
178 delZnorm = delR(1)/delR(k)
179 else
180 delZnorm = 1. _d 0
181 endif
182 cbuffindex = 0
183 do jp = 1,nPy
184 do bj = jtlo,jthi
185 do j = jmin,jmax
186 do ip = 1,nPx
187 do bi = itlo,ithi
188 do i = imin,imax
189 if (globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
190 cbuffindex = cbuffindex + 1
191 cph(
192 globfldtmp3(i,bi,ip,j,bj,jp) =
193 & globfld3d(i,bi,ip,j,bj,jp,k)
194 cph)
195 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
196 if (lxxadxx) then
197 cbuff(cbuffindex) = delZnorm
198 & * globfld3d(i,bi,ip,j,bj,jp,k)
199 # ifdef CTRL_PACK_PRECISE
200 & * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
201 # else
202 & * sqrt(weightfld(k,bi,bj))
203 # endif
204 else
205 cbuff(cbuffindex) = delZnorm
206 & * globfld3d(i,bi,ip,j,bj,jp,k)
207 # ifdef CTRL_PACK_PRECISE
208 & / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
209 # else
210 & / sqrt(weightfld(k,bi,bj))
211 # endif
212 endif
213 cph(
214 globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
215 cph)
216 #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
217 cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
218 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
219 #ifdef ALLOW_ADMTLM
220 nveccount = nveccount + 1
221 phtmpadmtlm(nveccount) = cbuff(cbuffindex)
222 #endif
223 endif
224 enddo
225 enddo
226 enddo
227 enddo
228 enddo
229 enddo
230 c --> check cbuffindex.
231 if ( cbuffindex .gt. 0) then
232 #ifndef ALLOW_ADMTLM
233 write(cunit) cbuffindex
234 write(cunit) k
235 #endif
236 write(cunit) (cbuff(ii), ii=1,cbuffindex)
237 endif
238 c
239 if ( doPackDiag ) then
240 write(cunit2,rec=irectrue) globfldtmp2
241 write(cunit3,rec=irectrue) globfldtmp3
242 endif
243 c
244 enddo
245 c
246 c -- end of irec loop --
247 enddo
248
249 if ( doPackDiag ) then
250 close ( cunit2 )
251 close ( cunit3 )
252 endif
253
254 _END_MASTER( mythid )
255
256 return
257 end
258

  ViewVC Help
Powered by ViewVC 1.1.22