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

Annotation of /MITgcm/pkg/ctrl/ctrl_init.F

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


Revision 1.2 - (hide annotations) (download)
Fri Jul 13 13:40:17 2001 UTC (22 years, 10 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint40pre3, checkpoint40pre2, checkpoint40pre4, checkpoint40pre5
Changes since 1.1: +57 -13 lines
o Added prototype routines to handle optimization
o Extended control vector to add passive tracer

1 heimbach 1.2 C $Header: /u/gcmpack/models/MITgcmUV/pkg/ctrl/ctrl_init.F,v 1.1 2001/03/25 22:33:55 heimbach Exp $
2 heimbach 1.1
3     #include "CTRL_CPPOPTIONS.h"
4    
5    
6     subroutine ctrl_Init(
7     I mythid
8     & )
9    
10     c ==================================================================
11     c SUBROUTINE ctrl_Init
12     c ==================================================================
13     c
14     c o Set parts of the vector of control variables and initialize the
15     c rest to zero.
16     c
17     c The vector of control variables is initialized here. The
18     c temperature and salinity contributions are read from file.
19     c Subsequently, the latter are dimensionalized and the tile
20     c edges are updated.
21     c
22     c started: Christian Eckert eckert@mit.edu 30-Jun-1999
23     c
24     c changed: Christian Eckert eckert@mit.edu 23-Feb-2000
25     c - Restructured the code in order to create a package
26     c for the MITgcmUV.
27     c
28     c Patrick Heimbach heimbach@mit.edu 30-May-2000
29     c - diffsec was falsely declared.
30     c
31     c Patrick Heimbach heimbach@mit.edu 06-Jun-2000
32     c - Transferred some filename declarations
33     c from ctrl_pack/ctrl_unpack to here
34     c - Transferred mask-per-tile to here
35     c - computation of control vector length here
36     c
37     c Patrick Heimbach heimbach@mit.edu 16-Jun-2000
38     c - Added call to ctrl_pack
39     c - Alternatively: transfer writing of scale files to
40     c ctrl_unpack
41     c
42     c ==================================================================
43     c SUBROUTINE ctrl_Init
44     c ==================================================================
45    
46     implicit none
47    
48     c == global variables ==
49    
50     #include "EEPARAMS.h"
51     #include "SIZE.h"
52     #include "PARAMS.h"
53     #include "GRID.h"
54     #include "ctrl.h"
55    
56     c == routine arguments ==
57    
58     integer mythid
59    
60     c == local variables ==
61    
62     integer bi,bj
63     integer i,j,k
64     integer itlo,ithi
65     integer jtlo,jthi
66     integer jmin,jmax
67     integer imin,imax
68     integer ntmp
69    
70     integer il
71     integer errio
72     integer startrec
73     integer endrec
74    
75 heimbach 1.2 _RL nwetc3d
76    
77 heimbach 1.1 character*(max_len_prec) record
78     character*(max_len_mbuf) msgbuf
79    
80     c == external ==
81    
82     integer ilnblnk
83     external ilnblnk
84    
85     c == end of interface ==
86    
87     c-- Read the namelist input.
88     namelist /ctrl_nml/
89     & xx_theta_file,
90     & xx_salt_file,
91 heimbach 1.2 & xx_tr1_file,
92 heimbach 1.1 & xx_tauu_file,
93     & xx_tauv_file,
94     & xx_sflux_file,
95     & xx_hflux_file,
96     & xx_sss_file,
97     & xx_sst_file,
98     & xx_diffkr_file,
99     & xx_kapgm_file
100    
101     namelist /ctrl_packnames/
102     & yadmark, expId,
103     & ctrlname, costname, scalname, maskname, metaname
104    
105     jtlo = mybylo(mythid)
106     jthi = mybyhi(mythid)
107     itlo = mybxlo(mythid)
108     ithi = mybxhi(mythid)
109     jmin = 1-oly
110     jmax = sny+oly
111     imin = 1-olx
112     imax = snx+olx
113    
114     _BEGIN_MASTER( myThid )
115    
116     c-- Set default values.
117     xx_theta_file = ' '
118     xx_salt_file = ' '
119 heimbach 1.2 xx_tr1_file = ' '
120 heimbach 1.1 xx_tauu_file = ' '
121     xx_tauv_file = ' '
122     xx_sflux_file = ' '
123     xx_hflux_file = ' '
124     xx_sss_file = ' '
125     xx_sst_file = ' '
126     xx_diffkr_file = ' '
127     xx_kapgm_file = ' '
128     yadmark = 'ad'
129     expId = ' '
130     ctrlname = ' '
131     costname = ' '
132     scalname = ' '
133     maskname = ' '
134     metaname = ' '
135    
136     c-- Check versions.
137    
138     open(unit=scrunit1,status='scratch')
139    
140     c-- Next, read the ecco data file.
141     open(unit = modeldataunit,file = 'data.ctrl',
142     & status = 'old', iostat = errio)
143     if ( errio .lt. 0 ) then
144     stop ' stopped in ctrl_init'
145     endif
146    
147     do while ( .true. )
148     read(modeldataunit, fmt='(a)', end=1001) record
149     il = max(ilnblnk(record),1)
150     if ( record(1:1) .ne. commentcharacter )
151     & write(unit=scrunit1, fmt='(a)') record(:il)
152     enddo
153     1001 continue
154     close( modeldataunit )
155    
156     rewind( scrunit1 )
157     read(unit = scrunit1, nml = ctrl_nml)
158     read(unit = scrunit1, nml = ctrl_packnames)
159     close( scrunit1 )
160    
161     c-- Set default values.
162     do i = 1,maxcvars
163     ncvarindex(i) = -1
164     ncvarrecs(i) = 0
165     ncvarxmax(i) = 0
166     ncvarymax(i) = 0
167     ncvarnrmax(i) = 0
168     ncvargrd(i) = '?'
169     enddo
170    
171     write(msgbuf,'(a)') ' '
172     call print_message( msgbuf, standardmessageunit,
173     & SQUEEZE_RIGHT , mythid)
174     write(msgbuf,'(a)')
175     & ' ctrl_init: Initializing temperature and salinity'
176     call print_message( msgbuf, standardmessageunit,
177     & SQUEEZE_RIGHT , mythid)
178     write(msgbuf,'(a)')
179     & ' part of the control vector.'
180     call print_message( msgbuf, standardmessageunit,
181     & SQUEEZE_RIGHT , mythid)
182     write(msgbuf,'(a,a)')
183     & ' The initial surface fluxes are set',
184     & ' to zero.'
185     call print_message( msgbuf, standardmessageunit,
186     & SQUEEZE_RIGHT , mythid)
187     write(msgbuf,'(a)') ' '
188     call print_message( msgbuf, standardmessageunit,
189     & SQUEEZE_RIGHT , mythid)
190     _END_MASTER( mythid )
191    
192     _BARRIER
193    
194     c-- =====================
195     c-- Initial state fields.
196     c-- =====================
197    
198     #ifdef ALLOW_THETA0_CONTROL
199     _BEGIN_MASTER( mythid )
200     ncvarindex(1) = 101
201     ncvarrecs(1) = 1
202     ncvarxmax(1) = snx
203     ncvarymax(1) = sny
204     ncvarnrmax(1) = nr
205     ncvargrd(1) = 'c'
206     _END_MASTER( mythid )
207     #endif /* ALLOW_THETA0_CONTROL */
208    
209     #ifdef ALLOW_SALT0_CONTROL
210     _BEGIN_MASTER( mythid )
211     ncvarindex(2) = 102
212     ncvarrecs(2) = 1
213     ncvarxmax(2) = snx
214     ncvarymax(2) = sny
215     ncvarnrmax(2) = nr
216     ncvargrd(2) = 'c'
217     _END_MASTER( mythid )
218     #endif /* ALLOW_SALT0_CONTROL */
219    
220     #ifdef ALLOW_HFLUX0_CONTROL
221     _BEGIN_MASTER( mythid )
222     ncvarindex(3) = 103
223     ncvarrecs(3) = 1
224     ncvarxmax(3) = snx
225     ncvarymax(3) = sny
226     ncvarnrmax(3) = 1
227     ncvargrd(3) = 'c'
228     _END_MASTER( mythid )
229     #endif /* ALLOW_HFLUX0_CONTROL */
230    
231     #ifdef ALLOW_SFLUX0_CONTROL
232     _BEGIN_MASTER( mythid )
233     ncvarindex(4) = 104
234     ncvarrecs(4) = 1
235     ncvarxmax(4) = snx
236     ncvarymax(4) = sny
237     ncvarnrmax(4) = 1
238     ncvargrd(4) = 'c'
239     _END_MASTER( mythid )
240     #endif /* ALLOW_SFLUX0_CONTROL */
241    
242     #ifdef ALLOW_TAUU0_CONTROL
243     _BEGIN_MASTER( mythid )
244     ncvarindex(5) = 105
245     ncvarrecs(5) = 1
246     ncvarxmax(5) = snx
247     ncvarymax(5) = sny
248     ncvarnrmax(5) = 1
249     ncvargrd(5) = 'w'
250     _END_MASTER( mythid )
251     #endif /* ALLOW_TAUU0_CONTROL */
252    
253     #ifdef ALLOW_TAUV0_CONTROL
254     _BEGIN_MASTER( mythid )
255     ncvarindex(6) = 106
256     ncvarrecs(6) = 1
257     ncvarxmax(6) = snx
258     ncvarymax(6) = sny
259     ncvarnrmax(6) = 1
260     ncvargrd(6) = 's'
261     _END_MASTER( mythid )
262     #endif /* ALLOW_TAUV0_CONTROL */
263    
264 heimbach 1.2 cph(
265     cph index 7-10 reserved for atmos. state,
266     cph index 11-14 reserved for open boundaries,
267     cph index 15-16 reserved for mixing coeff.
268     cph index 17 reserved for passive tracer TR1
269     cph index 18,19 reserved for sst, sss
270     cph)
271    
272     #ifdef ALLOW_TR10_CONTROL
273     _BEGIN_MASTER( mythid )
274     ncvarindex(17) = 117
275     ncvarrecs(17) = 1
276     ncvarxmax(17) = snx
277     ncvarymax(17) = sny
278     ncvarnrmax(17) = nr
279     ncvargrd(17) = 'c'
280     _END_MASTER( mythid )
281     #endif /* ALLOW_TR10_CONTROL */
282    
283 heimbach 1.1 #ifdef ALLOW_SST0_CONTROL
284     _BEGIN_MASTER( mythid )
285 heimbach 1.2 ncvarindex(18) = 118
286     ncvarrecs(18) = 1
287     ncvarxmax(18) = snx
288     ncvarymax(18) = sny
289     ncvarnrmax(18) = 1
290     ncvargrd(18) = 'c'
291 heimbach 1.1 _END_MASTER( mythid )
292     #endif /* ALLOW_SST0_CONTROL */
293    
294     #ifdef ALLOW_SSS0_CONTROL
295     _BEGIN_MASTER( mythid )
296 heimbach 1.2 ncvarindex(19) = 119
297     ncvarrecs(19) = 1
298     ncvarxmax(19) = snx
299     ncvarymax(19) = sny
300     ncvarnrmax(19) = 1
301     ncvargrd(19) = 'c'
302 heimbach 1.1 _END_MASTER( mythid )
303     #endif /* ALLOW_SSS0_CONTROL */
304    
305     c-- Determine the number of wet points in each tile:
306     c-- maskc, masks, and maskw.
307    
308     c-- Set loop ranges.
309     jmin = 1
310     jmax = sny
311     imin = 1
312     imax = snx
313    
314     c-- Initialise the counters.
315     do bj = jtlo,jthi
316     do bi = itlo,ithi
317     do k = 1,nr
318     nwetctile(bi,bj,k) = 0
319     nwetstile(bi,bj,k) = 0
320     nwetwtile(bi,bj,k) = 0
321     enddo
322     enddo
323     enddo
324    
325     c-- Count wet points on each tile.
326     do bj = jtlo,jthi
327     do bi = itlo,ithi
328     do k = 1,nr
329     do j = jmin,jmax
330     do i = imin,imax
331     c-- Center mask.
332     if (hFacC(i,j,k,bi,bj) .ne. 0.) then
333     nwetctile(bi,bj,k) = nwetctile(bi,bj,k) + 1
334     endif
335     c-- South mask.
336     if (maskS(i,j,k,bi,bj) .eq. 1.) then
337     nwetstile(bi,bj,k) = nwetstile(bi,bj,k) + 1
338     endif
339     c-- West mask.
340     if (maskW(i,j,k,bi,bj) .eq. 1.) then
341     nwetwtile(bi,bj,k) = nwetwtile(bi,bj,k) + 1
342     endif
343     enddo
344     enddo
345     enddo
346     enddo
347     enddo
348    
349    
350     _BEGIN_MASTER( mythid )
351     c-- Determine the total number of control variables.
352     nvartype = 0
353     nvarlength = 0
354     do i = 1,maxcvars
355     if ( ncvarindex(i) .ne. -1 ) then
356     nvartype = nvartype + 1
357     do bj = jtlo,jthi
358     do bi = itlo,ithi
359     if ( ncvargrd(i) .eq. 'c' ) then
360     do k = 1,ncvarnrmax(i)
361     nvarlength = nvarlength +
362     & ncvarrecs(i)*nwetctile(bi,bj,k)
363     enddo
364     else if ( ncvargrd(i) .eq. 's' ) then
365     do k = 1,ncvarnrmax(i)
366     nvarlength = nvarlength +
367     & ncvarrecs(i)*nwetstile(bi,bj,k)
368     enddo
369     else if ( ncvargrd(i) .eq. 'w' ) then
370     do k = 1,ncvarnrmax(i)
371     nvarlength = nvarlength +
372     & ncvarrecs(i)*nwetwtile(bi,bj,k)
373     enddo
374     else
375     print*,'ctrl_init: invalid grid location'
376     print*,' control variable = ',ncvarindex(i)
377     print*,' grid location = ',ncvargrd(i)
378     stop ' ... stopped in ctrl_init'
379     endif
380     enddo
381     enddo
382     endif
383     enddo
384    
385 heimbach 1.2 cph(
386     print *, 'ph-wet 1: nvarlength = ', nvarlength
387     print *, 'ph-wet 2: surface wet C = ', nwetctile(1,1,1)
388     print *, 'ph-wet 3: surface wet W = ', nwetwtile(1,1,1)
389     print *, 'ph-wet 4: surface wet S = ', nwetstile(1,1,1)
390     nwetc3d = 0
391     do k = 1, Nr
392     nwetc3d = nwetc3d + nwetctile(1,1,k)
393     end do
394     print *, 'ph-wet 5: 3D center wet points = ', nwetc3d
395     do i = 1, 6
396     print *, 'ph-wet 6: no recs for i = ', i, ncvarrecs(i)
397     end do
398     print *, 'ph-wet 7: ',
399     & 2*nwetc3d +
400     & ncvarrecs(3)*nwetctile(1,1,1) +
401     & ncvarrecs(4)*nwetctile(1,1,1) +
402     & ncvarrecs(5)*nwetwtile(1,1,1) +
403     & ncvarrecs(6)*nwetstile(1,1,1)
404     cph)
405    
406 heimbach 1.1 c
407     c Summation of wet point counters
408     c
409     CALL GLOBAL_SUM_INT( nvarlength, myThid )
410     ntmp=0
411     do bj=1,nSy
412     do bi=1,nSx
413     ntmp=ntmp+nWetcTile(bi,bj,k)
414     enddo
415     enddo
416     CALL GLOBAL_SUM_INT( ntmp, myThid )
417     nWetcTile(1,1,k)=ntmp
418     ntmp=0
419     do bj=1,nSy
420     do bi=1,nSx
421     ntmp=ntmp+nWetsTile(bi,bj,k)
422     enddo
423     enddo
424     CALL GLOBAL_SUM_INT( ntmp, myThid )
425     nWetsTile(1,1,k)=ntmp
426     ntmp=0
427     do bj=1,nSy
428     do bi=1,nSx
429     ntmp=ntmp+nWetwTile(bi,bj,k)
430     enddo
431     enddo
432     CALL GLOBAL_SUM_INT( ntmp, myThid )
433     nWetwTile(1,1,k)=ntmp
434    
435     print*, 'ctrl_init: no. of control variables: ', nvartype
436     print*, 'ctrl_init: control vector length: ', nvarlength
437     _END_MASTER( mythid )
438    
439     _BARRIER
440    
441     return
442     end
443    

  ViewVC Help
Powered by ViewVC 1.1.22