/[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.16 - (show annotations) (download)
Sat Jul 15 00:27:37 2006 UTC (17 years, 10 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint59, checkpoint58y_post, checkpoint58t_post, checkpoint58m_post, checkpoint58w_post, checkpoint58o_post, checkpoint58p_post, checkpoint58q_post, checkpoint58r_post, checkpoint58n_post, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint58v_post, checkpoint58x_post, checkpoint58u_post, checkpoint58s_post
Changes since 1.15: +1 -1 lines
bug, in preconditioning by delZnorm, at iteration 0

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 delZnorm = (delR(1)/delR(k))**delZexp
178 else
179 delZnorm = 1. _d 0
180 endif
181 cbuffindex = 0
182 do jp = 1,nPy
183 do bj = jtlo,jthi
184 do j = jmin,jmax
185 do ip = 1,nPx
186 do bi = itlo,ithi
187 do i = imin,imax
188 if (globmsk(i,bi,ip,j,bj,jp,k) .ne. 0. ) then
189 cbuffindex = cbuffindex + 1
190 cph(
191 globfldtmp3(i,bi,ip,j,bj,jp) =
192 & globfld3d(i,bi,ip,j,bj,jp,k)
193 cph)
194 #ifdef ALLOW_NONDIMENSIONAL_CONTROL_IO
195 if (lxxadxx) then
196 cbuff(cbuffindex) = 1/delZnorm
197 & * globfld3d(i,bi,ip,j,bj,jp,k)
198 # ifdef CTRL_PACK_PRECISE
199 & * sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
200 # else
201 & * sqrt(weightfld(k,bi,bj))
202 # endif
203 else
204 cbuff(cbuffindex) = delZnorm
205 & * globfld3d(i,bi,ip,j,bj,jp,k)
206 # ifdef CTRL_PACK_PRECISE
207 & / sqrt(weightfld3d(i,bi,ip,j,bj,jp,k))
208 # else
209 & / sqrt(weightfld(k,bi,bj))
210 # endif
211 endif
212 cph(
213 globfldtmp2(i,bi,ip,j,bj,jp) = cbuff(cbuffindex)
214 cph)
215 #else /* ALLOW_NONDIMENSIONAL_CONTROL_IO undef */
216 cbuff(cbuffindex) = globfld3d(i,bi,ip,j,bj,jp,k)
217 #endif /* ALLOW_NONDIMENSIONAL_CONTROL_IO */
218 #ifdef ALLOW_ADMTLM
219 nveccount = nveccount + 1
220 phtmpadmtlm(nveccount) = cbuff(cbuffindex)
221 #endif
222 endif
223 enddo
224 enddo
225 enddo
226 enddo
227 enddo
228 enddo
229 c --> check cbuffindex.
230 if ( cbuffindex .gt. 0) then
231 #ifndef ALLOW_ADMTLM
232 write(cunit) cbuffindex
233 write(cunit) k
234 cph#endif
235 write(cunit) (cbuff(ii), ii=1,cbuffindex)
236 #endif
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