13 use mom_coms,
only : pe_here, root_pe, num_pes, sum_across_pes
22 implicit none ;
private
83 #include "version_variable.h"
84 character(len=40) :: mdl =
"MOM_debugging"
88 "If true, write out verbose debugging data.", &
89 default=.false., debuggingparam=.true.)
91 "If true, checksums are performed on arrays in the "//&
92 "various vec_chksum routines.", default=
debug, &
93 debuggingparam=.true.)
95 "If true, debug redundant data points during calls to "//&
96 "the various vec_chksum routines.", default=
debug, &
97 debuggingparam=.true.)
106 character(len=*),
intent(in) :: mesg
108 real,
dimension(G%IsdB:,G%jsd:,:),
intent(in) :: u_comp
110 real,
dimension(G%isd:,G%JsdB:,:),
intent(in) :: v_comp
112 integer,
optional,
intent(in) :: is
113 integer,
optional,
intent(in) :: ie
114 integer,
optional,
intent(in) :: js
115 integer,
optional,
intent(in) :: je
116 integer,
optional,
intent(in) :: direction
119 character(len=24) :: mesg_k
122 do k=1,
size(u_comp,3)
123 if (k < 10)
then ;
write(mesg_k,
'(" Layer",i2," ")') k
124 elseif (k < 100)
then ;
write(mesg_k,
'(" Layer",i3," ")') k
125 elseif (k < 1000)
then ;
write(mesg_k,
'(" Layer",i4," ")') k
126 else ;
write(mesg_k,
'(" Layer",i9," ")') k ;
endif
129 v_comp(:,:,k), g, is, ie, js, je, direction)
136 character(len=*),
intent(in) :: mesg
138 real,
dimension(G%IsdB:,G%jsd:),
intent(in) :: u_comp
140 real,
dimension(G%isd:,G%JsdB:),
intent(in) :: v_comp
142 integer,
optional,
intent(in) :: is
143 integer,
optional,
intent(in) :: ie
144 integer,
optional,
intent(in) :: js
145 integer,
optional,
intent(in) :: je
146 integer,
optional,
intent(in) :: direction
149 real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed)
150 real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed)
151 real :: u_resym(G%IsdB:G%IedB,G%jsd:G%jed)
152 real :: v_resym(G%isd:G%ied,G%JsdB:G%JedB)
153 character(len=128) :: mesg2
154 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
155 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
157 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
158 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
159 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
161 if (.not.(
present(is) .or.
present(ie) .or.
present(js) .or.
present(je)))
then
163 if ((isd == isdb) .and. (jsd == jsdb))
return
166 do i=isd,ied ;
do j=jsd,jed
167 u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
170 if (.not.
associated(g%Domain_aux))
call mom_error(fatal,
" check_redundant"//&
171 " called with a non-associated auxiliary domain the grid type.")
172 call pass_vector(u_nonsym, v_nonsym, g%Domain_aux, direction)
174 do i=isdb,iedb ;
do j=jsd,jed ; u_resym(i,j) = u_comp(i,j) ;
enddo ;
enddo
175 do i=isd,ied ;
do j=jsdb,jedb ; v_resym(i,j) = v_comp(i,j) ;
enddo ;
enddo
176 do i=isd,ied ;
do j=jsd,jed
177 u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j)
179 call pass_vector(u_resym, v_resym, g%Domain, direction)
181 is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
182 if (
present(is)) is_ch = is ;
if (
present(ie)) ie_ch = ie
183 if (
present(js)) js_ch = js ;
if (
present(js)) je_ch = je
185 do i=is_ch,ie_ch ;
do j=js_ch+1,je_ch
186 if (u_resym(i,j) /= u_comp(i,j) .and. &
188 write(mesg2,
'(" redundant u-components",2(1pe12.4)," differ by ", &
189 & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
190 u_comp(i,j), u_resym(i,j),u_comp(i,j)-u_resym(i,j),i,j,pe_here()
191 write(0,
'(A130)') trim(mesg)//trim(mesg2)
195 do i=is_ch+1,ie_ch ;
do j=js_ch,je_ch
196 if (v_resym(i,j) /= v_comp(i,j) .and. &
198 write(mesg2,
'(" redundant v-comps",2(1pe12.4)," differ by ", &
199 & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') &
200 v_comp(i,j), v_resym(i,j),v_comp(i,j)-v_resym(i,j),i,j, &
201 g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
202 write(0,
'(A155)') trim(mesg)//trim(mesg2)
211 character(len=*),
intent(in) :: mesg
213 real,
dimension(G%IsdB:,G%JsdB:,:),
intent(in) :: array
214 integer,
optional,
intent(in) :: is
215 integer,
optional,
intent(in) :: ie
216 integer,
optional,
intent(in) :: js
217 integer,
optional,
intent(in) :: je
220 character(len=24) :: mesg_k
224 if (k < 10)
then ;
write(mesg_k,
'(" Layer",i2," ")') k
225 elseif (k < 100)
then ;
write(mesg_k,
'(" Layer",i3," ")') k
226 elseif (k < 1000)
then ;
write(mesg_k,
'(" Layer",i4," ")') k
227 else ;
write(mesg_k,
'(" Layer",i9," ")') k ;
endif
236 character(len=*),
intent(in) :: mesg
238 real,
dimension(G%IsdB:,G%JsdB:),
intent(in) :: array
239 integer,
optional,
intent(in) :: is
240 integer,
optional,
intent(in) :: ie
241 integer,
optional,
intent(in) :: js
242 integer,
optional,
intent(in) :: je
244 real :: a_nonsym(G%isd:G%ied,G%jsd:G%jed)
245 real :: a_resym(G%IsdB:G%IedB,G%JsdB:G%JedB)
246 character(len=128) :: mesg2
247 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
248 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
250 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
251 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
252 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
254 if (.not.(
present(is) .or.
present(ie) .or.
present(js) .or.
present(je)))
then
256 if ((isd == isdb) .and. (jsd == jsdb))
return
259 do i=isd,ied ;
do j=jsd,jed
260 a_nonsym(i,j) = array(i,j)
263 if (.not.
associated(g%Domain_aux))
call mom_error(fatal,
" check_redundant"//&
264 " called with a non-associated auxiliary domain the grid type.")
265 call pass_vector(a_nonsym, a_nonsym, g%Domain_aux, &
266 direction=
to_all+scalar_pair, stagger=bgrid_ne)
268 do i=isdb,iedb ;
do j=jsdb,jedb ; a_resym(i,j) = array(i,j) ;
enddo ;
enddo
269 do i=isd,ied ;
do j=jsd,jed
270 a_resym(i,j) = a_nonsym(i,j)
275 is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
276 if (
present(is)) is_ch = is ;
if (
present(ie)) ie_ch = ie
277 if (
present(js)) js_ch = js ;
if (
present(js)) je_ch = je
279 do i=is_ch,ie_ch ;
do j=js_ch,je_ch
280 if (a_resym(i,j) /= array(i,j) .and. &
282 write(mesg2,
'(" Redundant points",2(1pe12.4)," differ by ", &
283 & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
284 array(i,j), a_resym(i,j),array(i,j)-a_resym(i,j),i,j,pe_here()
285 write(0,
'(A130)') trim(mesg)//trim(mesg2)
295 character(len=*),
intent(in) :: mesg
297 real,
dimension(G%IsdB:,G%JsdB:,:),
intent(in) :: u_comp
299 real,
dimension(G%IsdB:,G%JsdB:,:),
intent(in) :: v_comp
301 integer,
optional,
intent(in) :: is
302 integer,
optional,
intent(in) :: ie
303 integer,
optional,
intent(in) :: js
304 integer,
optional,
intent(in) :: je
305 integer,
optional,
intent(in) :: direction
308 character(len=24) :: mesg_k
311 do k=1,
size(u_comp,3)
312 if (k < 10)
then ;
write(mesg_k,
'(" Layer",i2," ")') k
313 elseif (k < 100)
then ;
write(mesg_k,
'(" Layer",i3," ")') k
314 elseif (k < 1000)
then ;
write(mesg_k,
'(" Layer",i4," ")') k
315 else ;
write(mesg_k,
'(" Layer",i9," ")') k ;
endif
318 v_comp(:,:,k), g, is, ie, js, je, direction)
325 character(len=*),
intent(in) :: mesg
327 real,
dimension(G%IsdB:,G%JsdB:),
intent(in) :: u_comp
329 real,
dimension(G%IsdB:,G%JsdB:),
intent(in) :: v_comp
331 integer,
optional,
intent(in) :: is
332 integer,
optional,
intent(in) :: ie
333 integer,
optional,
intent(in) :: js
334 integer,
optional,
intent(in) :: je
335 integer,
optional,
intent(in) :: direction
338 real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed)
339 real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed)
340 real :: u_resym(G%IsdB:G%IedB,G%JsdB:G%JedB)
341 real :: v_resym(G%IsdB:G%IedB,G%JsdB:G%JedB)
342 character(len=128) :: mesg2
343 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
344 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
346 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
347 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
348 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
350 if (.not.(
present(is) .or.
present(ie) .or.
present(js) .or.
present(je)))
then
352 if ((isd == isdb) .and. (jsd == jsdb))
return
355 do i=isd,ied ;
do j=jsd,jed
356 u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
359 if (.not.
associated(g%Domain_aux))
call mom_error(fatal,
" check_redundant"//&
360 " called with a non-associated auxiliary domain the grid type.")
361 call pass_vector(u_nonsym, v_nonsym, g%Domain_aux, direction, stagger=bgrid_ne)
363 do i=isdb,iedb ;
do j=jsdb,jedb
364 u_resym(i,j) = u_comp(i,j) ; v_resym(i,j) = v_comp(i,j)
366 do i=isd,ied ;
do j=jsd,jed
367 u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j)
369 call pass_vector(u_resym, v_resym, g%Domain, direction, stagger=bgrid_ne)
371 is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
372 if (
present(is)) is_ch = is ;
if (
present(ie)) ie_ch = ie
373 if (
present(js)) js_ch = js ;
if (
present(js)) je_ch = je
375 do i=is_ch,ie_ch ;
do j=js_ch,je_ch
376 if (u_resym(i,j) /= u_comp(i,j) .and. &
378 write(mesg2,
'(" redundant u-components",2(1pe12.4)," differ by ", &
379 & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
380 u_comp(i,j), u_resym(i,j),u_comp(i,j)-u_resym(i,j),i,j,pe_here()
381 write(0,
'(A130)') trim(mesg)//trim(mesg2)
385 do i=is_ch,ie_ch ;
do j=js_ch,je_ch
386 if (v_resym(i,j) /= v_comp(i,j) .and. &
388 write(mesg2,
'(" redundant v-comps",2(1pe12.4)," differ by ", &
389 & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') &
390 v_comp(i,j), v_resym(i,j),v_comp(i,j)-v_resym(i,j),i,j, &
391 g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
392 write(0,
'(A155)') trim(mesg)//trim(mesg2)
401 character(len=*),
intent(in) :: mesg
403 real,
dimension(G%isd:,G%jsd:,:),
intent(in) :: array
404 integer,
optional,
intent(in) :: is
405 integer,
optional,
intent(in) :: ie
406 integer,
optional,
intent(in) :: js
407 integer,
optional,
intent(in) :: je
409 character(len=24) :: mesg_k
413 if (k < 10)
then ;
write(mesg_k,
'(" Layer",i2," ")') k
414 elseif (k < 100)
then ;
write(mesg_k,
'(" Layer",i3," ")') k
415 elseif (k < 1000)
then ;
write(mesg_k,
'(" Layer",i4," ")') k
416 else ;
write(mesg_k,
'(" Layer",i9," ")') k ;
endif
426 character(len=*),
intent(in) :: mesg
428 real,
dimension(G%isd:,G%jsd:),
intent(in) :: array
429 integer,
optional,
intent(in) :: is
430 integer,
optional,
intent(in) :: ie
431 integer,
optional,
intent(in) :: js
432 integer,
optional,
intent(in) :: je
434 real :: a_nonsym(G%isd:G%ied,G%jsd:G%jed)
435 character(len=128) :: mesg2
437 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
438 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed
439 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
441 is_ch = g%isc ; ie_ch = g%iec ; js_ch = g%jsc ; je_ch = g%jec
442 if (
present(is)) is_ch = is ;
if (
present(ie)) ie_ch = ie
443 if (
present(js)) js_ch = js ;
if (
present(js)) je_ch = je
446 if ((is_ch == g%isc) .and. (ie_ch == g%iec) .and. &
447 (js_ch == g%jsc) .and. (je_ch == g%jec))
return
449 do i=isd,ied ;
do j=jsd,jed
450 a_nonsym(i,j) = array(i,j)
455 do i=is_ch,ie_ch ;
do j=js_ch,je_ch
456 if (a_nonsym(i,j) /= array(i,j) .and. &
458 write(mesg2,
'(" Redundant points",2(1pe12.4)," differ by ", &
459 & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
460 array(i,j), a_nonsym(i,j),array(i,j)-a_nonsym(i,j),i,j,pe_here()
461 write(0,
'(A130)') trim(mesg)//trim(mesg2)
471 character(len=*),
intent(in) :: mesg
473 real,
dimension(G%isd:,G%jsd:,:),
intent(in) :: u_comp
475 real,
dimension(G%isd:,G%jsd:,:),
intent(in) :: v_comp
477 integer,
optional,
intent(in) :: is
478 integer,
optional,
intent(in) :: ie
479 integer,
optional,
intent(in) :: js
480 integer,
optional,
intent(in) :: je
481 integer,
optional,
intent(in) :: direction
484 character(len=24) :: mesg_k
487 do k=1,
size(u_comp,3)
488 if (k < 10)
then ;
write(mesg_k,
'(" Layer",i2," ")') k
489 elseif (k < 100)
then ;
write(mesg_k,
'(" Layer",i3," ")') k
490 elseif (k < 1000)
then ;
write(mesg_k,
'(" Layer",i4," ")') k
491 else ;
write(mesg_k,
'(" Layer",i9," ")') k ;
endif
494 v_comp(:,:,k), g, is, ie, js, je, direction)
501 character(len=*),
intent(in) :: mesg
503 real,
dimension(G%isd:,G%jsd:),
intent(in) :: u_comp
505 real,
dimension(G%isd:,G%jsd:),
intent(in) :: v_comp
507 integer,
optional,
intent(in) :: is
508 integer,
optional,
intent(in) :: ie
509 integer,
optional,
intent(in) :: js
510 integer,
optional,
intent(in) :: je
511 integer,
optional,
intent(in) :: direction
514 real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed)
515 real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed)
516 character(len=128) :: mesg2
518 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
519 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
520 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
521 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
522 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
524 is_ch = g%isc ; ie_ch = g%iec ; js_ch = g%jsc ; je_ch = g%jec
525 if (
present(is)) is_ch = is ;
if (
present(ie)) ie_ch = ie
526 if (
present(js)) js_ch = js ;
if (
present(js)) je_ch = je
529 if ((is_ch == g%isc) .and. (ie_ch == g%iec) .and. &
530 (js_ch == g%jsc) .and. (je_ch == g%jec))
return
532 do i=isd,ied ;
do j=jsd,jed
533 u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
536 call pass_vector(u_nonsym, v_nonsym, g%Domain, direction, stagger=agrid)
538 do i=is_ch,ie_ch ;
do j=js_ch+1,je_ch
539 if (u_nonsym(i,j) /= u_comp(i,j) .and. &
541 write(mesg2,
'(" redundant u-components",2(1pe12.4)," differ by ", &
542 & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
543 u_comp(i,j), u_nonsym(i,j),u_comp(i,j)-u_nonsym(i,j),i,j,pe_here()
544 write(0,
'(A130)') trim(mesg)//trim(mesg2)
548 do i=is_ch+1,ie_ch ;
do j=js_ch,je_ch
549 if (v_nonsym(i,j) /= v_comp(i,j) .and. &
551 write(mesg2,
'(" redundant v-comps",2(1pe12.4)," differ by ", &
552 & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') &
553 v_comp(i,j), v_nonsym(i,j),v_comp(i,j)-v_nonsym(i,j),i,j, &
554 g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
555 write(0,
'(A155)') trim(mesg)//trim(mesg2)
563 subroutine chksum_vec_c3d(mesg, u_comp, v_comp, G, halos, scalars)
564 character(len=*),
intent(in) :: mesg
566 real,
dimension(G%IsdB:,G%jsd:,:),
intent(in) :: u_comp
567 real,
dimension(G%isd:,G%JsdB:,:),
intent(in) :: v_comp
568 integer,
optional,
intent(in) :: halos
569 logical,
optional,
intent(in) :: scalars
572 logical :: are_scalars
573 are_scalars = .false. ;
if (
present(scalars)) are_scalars = scalars
576 call uvchksum(mesg, u_comp, v_comp, g%HI, halos)
579 if (are_scalars)
then
589 subroutine chksum_vec_c2d(mesg, u_comp, v_comp, G, halos, scalars)
590 character(len=*),
intent(in) :: mesg
592 real,
dimension(G%IsdB:,G%jsd:),
intent(in) :: u_comp
593 real,
dimension(G%isd:,G%JsdB:),
intent(in) :: v_comp
594 integer,
optional,
intent(in) :: halos
595 logical,
optional,
intent(in) :: scalars
598 logical :: are_scalars
599 are_scalars = .false. ;
if (
present(scalars)) are_scalars = scalars
602 call uvchksum(mesg, u_comp, v_comp, g%HI, halos)
605 if (are_scalars)
then
615 subroutine chksum_vec_b3d(mesg, u_comp, v_comp, G, halos, scalars)
616 character(len=*),
intent(in) :: mesg
618 real,
dimension(G%IsdB:,G%JsdB:,:),
intent(in) :: u_comp
619 real,
dimension(G%IsdB:,G%JsdB:,:),
intent(in) :: v_comp
620 integer,
optional,
intent(in) :: halos
621 logical,
optional,
intent(in) :: scalars
624 logical :: are_scalars
625 are_scalars = .false. ;
if (
present(scalars)) are_scalars = scalars
628 call bchksum(u_comp, mesg//
"(u)", g%HI, halos)
629 call bchksum(v_comp, mesg//
"(v)", g%HI, halos)
632 if (are_scalars)
then
642 subroutine chksum_vec_b2d(mesg, u_comp, v_comp, G, halos, scalars, symmetric)
643 character(len=*),
intent(in) :: mesg
645 real,
dimension(G%IsdB:,G%JsdB:),
intent(in) :: u_comp
646 real,
dimension(G%IsdB:,G%JsdB:),
intent(in) :: v_comp
647 integer,
optional,
intent(in) :: halos
648 logical,
optional,
intent(in) :: scalars
650 logical,
optional,
intent(in) :: symmetric
653 logical :: are_scalars
654 are_scalars = .false. ;
if (
present(scalars)) are_scalars = scalars
657 call bchksum(u_comp, mesg//
"(u)", g%HI, halos, symmetric=symmetric)
658 call bchksum(v_comp, mesg//
"(v)", g%HI, halos, symmetric=symmetric)
661 if (are_scalars)
then
671 subroutine chksum_vec_a3d(mesg, u_comp, v_comp, G, halos, scalars)
672 character(len=*),
intent(in) :: mesg
674 real,
dimension(G%isd:,G%jsd:,:),
intent(in) :: u_comp
675 real,
dimension(G%isd:,G%jsd:,:),
intent(in) :: v_comp
676 integer,
optional,
intent(in) :: halos
677 logical,
optional,
intent(in) :: scalars
680 logical :: are_scalars
681 are_scalars = .false. ;
if (
present(scalars)) are_scalars = scalars
684 call hchksum(u_comp, mesg//
"(u)", g%HI, halos)
685 call hchksum(v_comp, mesg//
"(v)", g%HI, halos)
688 if (are_scalars)
then
698 subroutine chksum_vec_a2d(mesg, u_comp, v_comp, G, halos, scalars)
699 character(len=*),
intent(in) :: mesg
701 real,
dimension(G%isd:,G%jsd:),
intent(in) :: u_comp
702 real,
dimension(G%isd:,G%jsd:),
intent(in) :: v_comp
703 integer,
optional,
intent(in) :: halos
704 logical,
optional,
intent(in) :: scalars
707 logical :: are_scalars
708 are_scalars = .false. ;
if (
present(scalars)) are_scalars = scalars
711 call hchksum(u_comp, mesg//
"(u)", g%HI, halos)
712 call hchksum(v_comp, mesg//
"(v)", g%HI, halos)
715 if (are_scalars)
then
728 real,
dimension(HI%isd:,HI%jsd:,:),
intent(in) :: hthick
729 real,
dimension(HI%isd:,HI%jsd:),
intent(in) :: areat
730 real,
dimension(HI%isd:,HI%jsd:,:),
intent(in) :: stuff
733 integer :: i, j, k, nz
737 do k = 1, nz ;
do j = hi%jsc, hi%jec ;
do i = hi%isc, hi%iec
739 enddo ;
enddo ;
enddo
746 subroutine totaltands(HI, hThick, areaT, temperature, salinity, mesg)
748 real,
dimension(HI%isd:,HI%jsd:,:),
intent(in) :: hthick
749 real,
dimension(HI%isd:,HI%jsd:),
intent(in) :: areat
750 real,
dimension(HI%isd:,HI%jsd:,:),
intent(in) :: temperature
751 real,
dimension(HI%isd:,HI%jsd:,:),
intent(in) :: salinity
752 character(len=*),
intent(in) :: mesg
755 real,
save :: totalh = 0., totalt = 0., totals = 0.
757 logical,
save :: firstcall = .true.
758 real :: thish, thist, thiss, delh, delt, dels
759 integer :: i, j, k, nz
763 do k = 1, nz ;
do j = hi%jsc, hi%jec ;
do i = hi%isc, hi%iec
764 thish = thish + hthick(i,j,k) * areat(i,j)
765 enddo ;
enddo ;
enddo
766 call sum_across_pes(thish)
767 thist =
totalstuff(hi, hthick, areat, temperature)
768 thiss =
totalstuff(hi, hthick, areat, salinity)
772 totalh = thish ; totalt = thist ; totals = thiss
773 write(0,*)
'Totals H,T,S:',thish,thist,thiss,
' ',mesg
776 delh = thish - totalh
777 delt = thist - totalt
778 dels = thiss - totals
779 totalh = thish ; totalt = thist ; totals = thiss
780 write(0,*)
'Tot/del H,T,S:',thish,thist,thiss,delh,delt,dels,
' ',mesg
788 integer,
intent(in) :: nk
789 real,
dimension(nk),
intent(in) :: field
790 real,
optional,
intent(in) :: known_answer
793 real :: u_sum, error, expected
801 u_sum = u_sum + field(k)
802 error = error + epsilon(u_sum)*max(abs(u_sum),abs(field(k)))
806 if (
present(known_answer))
then
807 expected = known_answer
813 if (abs(u_sum-expected) > error)
then
823 integer,
intent(in) :: nk_1
824 integer,
intent(in) :: nk_2
825 real,
dimension(nk_1),
intent(in) :: field_1
826 real,
dimension(nk_2),
intent(in) :: field_2
827 real,
optional,
intent(in) :: missing_value
830 real :: u1_sum, error1, u2_sum, error2, misval
834 if (
present(missing_value))
then
835 misval = missing_value
845 if (field_1(k)/=misval)
then
846 u1_sum = u1_sum + field_1(k)
847 error1 = error1 + epsilon(u1_sum)*max(abs(u1_sum),abs(field_1(k)))
856 if (field_2(k)/=misval)
then
857 u2_sum = u2_sum + field_2(k)
858 error2 = error2 + epsilon(u2_sum)*max(abs(u2_sum),abs(field_2(k)))
863 if (abs(u1_sum-u2_sum) > (error1+error2))
then