/[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.9 - (show annotations) (download)
Tue Jun 24 16:07:06 2003 UTC (20 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint51l_post, checkpoint51, checkpoint51f_post, checkpoint51d_post, checkpoint51j_post, checkpoint51l_pre, checkpoint51b_pre, checkpoint51h_pre, branchpoint-genmake2, checkpoint51i_post, checkpoint51b_post, checkpoint51c_post, checkpoint51i_pre, checkpoint51e_post, checkpoint51f_pre, checkpoint51g_post, checkpoint51m_post, checkpoint51a_post
Branch point for: branch-genmake2, tg2-branch
Changes since 1.8: +154 -90 lines
Merging for c51 vs. e34

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

  ViewVC Help
Powered by ViewVC 1.1.22