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

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

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


Revision 1.10 - (show annotations) (download)
Thu Oct 23 04:41:40 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51o_pre, checkpoint51n_post, checkpoint51n_pre, checkpoint51o_post, checkpoint51p_post
Branch point for: checkpoint51n_branch
Changes since 1.9: +4 -1 lines
 o added the [#include "AD_CONFIG.h"] statement to all files that need
   it for adjoint/tl #defines
 o re-worked the build logic in genmake2 to support AD_CONFIG.h
 o removed tools/genmake since it no longer works

1 C
2 C $Header: $
3 C $Name: $
4
5 #include "AD_CONFIG.h"
6 #include "CTRL_CPPOPTIONS.h"
7
8
9 subroutine ctrl_init( mythid )
10
11 c ==================================================================
12 c SUBROUTINE ctrl_init
13 c ==================================================================
14 c
15 c o Set parts of the vector of control variables and initialize the
16 c rest to zero.
17 c
18 c The vector of control variables is initialized here. The
19 c temperature and salinity contributions are read from file.
20 c Subsequently, the latter are dimensionalized and the tile
21 c edges are updated.
22 c
23 c started: Christian Eckert eckert@mit.edu 30-Jun-1999
24 c
25 c changed: Christian Eckert eckert@mit.edu 23-Feb-2000
26 c - Restructured the code in order to create a package
27 c for the MITgcmUV.
28 c
29 c Patrick Heimbach heimbach@mit.edu 30-May-2000
30 c - diffsec was falsely declared.
31 c
32 c Patrick Heimbach heimbach@mit.edu 06-Jun-2000
33 c - Transferred some filename declarations
34 c from ctrl_pack/ctrl_unpack to here
35 c - Transferred mask-per-tile to here
36 c - computation of control vector length here
37 c
38 c Patrick Heimbach heimbach@mit.edu 16-Jun-2000
39 c - Added call to ctrl_pack
40 c - Alternatively: transfer writing of scale files to
41 c ctrl_unpack
42 c
43 c Dimitris Menemenlis menemenlis@mit.edu 7-Mar-2003
44 c - To be consistent with usage in ctrl_getrec.F,
45 c startrec and endrec need to be referenced to
46 c model time = 0, not to startTime.
47 c Also "- modelstep" -> "+ modelstep/2":
48 c old: startrec = int((modelstart - diffsecs)/
49 c old: & xx_???period) + 1
50 c old: endrec = int((modelend - diffsecs - modelstep)/
51 c old: & xx_???period) + 2
52 c new: startrec = int((modelstart + startTime - diffsecs)/
53 c new: & xx_???period) + 1
54 c new: endrec = int((modelend + startTime - diffsecs + modelstep/2)/
55 c new: & xx_???period) + 2
56 c
57 c ==================================================================
58 c SUBROUTINE ctrl_init
59 c ==================================================================
60
61 implicit none
62
63 c == global variables ==
64
65 #include "EEPARAMS.h"
66 #include "SIZE.h"
67 #include "PARAMS.h"
68 #include "GRID.h"
69 #include "ctrl.h"
70
71 #ifdef ALLOW_CALENDAR
72 # include "cal.h"
73 #endif
74 #ifdef ALLOW_OBCS_CONTROL
75 # include "OBCS.h"
76 #endif
77 #ifdef ALLOW_ECCO_OPTIMIZATION
78 # include "optim.h"
79 #endif
80
81 c == routine arguments ==
82
83 integer mythid
84
85 c == local variables ==
86
87 integer bi,bj
88 integer i,j,k
89 integer itlo,ithi
90 integer jtlo,jthi
91 integer jmin,jmax
92 integer imin,imax
93 integer ntmp
94 integer ivarindex
95
96 integer iobcs
97 integer il
98 integer errio
99 integer startrec
100 integer endrec
101 integer difftime(4)
102 _RL diffsecs
103 _RL dummy
104
105 character*(80) ymaskobcs
106 character*(max_len_prec) record
107 character*(max_len_mbuf) msgbuf
108
109 integer nwetc3d
110
111 c == external ==
112
113 integer ilnblnk
114 external ilnblnk
115
116 c == end of interface ==
117
118 c-- Set loop ranges.
119 jtlo = mybylo(mythid)
120 jthi = mybyhi(mythid)
121 itlo = mybxlo(mythid)
122 ithi = mybxhi(mythid)
123 jmin = 1
124 jmax = sny
125 imin = 1
126 imax = snx
127
128
129 _BEGIN_MASTER( myThid )
130
131 #ifdef ALLOW_CALENDAR
132
133 c-- Get the complete dates of the control variables.
134 #if (defined (ALLOW_HFLUX_CONTROL))
135 c-- The heat flux contribution.
136 call cal_FullDate( xx_hfluxstartdate1, xx_hfluxstartdate2,
137 & xx_hfluxstartdate , mythid )
138 #elif (defined (ALLOW_ATEMP_CONTROL))
139 c-- Atmos. temperature contribution.
140 call cal_FullDate( xx_atempstartdate1, xx_atempstartdate2,
141 & xx_atempstartdate , mythid )
142 #endif
143
144 #if (defined (ALLOW_SFLUX_CONTROL))
145 c-- The salt flux contribution.
146 call cal_FullDate( xx_sfluxstartdate1, xx_sfluxstartdate2,
147 & xx_sfluxstartdate , mythid )
148 #elif (defined (ALLOW_AQH_CONTROL))
149 c-- Atmospheric humidity contribution.
150 call cal_FullDate( xx_aqhstartdate1, xx_aqhstartdate2,
151 & xx_aqhstartdate , mythid )
152 #endif
153
154 #if (defined (ALLOW_USTRESS_CONTROL))
155 c-- The zonal wind stress contribution.
156 call cal_FullDate( xx_tauustartdate1, xx_tauustartdate2,
157 & xx_tauustartdate, mythid )
158 #elif (defined (ALLOW_UWIND_CONTROL))
159 c-- Zonal wind speed contribution.
160 call cal_FullDate( xx_uwindstartdate1, xx_uwindstartdate2,
161 & xx_uwindstartdate , mythid )
162 #endif
163
164 #if (defined (ALLOW_VSTRESS_CONTROL))
165 c-- The merid. wind stress contribution.
166 call cal_FullDate( xx_tauvstartdate1, xx_tauvstartdate2,
167 & xx_tauvstartdate, mythid )
168 #elif (defined (ALLOW_VWIND_CONTROL))
169 c-- Merid. wind speed contribution.
170 call cal_FullDate( xx_vwindstartdate1, xx_vwindstartdate2,
171 & xx_vwindstartdate , mythid )
172 #endif
173
174 #ifdef ALLOW_OBCS_CONTROL
175 call cal_FullDate( xx_obcsnstartdate1, xx_obcsnstartdate2,
176 & xx_obcsnstartdate, mythid )
177 call cal_FullDate( xx_obcssstartdate1, xx_obcssstartdate2,
178 & xx_obcssstartdate, mythid )
179 call cal_FullDate( xx_obcswstartdate1, xx_obcswstartdate2,
180 & xx_obcswstartdate, mythid )
181 call cal_FullDate( xx_obcsestartdate1, xx_obcsestartdate2,
182 & xx_obcsestartdate, mythid )
183 #endif
184
185 #endif /* ALLOW_CALENDAR */
186
187 c-- Set default values.
188 do ivarindex = 1,maxcvars
189 ncvarindex(ivarindex) = -1
190 ncvarrecs(ivarindex) = 0
191 ncvarxmax(ivarindex) = 0
192 ncvarymax(ivarindex) = 0
193 ncvarnrmax(ivarindex) = 0
194 ncvargrd(ivarindex) = '?'
195 enddo
196
197 write(msgbuf,'(a)') ' '
198 call print_message( msgbuf, standardmessageunit,
199 & SQUEEZE_RIGHT , mythid)
200 write(msgbuf,'(a)')
201 & ' ctrl_init: Initializing temperature and salinity'
202 call print_message( msgbuf, standardmessageunit,
203 & SQUEEZE_RIGHT , mythid)
204 write(msgbuf,'(a)')
205 & ' part of the control vector.'
206 call print_message( msgbuf, standardmessageunit,
207 & SQUEEZE_RIGHT , mythid)
208 write(msgbuf,'(a,a)')
209 & ' The initial surface fluxes are set',
210 & ' to zero.'
211 call print_message( msgbuf, standardmessageunit,
212 & SQUEEZE_RIGHT , mythid)
213 write(msgbuf,'(a)') ' '
214 call print_message( msgbuf, standardmessageunit,
215 & SQUEEZE_RIGHT , mythid)
216 _END_MASTER( mythid )
217
218 _BARRIER
219
220 c-- =====================
221 c-- Initial state fields.
222 c-- =====================
223
224 cph(
225 cph index 7-10 reserved for atmos. state,
226 cph index 11-14 reserved for open boundaries,
227 cph index 15-16 reserved for mixing coeff.
228 cph index 17 reserved for passive tracer TR1
229 cph index 18,19 reserved for sst, sss
230 cph index 20 for hFacC
231 cph index 21-22 for efluxy, efluxp
232 cph index 23-24 for bottom drag
233 cph)
234
235 c-------------------------------------------------------------------------------
236 c--
237 #ifdef ALLOW_THETA0_CONTROL
238 c-- Initial state temperature contribution.
239
240 _BEGIN_MASTER( mythid )
241 ivarindex = 1
242 ncvarindex(ivarindex) = 101
243 ncvarrecs(ivarindex) = 1
244 ncvarxmax(ivarindex) = snx
245 ncvarymax(ivarindex) = sny
246 ncvarnrmax(ivarindex) = nr
247 ncvargrd(ivarindex) = 'c'
248 _END_MASTER( mythid )
249
250 #endif /* ALLOW_THETA0_CONTROL */
251
252 c-------------------------------------------------------------------------------
253 c--
254 #ifdef ALLOW_SALT0_CONTROL
255 c-- Initial state salinity contribution.
256
257 _BEGIN_MASTER( mythid )
258 ivarindex = 2
259 ncvarindex(ivarindex) = 102
260 ncvarrecs(ivarindex) = 1
261 ncvarxmax(ivarindex) = snx
262 ncvarymax(ivarindex) = sny
263 ncvarnrmax(ivarindex) = nr
264 ncvargrd(ivarindex) = 'c'
265 _END_MASTER( mythid )
266
267 #endif /* ALLOW_SALT0_CONTROL */
268
269 c-- ===========================
270 c-- Surface flux contributions.
271 c-- ===========================
272
273 c-------------------------------------------------------------------------------
274 c--
275 #if (defined (ALLOW_HFLUX_CONTROL))
276 c-- Heat flux.
277
278 _BEGIN_MASTER( mythid )
279 #ifdef ALLOW_CALENDAR
280 call cal_TimePassed( xx_hfluxstartdate, modelstartdate,
281 & difftime, mythid )
282 call cal_ToSeconds ( difftime, diffsecs, mythid )
283 startrec = int((modelstart + startTime - diffsecs)/
284 & xx_hfluxperiod) + 1
285 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
286 & xx_hfluxperiod) + 2
287 #else
288 startrec = 1
289 endrec = 1
290 #endif
291 ivarindex = 3
292 ncvarindex(ivarindex) = 103
293 ncvarrecs(ivarindex) = endrec - startrec + 1
294 ncvarrecstart(ivarindex) = startrec
295 ncvarrecsend(ivarindex) = endrec
296 ncvarxmax(ivarindex) = snx
297 ncvarymax(ivarindex) = sny
298 ncvarnrmax(ivarindex) = 1
299 ncvargrd(ivarindex) = 'c'
300 _END_MASTER( mythid )
301
302 #elif (defined (ALLOW_ATEMP_CONTROL))
303 c-- Atmos. temperature
304
305 _BEGIN_MASTER( mythid )
306 #ifdef ALLOW_CALENDAR
307 call cal_TimePassed( xx_atempstartdate, modelstartdate,
308 & difftime, mythid )
309 call cal_ToSeconds ( difftime, diffsecs, mythid )
310 startrec = int((modelstart + startTime - diffsecs)/
311 & xx_atempperiod) + 1
312 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
313 & xx_atempperiod) + 2
314 #else
315 startrec = 1
316 endrec = 1
317 #endif
318 ivarindex = 7
319 ncvarindex(ivarindex) = 107
320 ncvarrecs(ivarindex) = endrec - startrec + 1
321 ncvarrecstart(ivarindex) = startrec
322 ncvarrecsend(ivarindex) = endrec
323 ncvarxmax(ivarindex) = snx
324 ncvarymax(ivarindex) = sny
325 ncvarnrmax(ivarindex) = 1
326 ncvargrd(ivarindex) = 'c'
327 _END_MASTER( mythid )
328
329 #elif (defined (ALLOW_HFLUX0_CONTROL))
330 c-- initial forcing only
331 _BEGIN_MASTER( mythid )
332 ncvarindex(3) = 103
333 ncvarrecs(3) = 1
334 ncvarxmax(3) = snx
335 ncvarymax(3) = sny
336 ncvarnrmax(3) = 1
337 ncvargrd(3) = 'c'
338 _END_MASTER( mythid )
339
340 #endif /* ALLOW_HFLUX_CONTROL */
341
342 c-------------------------------------------------------------------------------
343 c--
344 #if (defined (ALLOW_SFLUX_CONTROL))
345 c-- Salt flux.
346
347 _BEGIN_MASTER( mythid )
348 #ifdef ALLOW_CALENDAR
349 call cal_TimePassed( xx_sfluxstartdate, modelstartdate,
350 & difftime, mythid )
351 call cal_ToSeconds ( difftime, diffsecs, mythid )
352 startrec = int((modelstart + startTime - diffsecs)/
353 & xx_sfluxperiod) + 1
354 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
355 & xx_sfluxperiod) + 2
356 #else
357 startrec = 1
358 endrec = 1
359 #endif
360 ivarindex = 4
361 ncvarindex(ivarindex) = 104
362 ncvarrecs(ivarindex) = endrec - startrec + 1
363 ncvarrecstart(ivarindex) = startrec
364 ncvarrecsend(ivarindex) = endrec
365 ncvarxmax(ivarindex) = snx
366 ncvarymax(ivarindex) = sny
367 ncvarnrmax(ivarindex) = 1
368 ncvargrd(ivarindex) = 'c'
369 _END_MASTER( mythid )
370
371 #elif (defined (ALLOW_AQH_CONTROL))
372 c-- Atmos. humidity
373
374 _BEGIN_MASTER( mythid )
375 #ifdef ALLOW_CALENDAR
376 call cal_TimePassed( xx_aqhstartdate, modelstartdate,
377 & difftime, mythid )
378 call cal_ToSeconds ( difftime, diffsecs, mythid )
379 startrec = int((modelstart + startTime - diffsecs)/
380 & xx_aqhperiod) + 1
381 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
382 & xx_aqhperiod) + 2
383 #else
384 startrec = 1
385 endrec = 1
386 #endif
387 ivarindex = 8
388 ncvarindex(ivarindex) = 108
389 ncvarrecs(ivarindex) = endrec - startrec + 1
390 ncvarrecstart(ivarindex) = startrec
391 ncvarrecsend(ivarindex) = endrec
392 ncvarxmax(ivarindex) = snx
393 ncvarymax(ivarindex) = sny
394 ncvarnrmax(ivarindex) = 1
395 ncvargrd(ivarindex) = 'c'
396 _END_MASTER( mythid )
397
398 #elif (defined (ALLOW_SFLUX0_CONTROL))
399 c-- initial forcing only
400 _BEGIN_MASTER( mythid )
401 ncvarindex(4) = 104
402 ncvarrecs(4) = 1
403 ncvarxmax(4) = snx
404 ncvarymax(4) = sny
405 ncvarnrmax(4) = 1
406 ncvargrd(4) = 'c'
407 _END_MASTER( mythid )
408
409 #endif /* ALLOW_SFLUX_CONTROL */
410
411 c-------------------------------------------------------------------------------
412 c--
413 #if (defined (ALLOW_USTRESS_CONTROL))
414 c-- Zonal wind stress.
415
416 _BEGIN_MASTER( mythid )
417 #ifdef ALLOW_CALENDAR
418 call cal_TimePassed( xx_tauustartdate, modelstartdate,
419 & difftime, mythid )
420 call cal_ToSeconds ( difftime, diffsecs, mythid )
421 startrec = int((modelstart + startTime - diffsecs)/
422 & xx_tauuperiod) + 1
423 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
424 & xx_tauuperiod) + 2
425 #else
426 startrec = 1
427 endrec = 1
428 #endif
429 ivarindex = 5
430 ncvarindex(ivarindex) = 105
431 ncvarrecs(ivarindex) = endrec - startrec + 1
432 ncvarrecstart(ivarindex) = startrec
433 ncvarrecsend(ivarindex) = endrec
434 ncvarxmax(ivarindex) = snx
435 ncvarymax(ivarindex) = sny
436 ncvarnrmax(ivarindex) = 1
437 ncvargrd(ivarindex) = 'w'
438 _END_MASTER( mythid )
439
440 #elif (defined (ALLOW_UWIND_CONTROL))
441 c-- Zonal wind speed.
442
443 _BEGIN_MASTER( mythid )
444 #ifdef ALLOW_CALENDAR
445 call cal_TimePassed( xx_uwindstartdate, modelstartdate,
446 & difftime, mythid )
447 call cal_ToSeconds ( difftime, diffsecs, mythid )
448 startrec = int((modelstart + startTime - diffsecs)/
449 & xx_uwindperiod) + 1
450 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
451 & xx_uwindperiod) + 2
452 #else
453 startrec = 1
454 endrec = 1
455 #endif
456 ivarindex = 9
457 ncvarindex(ivarindex) = 109
458 ncvarrecs(ivarindex) = endrec - startrec + 1
459 ncvarrecstart(ivarindex) = startrec
460 ncvarrecsend(ivarindex) = endrec
461 ncvarxmax(ivarindex) = snx
462 ncvarymax(ivarindex) = sny
463 ncvarnrmax(ivarindex) = 1
464 ncvargrd(ivarindex) = 'w'
465 _END_MASTER( mythid )
466
467 #elif (defined (ALLOW_TAUU0_CONTROL))
468 c-- initial forcing only
469 _BEGIN_MASTER( mythid )
470 ncvarindex(5) = 105
471 ncvarrecs(5) = 1
472 ncvarxmax(5) = snx
473 ncvarymax(5) = sny
474 ncvarnrmax(5) = 1
475 ncvargrd(5) = 'w'
476 _END_MASTER( mythid )
477
478 #endif /* ALLOW_USTRESS_CONTROL */
479
480 c-------------------------------------------------------------------------------
481 c--
482 #if (defined (ALLOW_VSTRESS_CONTROL))
483 c-- Meridional wind stress.
484
485 _BEGIN_MASTER( mythid )
486 #ifdef ALLOW_CALENDAR
487 call cal_TimePassed( xx_tauvstartdate, modelstartdate,
488 & difftime, mythid )
489 call cal_ToSeconds ( difftime, diffsecs, mythid )
490 startrec = int((modelstart + startTime - diffsecs)/
491 & xx_tauvperiod) + 1
492 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
493 & xx_tauvperiod) + 2
494 #else
495 startrec = 1
496 endrec = 1
497 #endif
498 ivarindex = 6
499 ncvarindex(ivarindex) = 106
500 ncvarrecs(ivarindex) = endrec - startrec + 1
501 ncvarrecstart(ivarindex) = startrec
502 ncvarrecsend(ivarindex) = endrec
503 ncvarxmax(ivarindex) = snx
504 ncvarymax(ivarindex) = sny
505 ncvarnrmax(ivarindex) = 1
506 ncvargrd(ivarindex) = 's'
507 _END_MASTER( mythid )
508
509 #elif (defined (ALLOW_VWIND_CONTROL))
510 c-- Meridional wind speed.
511
512 _BEGIN_MASTER( mythid )
513 #ifdef ALLOW_CALENDAR
514 call cal_TimePassed( xx_vwindstartdate, modelstartdate,
515 & difftime, mythid )
516 call cal_ToSeconds ( difftime, diffsecs, mythid )
517 startrec = int((modelstart + startTime - diffsecs)/
518 & xx_vwindperiod) + 1
519 endrec = int((modelend + startTime - diffsecs + modelstep/2)/
520 & xx_vwindperiod) + 2
521 #else
522 startrec = 1
523 endrec = 1
524 #endif
525 ivarindex = 10
526 ncvarindex(ivarindex) = 110
527 ncvarrecs(ivarindex) = endrec - startrec + 1
528 ncvarrecstart(ivarindex) = startrec
529 ncvarrecsend(ivarindex) = endrec
530 ncvarxmax(ivarindex) = snx
531 ncvarymax(ivarindex) = sny
532 ncvarnrmax(ivarindex) = 1
533 ncvargrd(ivarindex) = 's'
534 _END_MASTER( mythid )
535
536 #elif (defined (ALLOW_TAUV0_CONTROL))
537 c-- initial forcing only
538 _BEGIN_MASTER( mythid )
539 ncvarindex(6) = 106
540 ncvarrecs(6) = 1
541 ncvarxmax(6) = snx
542 ncvarymax(6) = sny
543 ncvarnrmax(6) = 1
544 ncvargrd(6) = 's'
545 _END_MASTER( mythid )
546
547 #endif /* ALLOW_VSTRESS_CONTROL */
548
549 c-------------------------------------------------------------------------------
550 c--
551 #ifdef ALLOW_OBCSN_CONTROL
552 c-- Northern obc.
553
554 _BEGIN_MASTER( mythid )
555 #ifdef ALLOW_CALENDAR
556 call cal_TimePassed( xx_obcsnstartdate, modelstartdate,
557 & difftime, mythid )
558 call cal_ToSeconds ( difftime, diffsecs, mythid )
559 cgg O.B. future values are needed at the last timestep, so lets
560 cgg take this into account.
561 startrec = int((modelstart - diffsecs)/xx_obcsnperiod) + 1
562 endrec = int((modelend - diffsecs)/xx_obcsnperiod) + 2
563 #else
564 startrec = 1
565 endrec = 1
566 #endif
567 ivarindex = 11
568 ncvarindex(ivarindex) = 111
569 cgg( Implement heimbach fix for nobcs.
570 ncvarrecs(ivarindex) = (endrec - startrec + 1)*nobcs
571 ncvarrecstart(ivarindex) = (startrec - 1)*nobcs + 1
572 ncvarrecsend(ivarindex) = endrec*nobcs
573 cgg)
574 ncvarxmax(ivarindex) = snx
575 ncvarymax(ivarindex) = 1
576 ncvarnrmax(ivarindex) = nr
577 ncvargrd(ivarindex) = 'm'
578 _END_MASTER( mythid )
579
580 #endif /* ALLOW_OBCSN_CONTROL */
581
582 #ifdef ALLOW_OBCSS_CONTROL
583 c-- Southern obc.
584
585 c-------------------------------------------------------------------------------
586 c--
587 _BEGIN_MASTER( mythid )
588 #ifdef ALLOW_CALENDAR
589 call cal_TimePassed( xx_obcssstartdate, modelstartdate,
590 & difftime, mythid )
591 call cal_ToSeconds ( difftime, diffsecs, mythid )
592 cgg O.B. future values are needed at the last timestep, so lets
593 cgg take this into account.
594 startrec = int((modelstart - diffsecs)/xx_obcssperiod) + 1
595 endrec = int((modelend - diffsecs)/xx_obcssperiod) + 2
596 #else
597 startrec = 1
598 endrec = 1
599 #endif
600 ivarindex = 12
601 ncvarindex(ivarindex) = 112
602 cgg( Implement heimbach fix for nobcs.
603 ncvarrecs(ivarindex) = (endrec - startrec + 1)*nobcs
604 ncvarrecstart(ivarindex) = (startrec - 1)*nobcs + 1
605 ncvarrecsend(ivarindex) = endrec*nobcs
606 cph)
607 ncvarxmax(ivarindex) = snx
608 ncvarymax(ivarindex) = 1
609 ncvarnrmax(ivarindex) = nr
610 ncvargrd(ivarindex) = 'm'
611 _END_MASTER( mythid )
612
613 #endif /* ALLOW_OBCSS_CONTROL */
614
615 c-------------------------------------------------------------------------------
616 c--
617 #ifdef ALLOW_OBCSW_CONTROL
618 c-- Western obc.
619
620 _BEGIN_MASTER( mythid )
621 #ifdef ALLOW_CALENDAR
622 call cal_TimePassed( xx_obcswstartdate, modelstartdate,
623 & difftime, mythid )
624 call cal_ToSeconds ( difftime, diffsecs, mythid )
625 startrec = int((modelstart - diffsecs)/xx_obcswperiod) + 1
626 endrec = int((modelend - diffsecs)/xx_obcswperiod) + 2
627 #else
628 startrec = 1
629 endrec = 1
630 #endif
631 ivarindex = 13
632 ncvarindex(ivarindex) = 113
633 cgg( Implement heimbach fix for nobcs.
634 ncvarrecs(ivarindex) = (endrec - startrec + 1)*nobcs
635 ncvarrecstart(ivarindex) = (startrec - 1)*nobcs + 1
636 ncvarrecsend(ivarindex) = endrec*nobcs
637 cgg)
638 ncvarxmax(ivarindex) = 1
639 ncvarymax(ivarindex) = sny
640 ncvarnrmax(ivarindex) = nr
641 ncvargrd(ivarindex) = 'm'
642 _END_MASTER( mythid )
643
644 #endif /* ALLOW_OBCSW_CONTROL */
645
646 c-------------------------------------------------------------------------------
647 c--
648 #ifdef ALLOW_OBCSE_CONTROL
649 c-- Eastern obc.
650
651 _BEGIN_MASTER( mythid )
652 #ifdef ALLOW_CALENDAR
653 call cal_TimePassed( xx_obcsestartdate, modelstartdate,
654 & difftime, mythid )
655 call cal_ToSeconds ( difftime, diffsecs, mythid )
656 startrec = int((modelstart - diffsecs)/xx_obcseperiod) + 1
657 endrec = int((modelend - diffsecs)/xx_obcseperiod) + 2
658 #else
659 startrec = 1
660 endrec = 1
661 #endif
662 ivarindex = 14
663 ncvarindex(ivarindex) = 114
664 cgg( Implement heimbach fix for nobcs.
665 ncvarrecs(ivarindex) = (endrec - startrec + 1)*nobcs
666 ncvarrecstart(ivarindex) = (startrec - 1)*nobcs + 1
667 ncvarrecsend(ivarindex) = endrec*nobcs
668 cgg)
669 ncvarxmax(ivarindex) = 1
670 ncvarymax(ivarindex) = sny
671 ncvarnrmax(ivarindex) = nr
672 ncvargrd(ivarindex) = 'm'
673 _END_MASTER( mythid )
674
675 #endif /* ALLOW_OBCSE_CONTROL */
676
677 c-------------------------------------------------------------------------------
678 c--
679 #ifdef ALLOW_DIFFKR_CONTROL
680 _BEGIN_MASTER( mythid )
681 ivarindex = 15
682 ncvarindex(ivarindex) = 115
683 ncvarrecs (ivarindex) = 1
684 ncvarxmax (ivarindex) = snx
685 ncvarymax (ivarindex) = sny
686 ncvarnrmax(ivarindex) = nr
687 ncvargrd (ivarindex) = 'c'
688 _END_MASTER( mythid )
689 #endif /* ALLOW_DIFFKR_CONTROL */
690
691 c-------------------------------------------------------------------------------
692 c--
693 #ifdef ALLOW_KAPGM_CONTROL
694 _BEGIN_MASTER( mythid )
695 ivarindex = 16
696 ncvarindex(ivarindex) = 116
697 ncvarrecs (ivarindex) = 1
698 ncvarxmax (ivarindex) = snx
699 ncvarymax (ivarindex) = sny
700 ncvarnrmax(ivarindex) = nr
701 ncvargrd (ivarindex) = 'c'
702 _END_MASTER( mythid )
703 #endif /* ALLOW_KAPGM_CONTROL */
704
705 c-------------------------------------------------------------------------------
706 c--
707 #ifdef ALLOW_TR10_CONTROL
708 _BEGIN_MASTER( mythid )
709 ivarindex = 17
710 ncvarindex(ivarindex) = 117
711 ncvarrecs (ivarindex) = 1
712 ncvarxmax (ivarindex) = snx
713 ncvarymax (ivarindex) = sny
714 ncvarnrmax(ivarindex) = nr
715 ncvargrd (ivarindex) = 'c'
716 _END_MASTER( mythid )
717 #endif /* ALLOW_TR10_CONTROL */
718
719 c-------------------------------------------------------------------------------
720 c--
721 #ifdef ALLOW_SST0_CONTROL
722 _BEGIN_MASTER( mythid )
723 ivarindex = 18
724 ncvarindex(ivarindex) = 118
725 ncvarrecs (ivarindex) = 1
726 ncvarxmax (ivarindex) = snx
727 ncvarymax (ivarindex) = sny
728 ncvarnrmax(ivarindex) = 1
729 ncvargrd (ivarindex) = 'c'
730 _END_MASTER( mythid )
731 #endif /* ALLOW_SST0_CONTROL */
732
733 c-------------------------------------------------------------------------------
734 c--
735 #ifdef ALLOW_SSS0_CONTROL
736 _BEGIN_MASTER( mythid )
737 ivarindex = 19
738 ncvarindex(ivarindex) = 119
739 ncvarrecs (ivarindex) = 1
740 ncvarxmax (ivarindex) = snx
741 ncvarymax (ivarindex) = sny
742 ncvarnrmax(ivarindex) = 1
743 ncvargrd (ivarindex) = 'c'
744 _END_MASTER( mythid )
745 #endif /* ALLOW_SSS0_CONTROL */
746
747 c-------------------------------------------------------------------------------
748 c--
749 #ifdef ALLOW_HFACC_CONTROL
750 _BEGIN_MASTER( mythid )
751 ivarindex = 20
752 ncvarindex(ivarindex) = 120
753 ncvarrecs (ivarindex) = 1
754 ncvarxmax (ivarindex) = snx
755 ncvarymax (ivarindex) = sny
756 ncvargrd (ivarindex) = 'c'
757 #ifdef ALLOW_HFACC3D_CONTROL
758 ncvarnrmax(ivarindex) = nr
759 #else
760 ncvarnrmax(ivarindex) = 1
761 #endif /*ALLOW_HFACC3D_CONTROL*/
762 _END_MASTER( mythid )
763 #endif /* ALLOW_HFACC_CONTROL */
764
765 c-------------------------------------------------------------------------------
766 c--
767 #ifdef ALLOW_EFLUXY0_CONTROL
768 _BEGIN_MASTER( mythid )
769 ivarindex = 21
770 ncvarindex(ivarindex) = 121
771 ncvarrecs(ivarindex) = 1
772 ncvarxmax(ivarindex) = snx
773 ncvarymax(ivarindex) = sny
774 ncvarnrmax(ivarindex) = nr
775 ncvargrd(ivarindex) = 's'
776 _END_MASTER( mythid )
777 #endif /* ALLOW_EFLUXY0_CONTROL */
778
779 c-------------------------------------------------------------------------------
780 c--
781 #ifdef ALLOW_EFLUXP0_CONTROL
782 _BEGIN_MASTER( mythid )
783 ivarindex = 22
784 ncvarindex(ivarindex) = 122
785 ncvarrecs(ivarindex) = 1
786 ncvarxmax(ivarindex) = snx
787 ncvarymax(ivarindex) = sny
788 ncvarnrmax(ivarindex) = nr
789 ncvargrd(ivarindex) = 'v'
790 _END_MASTER( mythid )
791 #endif /* ALLOW_EFLUXP0_CONTROL */
792
793 c-------------------------------------------------------------------------------
794 c--
795 #ifdef ALLOW_BOTTOMDRAG_CONTROL
796 _BEGIN_MASTER( mythid )
797 ivarindex = 23
798 ncvarindex(ivarindex) = 123
799 ncvarrecs (ivarindex) = 1
800 ncvarxmax (ivarindex) = snx
801 ncvarymax (ivarindex) = sny
802 ncvarnrmax(ivarindex) = 1
803 ncvargrd (ivarindex) = 'c'
804 _END_MASTER( mythid )
805 #endif /* ALLOW_BOTTOMDRAG_CONTROL */
806
807 c-------------------------------------------------------------------------------
808 c-------------------------------------------------------------------------------
809 c-------------------------------------------------------------------------------
810
811 c-- Determine the number of wet points in each tile:
812 c-- maskc, masks, and maskw.
813
814 c-- Initialise the counters.
815 do bj = jtlo,jthi
816 do bi = itlo,ithi
817 do k = 1,nr
818 nwetctile(bi,bj,k) = 0
819 nwetstile(bi,bj,k) = 0
820 nwetwtile(bi,bj,k) = 0
821 nwetvtile(bi,bj,k) = 0
822 enddo
823 enddo
824 enddo
825
826 #ifdef ALLOW_OBCS_CONTROL
827 c-- Initialise obcs counters.
828 do bj = jtlo,jthi
829 do bi = itlo,ithi
830 do k = 1,nr
831 do iobcs = 1,nobcs
832 #ifdef ALLOW_OBCSN_CONTROL
833 nwetobcsn(bi,bj,k,iobcs) = 0
834 #endif
835 #ifdef ALLOW_OBCSS_CONTROL
836 nwetobcss(bi,bj,k,iobcs) = 0
837 #endif
838 #ifdef ALLOW_OBCSW_CONTROL
839 nwetobcsw(bi,bj,k,iobcs) = 0
840 #endif
841 #ifdef ALLOW_OBCSE_CONTROL
842 nwetobcse(bi,bj,k,iobcs) = 0
843 #endif
844 enddo
845 enddo
846 enddo
847 enddo
848 #endif
849
850 c-- Count wet points on each tile.
851 do bj = jtlo,jthi
852 do bi = itlo,ithi
853 do k = 1,nr
854 do j = jmin,jmax
855 do i = imin,imax
856 c-- Center mask.
857 if (hFacC(i,j,k,bi,bj) .ne. 0.) then
858 nwetctile(bi,bj,k) = nwetctile(bi,bj,k) + 1
859 endif
860 c-- South mask.
861 if (maskS(i,j,k,bi,bj) .eq. 1.) then
862 nwetstile(bi,bj,k) = nwetstile(bi,bj,k) + 1
863 endif
864 c-- West mask.
865 if (maskW(i,j,k,bi,bj) .eq. 1.) then
866 nwetwtile(bi,bj,k) = nwetwtile(bi,bj,k) + 1
867 endif
868 #if (defined (ALLOW_EFLUXP0_CONTROL))
869 c-- Vertical mask.
870 if (hFacV(i,j,k,bi,bj) .ne. 0.) then
871 nwetvtile(bi,bj,k) = nwetvtile(bi,bj,k) + 1
872 endif
873 #endif
874 enddo
875 enddo
876 enddo
877 enddo
878 enddo
879
880 #ifdef ALLOW_OBCSN_CONTROL
881 c-- Count wet points at Northern boundary.
882 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
883 ymaskobcs = 'maskobcsn'
884 call ctrl_mask_set_xz( 0, OB_Jn, nwetobcsn, ymaskobcs, mythid )
885 #endif
886
887 #ifdef ALLOW_OBCSS_CONTROL
888 c-- Count wet points at Southern boundary.
889 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
890 ymaskobcs = 'maskobcss'
891 call ctrl_mask_set_xz( 1, OB_Js, nwetobcss, ymaskobcs, mythid )
892 #endif
893
894 #ifdef ALLOW_OBCSW_CONTROL
895 c-- Count wet points at Western boundary.
896 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
897 ymaskobcs = 'maskobcsw'
898 call ctrl_mask_set_yz( 1, OB_Iw, nwetobcsw, ymaskobcs, mythid )
899 #endif
900
901 #ifdef ALLOW_OBCSE_CONTROL
902 c-- Count wet points at Eastern boundary.
903 c-- mask conventions are adopted from obcs_apply_ts, obcs_apply_uv
904 ymaskobcs = 'maskobcse'
905 call ctrl_mask_set_yz( 0, OB_Ie, nwetobcse, ymaskobcs, mythid )
906 #endif
907
908 _BEGIN_MASTER( mythid )
909 c-- Determine the total number of control variables.
910 nvartype = 0
911 nvarlength = 0
912 do i = 1,maxcvars
913 c
914 if ( ncvarindex(i) .ne. -1 ) then
915 nvartype = nvartype + 1
916 do bj = jtlo,jthi
917 do bi = itlo,ithi
918 do k = 1,ncvarnrmax(i)
919 if ( ncvargrd(i) .eq. 'c' ) then
920 nvarlength = nvarlength +
921 & ncvarrecs(i)*nwetctile(bi,bj,k)
922 else if ( ncvargrd(i) .eq. 's' ) then
923 nvarlength = nvarlength +
924 & ncvarrecs(i)*nwetstile(bi,bj,k)
925 else if ( ncvargrd(i) .eq. 'w' ) then
926 nvarlength = nvarlength +
927 & ncvarrecs(i)*nwetwtile(bi,bj,k)
928 else if ( ncvargrd(i) .eq. 'v' ) then
929 nvarlength = nvarlength +
930 & ncvarrecs(i)*nwetvtile(bi,bj,k)
931 else if ( ncvargrd(i) .eq. 'm' ) then
932 #ifdef ALLOW_OBCS_CONTROL
933 do iobcs = 1, nobcs
934 cgg This overcounts the number of o.b. control points by a factor of "nobcs".
935 cgg As an ad-hoc solution I've divided by nobcs everywhere.
936 if ( i .eq. 11 ) then
937 #ifdef ALLOW_OBCSN_CONTROL
938 nvarlength = nvarlength +
939 & (ncvarrecs(i)/nobcs)
940 & *nwetobcsn(bi,bj,k,iobcs)
941 #endif
942 else if ( i .eq. 12 ) then
943 #ifdef ALLOW_OBCSS_CONTROL
944 nvarlength = nvarlength +
945 & (ncvarrecs(i)/nobcs)
946 & *nwetobcss(bi,bj,k,iobcs)
947 #endif
948 else if ( i .eq. 13 ) then
949 #ifdef ALLOW_OBCSW_CONTROL
950 nvarlength = nvarlength +
951 & (ncvarrecs(i)/nobcs)
952 & *nwetobcsw(bi,bj,k,iobcs)
953 #endif
954 else if ( i .eq. 14 ) then
955 #ifdef ALLOW_OBCSE_CONTROL
956 nvarlength = nvarlength +
957 & (ncvarrecs(i)/nobcs)
958 & *nwetobcse(bi,bj,k,iobcs)
959 #endif
960 end if
961 enddo
962 #endif
963 else
964 print*,'ctrl_init: invalid grid location'
965 print*,' control variable = ',ncvarindex(i)
966 print*,' grid location = ',ncvargrd(i)
967 stop ' ... stopped in ctrl_init'
968 endif
969 enddo
970 enddo
971 enddo
972 endif
973 enddo
974
975 cph(
976 print *, 'ph-wet 1: nvarlength = ', nvarlength
977 print *, 'ph-wet 2: surface wet C = ', nwetctile(1,1,1)
978 print *, 'ph-wet 3: surface wet W = ', nwetwtile(1,1,1)
979 print *, 'ph-wet 4: surface wet S = ', nwetstile(1,1,1)
980 print *, 'ph-wet 4a:surface wet V = ', nwetvtile(1,1,1)
981 nwetc3d = 0
982 do k = 1, Nr
983 nwetc3d = nwetc3d + nwetctile(1,1,k)
984 end do
985 print *, 'ph-wet 5: 3D wet points = ', nwetc3d
986 do i = 1, maxcvars
987 print *, 'ph-wet 6: no recs for i = ', i, ncvarrecs(i)
988 end do
989 print *, 'ph-wet 7: ',
990 & 2*nwetc3d +
991 & ncvarrecs(3)*nwetctile(1,1,1) +
992 & ncvarrecs(4)*nwetctile(1,1,1) +
993 & ncvarrecs(5)*nwetwtile(1,1,1) +
994 & ncvarrecs(6)*nwetstile(1,1,1)
995 print *, 'ph-wet 8: ',
996 & 2*nwetc3d +
997 & ncvarrecs(7)*nwetctile(1,1,1) +
998 & ncvarrecs(8)*nwetctile(1,1,1) +
999 & ncvarrecs(9)*nwetwtile(1,1,1) +
1000 & ncvarrecs(10)*nwetstile(1,1,1)
1001 #ifdef ALLOW_OBCSN_CONTROL
1002 print *, 'ph-wet 9: surface wet obcsn = '
1003 & , nwetobcsn(1,1,1,1), nwetobcsn(1,1,1,2)
1004 & , nwetobcsn(1,1,1,3), nwetobcsn(1,1,1,4)
1005 #endif
1006 #ifdef ALLOW_OBCSS_CONTROL
1007 print *, 'ph-wet 10: surface wet obcss = '
1008 & , nwetobcss(1,1,1,1), nwetobcss(1,1,1,2)
1009 & , nwetobcss(1,1,1,3), nwetobcss(1,1,1,4)
1010 #endif
1011 #ifdef ALLOW_OBCSW_CONTROL
1012 print *, 'ph-wet 11: surface wet obcsw = '
1013 & , nwetobcsw(1,1,1,1), nwetobcsw(1,1,1,2)
1014 & , nwetobcsw(1,1,1,3), nwetobcsw(1,1,1,4)
1015 #endif
1016 #ifdef ALLOW_OBCSE_CONTROL
1017 print *, 'ph-wet 12: surface wet obcse = '
1018 & , nwetobcse(1,1,1,1), nwetobcse(1,1,1,2)
1019 & , nwetobcse(1,1,1,3), nwetobcse(1,1,1,4)
1020 #endif
1021 cph)
1022
1023 CALL GLOBAL_SUM_INT( nvarlength, myThid )
1024
1025 print *, 'ph-wet 13: global nvarlength vor k=', k, nvarlength
1026
1027 c
1028 c Summation of wet point counters
1029 c
1030 do k = 1, nr
1031
1032 ntmp=0
1033 do bj=1,nSy
1034 do bi=1,nSx
1035 ntmp=ntmp+nWetcTile(bi,bj,k)
1036 enddo
1037 enddo
1038 CALL GLOBAL_SUM_INT( ntmp, myThid )
1039 nWetcGlobal(k)=ntmp
1040 print *, 'ph-wet 14a: global nWet... k=', k, ntmp
1041
1042 ntmp=0
1043 do bj=1,nSy
1044 do bi=1,nSx
1045 ntmp=ntmp+nWetsTile(bi,bj,k)
1046 enddo
1047 enddo
1048 CALL GLOBAL_SUM_INT( ntmp, myThid )
1049 nWetsGlobal(k)=ntmp
1050 print *, 'ph-wet 14b: global nWet... k=', k, ntmp
1051
1052 ntmp=0
1053 do bj=1,nSy
1054 do bi=1,nSx
1055 ntmp=ntmp+nWetwTile(bi,bj,k)
1056 enddo
1057 enddo
1058 CALL GLOBAL_SUM_INT( ntmp, myThid )
1059 nWetwGlobal(k)=ntmp
1060 print *, 'ph-wet 14c: global nWet... k=', k, ntmp
1061
1062 ntmp=0
1063 do bj=1,nSy
1064 do bi=1,nSx
1065 ntmp=ntmp+nWetvTile(bi,bj,k)
1066 enddo
1067 enddo
1068 CALL GLOBAL_SUM_INT( ntmp, myThid )
1069 nWetvGlobal(k)=ntmp
1070 print *, 'ph-wet 14d: global nWet... k=', k, ntmp
1071
1072 #ifdef ALLOW_OBCSN_CONTROL
1073 do iobcs = 1, nobcs
1074 ntmp=0
1075 do bj=1,nSy
1076 do bi=1,nSx
1077 ntmp=ntmp+nwetobcsn(bi,bj,k,iobcs)
1078 enddo
1079 enddo
1080 CALL GLOBAL_SUM_INT( ntmp, myThid )
1081 nwetobcsnglo(k,iobcs)=ntmp
1082 print *, 'ph-wet 15a: global nWet... k=', k, iobcs, ntmp
1083 enddo
1084 #endif
1085 #ifdef ALLOW_OBCSS_CONTROL
1086 do iobcs = 1, nobcs
1087 ntmp=0
1088 do bj=1,nSy
1089 do bi=1,nSx
1090 ntmp=ntmp+nwetobcss(bi,bj,k,iobcs)
1091 enddo
1092 enddo
1093 CALL GLOBAL_SUM_INT( ntmp, myThid )
1094 nwetobcssglo(k,iobcs)=ntmp
1095 print *, 'ph-wet 15b: global nWet... k=', k, iobcs, ntmp
1096 enddo
1097 #endif
1098 #ifdef ALLOW_OBCSW_CONTROL
1099 do iobcs = 1, nobcs
1100 ntmp=0
1101 do bj=1,nSy
1102 do bi=1,nSx
1103 ntmp=ntmp+nwetobcsw(bi,bj,k,iobcs)
1104 enddo
1105 enddo
1106 CALL GLOBAL_SUM_INT( ntmp, myThid )
1107 nwetobcswglo(k,iobcs)=ntmp
1108 print *, 'ph-wet 15c: global nWet... k=', k, iobcs, ntmp
1109 enddo
1110 #endif
1111 #ifdef ALLOW_OBCSE_CONTROL
1112 do iobcs = 1, nobcs
1113 ntmp=0
1114 do bj=1,nSy
1115 do bi=1,nSx
1116 ntmp=ntmp+nwetobcse(bi,bj,k,iobcs)
1117 enddo
1118 enddo
1119 CALL GLOBAL_SUM_INT( ntmp, myThid )
1120 nwetobcseglo(k,iobcs)=ntmp
1121 print *, 'ph-wet 15d: global nWet... k=', k, iobcs, ntmp
1122 enddo
1123 #endif
1124
1125 enddo
1126
1127 print*, 'ctrl_init: no. of control variables: ', nvartype
1128 print*, 'ctrl_init: control vector length: ', nvarlength
1129 _END_MASTER( mythid )
1130
1131 c write masks and weights to files to be read by a master process
1132 c
1133 call active_write_xyz( 'hFacC', hFacC, 1, 0, mythid, dummy)
1134 call active_write_xyz( 'maskW', maskW, 1, 0, mythid, dummy)
1135 call active_write_xyz( 'maskS', maskS, 1, 0, mythid, dummy)
1136 #if (defined (ALLOW_EFLUXP0_CONTROL))
1137 call active_write_xyz( 'hFacV', hFacV, 1, 0, mythid, dummy)
1138 #endif
1139
1140 c-- Summarize the control vector's setup.
1141 _BEGIN_MASTER( mythid )
1142 cph call ctrl_Summary( mythid )
1143 _END_MASTER( mythid )
1144
1145 _BARRIER
1146
1147 return
1148 end
1149

  ViewVC Help
Powered by ViewVC 1.1.22