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

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

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


Revision 1.26 - (show annotations) (download)
Tue Oct 30 20:19:13 2007 UTC (16 years, 7 months ago) by heimbach
Branch: MAIN
Changes since 1.25: +3 -1 lines
Add flag to disable theta threshold after ctrl update.

1 C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_map_ini.F,v 1.25 2007/10/09 00:00:00 jmc Exp $
2 C $Name: $
3
4 #include "CTRL_CPPOPTIONS.h"
5
6 CBOP
7 C !ROUTINE: ctrl_map_ini
8 C !INTERFACE:
9 subroutine ctrl_map_ini( mythid )
10
11 C !DESCRIPTION: \bv
12 c *=================================================================
13 c | SUBROUTINE ctrl_map_ini
14 c | Add the temperature, salinity, and diffusivity parts of the
15 c | control vector to the model state and update the tile halos.
16 c | The control vector is defined in the header file "ctrl.h".
17 c *=================================================================
18 C \ev
19
20 C !USES:
21 implicit none
22
23 c == global variables ==
24 #include "SIZE.h"
25 #include "EEPARAMS.h"
26 #include "PARAMS.h"
27 #include "GRID.h"
28 #include "DYNVARS.h"
29 #include "FFIELDS.h"
30 #include "ctrl.h"
31 #include "ctrl_dummy.h"
32 #include "optim.h"
33 #ifdef ALLOW_PTRACERS
34 # include "PTRACERS_SIZE.h"
35 # include "PTRACERS.h"
36 #endif
37 #ifdef ALLOW_ECCO
38 # include "ecco_cost.h"
39 #endif
40
41 C !INPUT/OUTPUT PARAMETERS:
42 c == routine arguments ==
43 integer mythid
44
45 C !LOCAL VARIABLES:
46 c == local variables ==
47
48 integer bi,bj
49 integer i,j,k
50 integer itlo,ithi
51 integer jtlo,jthi
52 integer jmin,jmax
53 integer imin,imax
54 integer il
55
56 logical equal
57 logical doglobalread
58 logical ladinit
59
60 character*( 80) fnamegeneric
61
62 _RL fac
63 _RL tmptest
64
65 c == external ==
66 integer ilnblnk
67 external ilnblnk
68
69 c == end of interface ==
70 CEOP
71
72 jtlo = mybylo(mythid)
73 jthi = mybyhi(mythid)
74 itlo = mybxlo(mythid)
75 ithi = mybxhi(mythid)
76 jmin = 1
77 jmax = sny
78 imin = 1
79 imax = snx
80
81 doglobalread = .false.
82 ladinit = .false.
83
84 equal = .true.
85
86 if ( equal ) then
87 fac = 1. _d 0
88 else
89 fac = 0. _d 0
90 endif
91
92 #ifdef ALLOW_THETA0_CONTROL
93 c-- Temperature field.
94 il=ilnblnk( xx_theta_file )
95 write(fnamegeneric(1:80),'(2a,i10.10)')
96 & xx_theta_file(1:il),'.',optimcycle
97 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
98 & doglobalread, ladinit, optimcycle,
99 & mythid, xx_theta_dummy )
100
101 do bj = jtlo,jthi
102 do bi = itlo,ithi
103 do k = 1,nr
104 do j = jmin,jmax
105 do i = imin,imax
106 #ifdef ALLOW_ECCO
107 IF (abs(tmpfld3d(i,j,k,bi,bj)).gt.
108 $ 2.0/sqrt(wtheta(k,bi,bj)))
109 $ tmpfld3d(i,j,k,bi,bj)=
110 $ sign(2.0/sqrt(wtheta(k,bi,bj)),tmpfld3d(i,j,k,bi,bj))
111 #endif
112 #ifdef ALLOW_AUTODIFF_OPENAD
113 theta(i,j,k,bi,bj) = theta(i,j,k,bi,bj) +
114 & fac*xx_theta(i,j,k,bi,bj) +
115 & fac*tmpfld3d(i,j,k,bi,bj)
116 #else
117 theta(i,j,k,bi,bj) = theta(i,j,k,bi,bj) +
118 & fac*tmpfld3d(i,j,k,bi,bj)
119 #endif
120 #ifndef DISABLE_CTRL_THETA_LIMIT
121 if(theta(i,j,k,bi,bj).lt.-2.0)
122 & theta(i,j,k,bi,bj)= -2.0
123 #endif
124 enddo
125 enddo
126 enddo
127 enddo
128 enddo
129
130 #endif
131
132 #ifdef ALLOW_SALT0_CONTROL
133 c-- Temperature field.
134 il=ilnblnk( xx_salt_file )
135 write(fnamegeneric(1:80),'(2a,i10.10)')
136 & xx_salt_file(1:il),'.',optimcycle
137 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
138 & doglobalread, ladinit, optimcycle,
139 & mythid, xx_salt_dummy )
140
141 do bj = jtlo,jthi
142 do bi = itlo,ithi
143 do k = 1,nr
144 do j = jmin,jmax
145 do i = imin,imax
146 #ifdef ALLOW_ECCO
147 IF (abs(tmpfld3d(i,j,k,bi,bj)).gt.
148 $ 2.0/sqrt(wsalt(k,bi,bj)))
149 $ tmpfld3d(i,j,k,bi,bj)=
150 $ sign(2.0/sqrt(wsalt(k,bi,bj)),tmpfld3d(i,j,k,bi,bj))
151 #endif
152 #ifdef ALLOW_AUTODIFF_OPENAD
153 salt(i,j,k,bi,bj) = salt(i,j,k,bi,bj) +
154 & fac*xx_salt(i,j,k,bi,bj) +
155 & fac*tmpfld3d(i,j,k,bi,bj)
156 #else
157 salt(i,j,k,bi,bj) = salt(i,j,k,bi,bj) +
158 & fac*tmpfld3d(i,j,k,bi,bj)
159 #endif
160
161 enddo
162 enddo
163 enddo
164 enddo
165 enddO
166 #endif
167
168 #ifdef ALLOW_TR10_CONTROL
169 #ifdef ALLOW_PTRACERS
170 c-- Temperature field.
171 il=ilnblnk( xx_tr1_file )
172 write(fnamegeneric(1:80),'(2a,i10.10)')
173 & xx_tr1_file(1:il),'.',optimcycle
174 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
175 & doglobalread, ladinit, optimcycle,
176 & mythid, xx_tr1_dummy )
177
178 do bj = jtlo,jthi
179 do bi = itlo,ithi
180 do k = 1,nr
181 do j = jmin,jmax
182 do i = imin,imax
183 ptracer(i,j,k,bi,bj,1) = ptracer(i,j,k,bi,bj,1) +
184 & fac*tmpfld3d(i,j,k,bi,bj)
185 enddo
186 enddo
187 enddo
188 enddo
189 enddo
190 #endif
191 #endif
192
193 #ifdef ALLOW_SST0_CONTROL
194 c-- sst0.
195 il=ilnblnk( xx_sst_file )
196 write(fnamegeneric(1:80),'(2a,i10.10)')
197 & xx_sst_file(1:il),'.',optimcycle
198 call active_read_xy ( fnamegeneric, tmpfld2d, 1,
199 & doglobalread, ladinit, optimcycle,
200 & mythid, xx_sst_dummy )
201 do bj = jtlo,jthi
202 do bi = itlo,ithi
203 do j = jmin,jmax
204 do i = imin,imax
205 cph sst(i,j,bi,bj) = sst(i,j,bi,bj) + tmpfld2d(i,j,bi,bj)
206 theta(i,j,1,bi,bj) = theta(i,j,1,bi,bj)
207 & + tmpfld2d(i,j,bi,bj)
208 enddo
209 enddo
210 enddo
211 enddo
212 #endif
213
214 #ifdef ALLOW_SSS0_CONTROL
215 c-- sss0.
216 il=ilnblnk( xx_sss_file )
217 write(fnamegeneric(1:80),'(2a,i10.10)')
218 & xx_sss_file(1:il),'.',optimcycle
219 call active_read_xy ( fnamegeneric, tmpfld2d, 1,
220 & doglobalread, ladinit, optimcycle,
221 & mythid, xx_sss_dummy )
222 do bj = jtlo,jthi
223 do bi = itlo,ithi
224 do j = jmin,jmax
225 do i = imin,imax
226 cph sss(i,j,bi,bj) = sss(i,j,bi,bj) + tmpfld2d(i,j,bi,bj)
227 salt(i,j,1,bi,bj) = salt(i,j,1,bi,bj)
228 & + tmpfld2d(i,j,bi,bj)
229 enddo
230 enddo
231 enddo
232 enddo
233 #endif
234
235 #ifdef ALLOW_DIFFKR_CONTROL
236 c-- diffkr.
237 il=ilnblnk( xx_diffkr_file )
238 write(fnamegeneric(1:80),'(2a,i10.10)')
239 & xx_diffkr_file(1:il),'.',optimcycle
240 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
241 & doglobalread, ladinit, optimcycle,
242 & mythid, xx_diffkr_dummy )
243 do bj = jtlo,jthi
244 do bi = itlo,ithi
245 do k = 1,nr
246 do j = jmin,jmax
247 do i = imin,imax
248 diffkr(i,j,k,bi,bj) = diffkr(i,j,k,bi,bj) +
249 & tmpfld3d(i,j,k,bi,bj)
250 enddo
251 enddo
252 enddo
253 enddo
254 enddo
255 #endif
256
257 #ifdef ALLOW_KAPGM_CONTROL
258 c-- kapgm.
259 il=ilnblnk( xx_kapgm_file )
260 write(fnamegeneric(1:80),'(2a,i10.10)')
261 & xx_kapgm_file(1:il),'.',optimcycle
262 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
263 & doglobalread, ladinit, optimcycle,
264 & mythid, xx_kapgm_dummy )
265 do bj = jtlo,jthi
266 do bi = itlo,ithi
267 do k = 1,nr
268 do j = jmin,jmax
269 do i = imin,imax
270 kapgm(i,j,k,bi,bj) = kapgm(i,j,k,bi,bj) +
271 & tmpfld3d(i,j,k,bi,bj)
272 enddo
273 enddo
274 enddo
275 enddo
276 enddo
277 #endif
278
279 #ifdef ALLOW_EFLUXY0_CONTROL
280 c-- y-component EP-flux field.
281 il=ilnblnk( xx_efluxy_file )
282 write(fnamegeneric(1:80),'(2a,i10.10)')
283 & xx_efluxy_file(1:il),'.',optimcycle
284 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
285 & doglobalread, ladinit, optimcycle,
286 & mythid, xx_efluxy_dummy )
287
288 do bj = jtlo,jthi
289 do bi = itlo,ithi
290 do k = 1,nr
291 do j = jmin,jmax
292 do i = imin,imax
293 EfluxY(i,j,k,bi,bj) = EfluxY(i,j,k,bi,bj)
294 & - fac*tmpfld3d(i,j,k,bi,bj)
295 & *maskS(i,j,k,bi,bj)
296 cph EfluxY(i,j,k,bi,bj) = EfluxY(i,j,k,bi,bj)
297 cph & - rSphere*cosFacU(J,bi,bj)
298 cph & *fac*tmpfld3d(i,j,k,bi,bj)
299 enddo
300 enddo
301 enddo
302 enddo
303 enddo
304 #endif
305
306 #ifdef ALLOW_EFLUXP0_CONTROL
307 c-- p-component EP-flux field.
308 il=ilnblnk( xx_efluxp_file )
309 write(fnamegeneric(1:80),'(2a,i10.10)')
310 & xx_efluxp_file(1:il),'.',optimcycle
311 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
312 & doglobalread, ladinit, optimcycle,
313 & mythid, xx_efluxp_dummy )
314
315 do bj = jtlo,jthi
316 do bi = itlo,ithi
317 do k = 1,nr
318 do j = jmin,jmax
319 do i = imin,imax
320 EfluxP(i,j,k,bi,bj) = EfluxP(i,j,k,bi,bj)
321 & + fCori(i,j,bi,bj)
322 & *fac*tmpfld3d(i,j,k,bi,bj)
323 & *hFacV(i,j,k,bi,bj)
324 cph EfluxP(i,j,k,bi,bj) = EfluxP(i,j,k,bi,bj)
325 cph & + fCori(i,j,bi,bj)
326 cph & *rSphere*cosFacU(J,bi,bj)
327 cph & *fac*tmpfld3d(i,j,k,bi,bj)
328 enddo
329 enddo
330 enddo
331 enddo
332 enddo
333 #endif
334
335 #ifdef ALLOW_BOTTOMDRAG_CONTROL
336 c-- bottom drag
337 il=ilnblnk( xx_bottomdrag_file )
338 write(fnamegeneric(1:80),'(2a,i10.10)')
339 & xx_bottomdrag_file(1:il),'.',optimcycle
340 call active_read_xy ( fnamegeneric, tmpfld2d, 1,
341 & doglobalread, ladinit, optimcycle,
342 & mythid, xx_bottomdrag_dummy )
343 do bj = jtlo,jthi
344 do bi = itlo,ithi
345 do j = jmin,jmax
346 do i = imin,imax
347 bottomdragfld(i,j,bi,bj) = bottomdragfld(i,j,bi,bj)
348 & + tmpfld2d(i,j,bi,bj)
349 enddo
350 enddo
351 enddo
352 enddo
353 #endif
354
355 #ifdef ALLOW_EDTAUX_CONTROL
356 c-- zonal eddy stress : edtaux
357 il=ilnblnk( xx_edtaux_file )
358 write(fnamegeneric(1:80),'(2a,i10.10)')
359 & xx_edtaux_file(1:il),'.',optimcycle
360 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
361 & doglobalread, ladinit, optimcycle,
362 & mythid, xx_edtaux_dummy )
363 do bj = jtlo,jthi
364 do bi = itlo,ithi
365 do k = 1,nr
366 do j = jmin,jmax
367 do i = imin,imax
368 eddyTauX(i,j,k,bi,bj) = eddyTauX(i,j,k,bi,bj) +
369 & fCori(i,j,bi,bj)*tmpfld3d(i,j,k,bi,bj)
370 enddo
371 enddo
372 enddo
373 enddo
374 enddo
375 #endif
376
377 #ifdef ALLOW_EDTAUY_CONTROL
378 c-- meridional eddy stress : edtauy
379 il=ilnblnk( xx_edtauy_file )
380 write(fnamegeneric(1:80),'(2a,i10.10)')
381 & xx_edtauy_file(1:il),'.',optimcycle
382 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
383 & doglobalread, ladinit, optimcycle,
384 & mythid, xx_edtauy_dummy )
385 do bj = jtlo,jthi
386 do bi = itlo,ithi
387 do k = 1,nr
388 do j = jmin,jmax
389 do i = imin,imax
390 eddyTauY(i,j,k,bi,bj) = eddyTauY(i,j,k,bi,bj) +
391 & fCoriG(i,j,bi,bj)*tmpfld3d(i,j,k,bi,bj)
392 enddo
393 enddo
394 enddo
395 enddo
396 enddo
397 #endif
398
399 #ifdef ALLOW_UVEL0_CONTROL
400 c-- initial zonal velocity
401 il=ilnblnk( xx_uvel_file )
402 write(fnamegeneric(1:80),'(2a,i10.10)')
403 & xx_uvel_file(1:il),'.',optimcycle
404 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
405 & doglobalread, ladinit, optimcycle,
406 & mythid, xx_uvel_dummy )
407 do bj = jtlo,jthi
408 do bi = itlo,ithi
409 do k = 1,nr
410 do j = jmin,jmax
411 do i = imin,imax
412 #ifdef ALLOW_AUTODIFF_OPENAD
413 uVel(i,j,k,bi,bj) = uVel(i,j,k,bi,bj) +
414 & fac*xx_uvel(i,j,k,bi,bj)
415 #else
416 uVel(i,j,k,bi,bj) = uVel(i,j,k,bi,bj) +
417 & fac*tmpfld3d(i,j,k,bi,bj)
418 #endif
419 enddo
420 enddo
421 enddo
422 enddo
423 enddo
424 #endif
425
426 #ifdef ALLOW_VVEL0_CONTROL
427 c-- initial merid. velocity
428 il=ilnblnk( xx_vvel_file )
429 write(fnamegeneric(1:80),'(2a,i10.10)')
430 & xx_vvel_file(1:il),'.',optimcycle
431 call active_read_xyz( fnamegeneric, tmpfld3d, 1,
432 & doglobalread, ladinit, optimcycle,
433 & mythid, xx_vvel_dummy )
434 do bj = jtlo,jthi
435 do bi = itlo,ithi
436 do k = 1,nr
437 do j = jmin,jmax
438 do i = imin,imax
439 #ifdef ALLOW_AUTODIFF_OPENAD
440 vVel(i,j,k,bi,bj) = vVel(i,j,k,bi,bj) +
441 & fac*xx_vvel(i,j,k,bi,bj)
442 #else
443 vVel(i,j,k,bi,bj) = vVel(i,j,k,bi,bj) +
444 & fac*tmpfld3d(i,j,k,bi,bj)
445 #endif
446 enddo
447 enddo
448 enddo
449 enddo
450 enddo
451 #endif
452
453 #ifdef ALLOW_ETAN0_CONTROL
454 c-- initial Eta.
455 il=ilnblnk( xx_etan_file )
456 write(fnamegeneric(1:80),'(2a,i10.10)')
457 & xx_etan_file(1:il),'.',optimcycle
458 call active_read_xy ( fnamegeneric, tmpfld2d, 1,
459 & doglobalread, ladinit, optimcycle,
460 & mythid, xx_etan_dummy )
461 do bj = jtlo,jthi
462 do bi = itlo,ithi
463 do j = jmin,jmax
464 do i = imin,imax
465 #ifdef ALLOW_AUTODIFF_OPENAD
466 etaN(i,j,bi,bj) = etaN(i,j,bi,bj) +
467 & fac*xx_etan(i,j,bi,bj)
468 #else
469 etaN(i,j,bi,bj) = etaN(i,j,bi,bj) +
470 & fac*tmpfld2d(i,j,bi,bj)
471 #endif
472 enddo
473 enddo
474 enddo
475 enddo
476 #endif
477
478 #ifdef ALLOW_RELAXSST_CONTROL
479 c-- SST relaxation coefficient.
480 il=ilnblnk( xx_relaxsst_file )
481 write(fnamegeneric(1:80),'(2a,i10.10)')
482 & xx_relaxsst_file(1:il),'.',optimcycle
483 call active_read_xy ( fnamegeneric, tmpfld2d, 1,
484 & doglobalread, ladinit, optimcycle,
485 & mythid, xx_relaxsst_dummy )
486 do bj = jtlo,jthi
487 do bi = itlo,ithi
488 do j = jmin,jmax
489 do i = imin,imax
490 lambdaThetaClimRelax(i,j,bi,bj) =
491 & lambdaThetaClimRelax(i,j,bi,bj)
492 & + tmpfld2d(i,j,bi,bj)
493 enddo
494 enddo
495 enddo
496 enddo
497 #endif
498
499 #ifdef ALLOW_RELAXSSS_CONTROL
500 c-- SSS relaxation coefficient.
501 il=ilnblnk( xx_relaxsss_file )
502 write(fnamegeneric(1:80),'(2a,i10.10)')
503 & xx_relaxsss_file(1:il),'.',optimcycle
504 call active_read_xy ( fnamegeneric, tmpfld2d, 1,
505 & doglobalread, ladinit, optimcycle,
506 & mythid, xx_relaxsss_dummy )
507 do bj = jtlo,jthi
508 do bi = itlo,ithi
509 do j = jmin,jmax
510 do i = imin,imax
511 lambdaSaltClimRelax(i,j,bi,bj) =
512 & lambdaSaltClimRelax(i,j,bi,bj)
513 & + tmpfld2d(i,j,bi,bj)
514 enddo
515 enddo
516 enddo
517 enddo
518 #endif
519
520 #ifdef ALLOW_SEAICE
521 call seaice_ctrl_map_ini( mythid )
522 #endif
523
524 c-- Update the tile edges.
525
526 #if (defined (ALLOW_THETA0_CONTROL) || defined (ALLOW_SST0_CONTROL))
527 _EXCH_XYZ_R8( theta, mythid )
528 #endif
529 #if (defined (ALLOW_SALT0_CONTROL) || defined (ALLOW_SSS0_CONTROL))
530 _EXCH_XYZ_R8( salt, mythid )
531 #endif
532 #ifdef ALLOW_TR10_CONTROL
533 #ifdef ALLOW_PTRACERS
534 _EXCH_XYZ_R8(pTracer(1-Olx,1-Oly,1,1,1,1),myThid)
535 #endif
536 #endif
537 #ifdef ALLOW_DIFFKR_CONTROL
538 _EXCH_XYZ_R8( diffkr, mythid)
539 #endif
540 #ifdef ALLOW_KAPGM_CONTROL
541 _EXCH_XYZ_R8( kapgm, mythid)
542 #endif
543 #ifdef ALLOW_EFLUXY0_CONTROL
544 _EXCH_XYZ_R8( EfluxY, mythid )
545 #endif
546 #ifdef ALLOW_EFLUXP0_CONTROL
547 _EXCH_XYZ_R8( EfluxP, mythid )
548 #endif
549 #ifdef ALLOW_BOTTOMDRAG_CONTROL
550 _EXCH_XY_R8( bottomdragfld, mythid )
551 #endif
552
553 #if (defined (ALLOW_EDTAUX_CONTROL) && defined (ALLOW_EDTAUY_CONTROL))
554 CALL EXCH_UV_XYZ_RS(eddyTauX,eddyTauY,.TRUE.,myThid)
555 #elif (defined (ALLOW_EDTAUX_CONTROL) || defined (ALLOW_EDTAUY_CONTROL))
556 STOP 'ctrl_map_forcing: need BOTH ALLOW_EDTAU[X,Y]_CONTROL'
557 #endif
558
559 #ifdef ALLOW_UVEL0_CONTROL
560 _EXCH_XYZ_R8( uVel, mythid)
561 #endif
562
563 #ifdef ALLOW_VVEL0_CONTROL
564 _EXCH_XYZ_R8( vVel, mythid)
565 #endif
566
567 #ifdef ALLOW_ETAN0_CONTROL
568 _EXCH_XY_R8( etaN, mythid )
569 #endif
570
571 #ifdef ALLOW_RELAXSST_CONTROL
572 _EXCH_XY_R4( lambdaThetaClimRelax, mythid )
573 #endif
574
575 #ifdef ALLOW_RELAXSSS_CONTROL
576 _EXCH_XY_R4( lambdaThetaClimRelax, mythid )
577 #endif
578
579 return
580 end
581

  ViewVC Help
Powered by ViewVC 1.1.22