MOM_checksum_packages.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Provides routines that do checksums of groups of MOM variables
7
8! This module provides several routines that do check-sums of groups
9! of variables in the various dynamic solver routines.
10
11use mom_coms, only : min_across_pes, max_across_pes, reproducing_sum
12use mom_debugging, only : hchksum, uvchksum
13use mom_error_handler, only : mom_mesg, is_root_pe
14use mom_grid, only : ocean_grid_type
18
19implicit none ; private
20
21public mom_state_chksum, mom_thermo_chksum, mom_accel_chksum
23
24!> Write out checksums of the MOM6 state variables
25interface mom_state_chksum
26 module procedure mom_state_chksum_5arg
27 module procedure mom_state_chksum_3arg
28end interface
29
30#include <MOM_memory.h>
31
32!> A type for storing statistica about a variable
33type :: stats ; private
34 real :: minimum = 1.e34 !< The minimum value [degC] or [ppt] or other units
35 real :: maximum = -1.e34 !< The maximum value [degC] or [ppt] or other units
36 real :: average = 0. !< The average value [degC] or [ppt] or other units
37end type stats
38
39contains
40
41! =============================================================================
42
43!> Write out chksums for the model's basic state variables, including transports.
44subroutine mom_state_chksum_5arg(mesg, u, v, h, uh, vh, G, GV, US, haloshift, symmetric, omit_corners, vel_scale)
45 character(len=*), &
46 intent(in) :: mesg !< A message that appears on the chksum lines.
47 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
48 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
49 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
50 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1] or other units.
51 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
52 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1] or other units.
53 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
54 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
55 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
56 intent(in) :: uh !< Volume flux through zonal faces = u*h*dy
57 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
58 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
59 intent(in) :: vh !< Volume flux through meridional faces = v*h*dx
60 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
61 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
62 integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0).
63 logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully symmetric
64 !! computational domain.
65 logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
66 real, optional, intent(in) :: vel_scale !< The scaling factor to convert velocities to [T m L-1 s-1 ~> 1]
67
68 real :: scale_vel ! The scaling factor to convert velocities to mks units [T m L-1 s-1 ~> 1]
69 logical :: sym
70 integer :: hs
71
72 ! Note that for the chksum calls to be useful for reproducing across PE
73 ! counts, there must be no redundant points, so all variables use is..ie
74 ! and js...je as their extent.
75 hs = 1 ; if (present(haloshift)) hs=haloshift
76 sym = .false. ; if (present(symmetric)) sym=symmetric
77 scale_vel = us%L_T_to_m_s ; if (present(vel_scale)) scale_vel = vel_scale
78
79 call uvchksum(mesg//" [uv]", u, v, g%HI, haloshift=hs, symmetric=sym, &
80 omit_corners=omit_corners, unscale=scale_vel)
81 call hchksum(h, mesg//" h", g%HI, haloshift=hs, omit_corners=omit_corners, unscale=gv%H_to_MKS)
82 call uvchksum(mesg//" [uv]h", uh, vh, g%HI, haloshift=hs, symmetric=sym, &
83 omit_corners=omit_corners, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
84end subroutine mom_state_chksum_5arg
85
86! =============================================================================
87
88!> Write out chksums for the model's basic state variables.
89subroutine mom_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric, omit_corners)
90 character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines.
91 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
92 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
93 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
94 intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1] or [m s-1].
95 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
96 intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1] or [m s-1]..
97 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
98 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
99 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type, which is
100 !! used to rescale u and v if present.
101 integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0).
102 logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully
103 !! symmetric computational domain.
104 logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
105
106 integer :: hs
107 logical :: sym
108
109 ! Note that for the chksum calls to be useful for reproducing across PE
110 ! counts, there must be no redundant points, so all variables use is..ie
111 ! and js...je as their extent.
112 hs = 1 ; if (present(haloshift)) hs = haloshift
113 sym = .false. ; if (present(symmetric)) sym = symmetric
114 call uvchksum(mesg//" u", u, v, g%HI, haloshift=hs, symmetric=sym, &
115 omit_corners=omit_corners, unscale=us%L_T_to_m_s)
116 call hchksum(h, mesg//" h",g%HI, haloshift=hs, omit_corners=omit_corners, unscale=gv%H_to_MKS)
117end subroutine mom_state_chksum_3arg
118
119! =============================================================================
120
121!> Write out chksums for the model's thermodynamic state variables.
122subroutine mom_thermo_chksum(mesg, tv, G, US, haloshift, omit_corners)
123 character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines.
124 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
125 !! thermodynamic variables.
126 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
127 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
128 integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0).
129 logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
130
131 integer :: hs
132 hs=1 ; if (present(haloshift)) hs=haloshift
133
134 if (associated(tv%T)) &
135 call hchksum(tv%T, mesg//" T", g%HI, haloshift=hs, omit_corners=omit_corners, unscale=us%C_to_degC)
136 if (associated(tv%S)) &
137 call hchksum(tv%S, mesg//" S", g%HI, haloshift=hs, omit_corners=omit_corners, unscale=us%S_to_ppt)
138 if (associated(tv%frazil)) &
139 call hchksum(tv%frazil, mesg//" frazil", g%HI, haloshift=hs, omit_corners=omit_corners, &
140 unscale=us%Q_to_J_kg*us%R_to_kg_m3*us%Z_to_m)
141 if (associated(tv%salt_deficit)) &
142 call hchksum(tv%salt_deficit, mesg//" salt deficit", g%HI, haloshift=hs, omit_corners=omit_corners, &
143 unscale=us%S_to_ppt*us%RZ_to_kg_m2)
144 if (associated(tv%varT)) &
145 call hchksum(tv%varT, mesg//" varT", g%HI, haloshift=hs, omit_corners=omit_corners, unscale=us%C_to_degC**2)
146 if (associated(tv%varS)) &
147 call hchksum(tv%varS, mesg//" varS", g%HI, haloshift=hs, omit_corners=omit_corners, unscale=us%S_to_ppt**2)
148 if (associated(tv%covarTS)) &
149 call hchksum(tv%covarTS, mesg//" covarTS", g%HI, haloshift=hs, omit_corners=omit_corners, &
150 unscale=us%S_to_ppt*us%C_to_degC)
151
152end subroutine mom_thermo_chksum
153
154! =============================================================================
155
156!> Write out chksums for the ocean surface variables.
157subroutine mom_surface_chksum(mesg, sfc_state, G, US, haloshift, symmetric)
158 character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines.
159 type(surface), intent(inout) :: sfc_state !< transparent ocean surface state structure
160 !! shared with the calling routine data in this
161 !! structure is intent out.
162 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
163 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
164 integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0).
165 logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully symmetric
166 !! computational domain.
167
168 integer :: hs
169 logical :: sym
170
171 sym = .false. ; if (present(symmetric)) sym = symmetric
172 hs = 0 ; if (present(haloshift)) hs = haloshift
173
174 if (allocated(sfc_state%SST)) call hchksum(sfc_state%SST, mesg//" SST", g%HI, haloshift=hs, &
175 unscale=us%C_to_degC)
176 if (allocated(sfc_state%SSS)) call hchksum(sfc_state%SSS, mesg//" SSS", g%HI, haloshift=hs, &
177 unscale=us%S_to_ppt)
178 if (allocated(sfc_state%sea_lev)) call hchksum(sfc_state%sea_lev, mesg//" sea_lev", g%HI, &
179 haloshift=hs, unscale=us%Z_to_m)
180 if (allocated(sfc_state%Hml)) call hchksum(sfc_state%Hml, mesg//" Hml", g%HI, haloshift=hs, &
181 unscale=us%Z_to_m)
182 if (allocated(sfc_state%u) .and. allocated(sfc_state%v)) &
183 call uvchksum(mesg//" SSU", sfc_state%u, sfc_state%v, g%HI, haloshift=hs, symmetric=sym, &
184 unscale=us%L_T_to_m_s)
185 if (allocated(sfc_state%frazil)) call hchksum(sfc_state%frazil, mesg//" frazil", g%HI, &
186 haloshift=hs, unscale=us%Q_to_J_kg*us%RZ_to_kg_m2)
187 if (allocated(sfc_state%melt_potential)) call hchksum(sfc_state%melt_potential, mesg//" melt_potential", &
188 g%HI, haloshift=hs, unscale=us%Q_to_J_kg*us%RZ_to_kg_m2)
189 if (allocated(sfc_state%ocean_mass)) call hchksum(sfc_state%ocean_mass, mesg//" ocean_mass", &
190 g%HI, haloshift=hs, unscale=us%RZ_to_kg_m2)
191 if (allocated(sfc_state%ocean_heat)) call hchksum(sfc_state%ocean_heat, mesg//" ocean_heat", &
192 g%HI, haloshift=hs, unscale=us%C_to_degC*us%RZ_to_kg_m2)
193 if (allocated(sfc_state%ocean_salt)) call hchksum(sfc_state%ocean_salt, mesg//" ocean_salt", &
194 g%HI, haloshift=hs, unscale=us%S_to_ppt*us%RZ_to_kg_m2)
195
196end subroutine mom_surface_chksum
197
198! =============================================================================
199
200!> Write out chksums for the model's accelerations
201subroutine mom_accel_chksum(mesg, CAu, CAv, PFu, PFv, diffu, diffv, G, GV, US, pbce, &
202 u_accel_bt, v_accel_bt, symmetric)
203 character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines.
204 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
205 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
206 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
207 intent(in) :: cau !< Zonal acceleration due to Coriolis
208 !! and momentum advection terms [L T-2 ~> m s-2].
209 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
210 intent(in) :: cav !< Meridional acceleration due to Coriolis
211 !! and momentum advection terms [L T-2 ~> m s-2].
212 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
213 intent(in) :: pfu !< Zonal acceleration due to pressure gradients
214 !! (equal to -dM/dx) [L T-2 ~> m s-2].
215 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
216 intent(in) :: pfv !< Meridional acceleration due to pressure gradients
217 !! (equal to -dM/dy) [L T-2 ~> m s-2].
218 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
219 intent(in) :: diffu !< Zonal acceleration due to convergence of the
220 !! along-isopycnal stress tensor [L T-2 ~> m s-2].
221 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
222 intent(in) :: diffv !< Meridional acceleration due to convergence of
223 !! the along-isopycnal stress tensor [L T-2 ~> m s-2].
224 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
225 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
226 optional, intent(in) :: pbce !< The baroclinic pressure anomaly in each layer
227 !! due to free surface height anomalies
228 !! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
229 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
230 optional, intent(in) :: u_accel_bt !< The zonal acceleration from terms in the
231 !! barotropic solver [L T-2 ~> m s-2].
232 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
233 optional, intent(in) :: v_accel_bt !< The meridional acceleration from terms in
234 !! the barotropic solver [L T-2 ~> m s-2].
235 logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully symmetric
236 !! computational domain.
237
238 logical :: sym
239
240 sym = .false. ; if (present(symmetric)) sym = symmetric
241
242 ! Note that for the chksum calls to be useful for reproducing across PE
243 ! counts, there must be no redundant points, so all variables use is..ie
244 ! and js...je as their extent.
245 call uvchksum(mesg//" CA[uv]", cau, cav, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T2_to_m_s2)
246 call uvchksum(mesg//" PF[uv]", pfu, pfv, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T2_to_m_s2)
247 call uvchksum(mesg//" diffu", diffu, diffv, g%HI,haloshift=0, symmetric=sym, unscale=us%L_T2_to_m_s2)
248 if (present(pbce)) &
249 call hchksum(pbce, mesg//" pbce",g%HI,haloshift=0, unscale=gv%m_to_H*us%L_T_to_m_s**2)
250 if (present(u_accel_bt) .and. present(v_accel_bt)) &
251 call uvchksum(mesg//" [uv]_accel_bt", u_accel_bt, v_accel_bt, g%HI,haloshift=0, symmetric=sym, &
252 unscale=us%L_T2_to_m_s2)
253end subroutine mom_accel_chksum
254
255! =============================================================================
256
257!> Monitor and write out statistics for the model's state variables.
258subroutine mom_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, permitDiminishing)
259 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
260 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
261 character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines.
262 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
263 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
264 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
265 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
266 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
267 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
268 real, pointer, dimension(:,:,:), &
269 intent(in) :: temp !< Temperature [C ~> degC].
270 real, pointer, dimension(:,:,:), &
271 intent(in) :: salt !< Salinity [S ~> ppt].
272 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
273 logical, optional, intent(in) :: allowchange !< do not flag an error
274 !! if the statistics change.
275 logical, optional, intent(in) :: permitdiminishing !< do not flag error if the
276 !! extrema are diminishing.
277
278 ! Local variables
279 real, dimension(G%isc:G%iec, G%jsc:G%jec) :: &
280 tmp_a, & ! The area per cell [L2 ~> m2]
281 tmp_v, & ! The column-integrated volume or mass [H L2 ~> m3 or kg],
282 ! depending on whether the Boussinesq approximation is used
283 tmp_t, & ! The column-integrated temperature [C H L2 ~> degC m3 or degC kg]
284 tmp_s ! The column-integrated salinity [S H L2 ~> ppt m3 or ppt kg]
285 real :: vol, dv ! The total ocean volume or mass and its change [H L2 ~> m3 or kg]
286 real :: area ! The total ocean surface area [L2 ~> m2].
287 real :: h_minimum ! The minimum layer thicknesses [H ~> m or kg m-2]
288 real :: t_scale ! The scaling conversion factor for temperatures [degC C-1 ~> 1]
289 real :: s_scale ! The scaling conversion factor for salinities [ppt S-1 ~> 1]
290 logical :: do_ts ! If true, evaluate statistics for temperature and salinity
291 type(stats) :: t, delt ! Temperature statistics in unscaled units [degC]
292 type(stats) :: s, dels ! Salinity statistics in unscaled units [ppt]
293
294 ! NOTE: save data is not normally allowed but we use it for debugging purposes here on the
295 ! assumption we will not turn this on with threads
296 type(stats), save :: oldt, olds
297 logical, save :: firstcall = .true.
298 real, save :: oldvol ! The previous total ocean volume or mass [H L2 ~> m3 or kg]
299
300 character(len=80) :: lmsg
301 integer :: is, ie, js, je, nz, i, j, k
302
303 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
304 do_ts = associated(temp) .and. associated(salt)
305
306 tmp_a(:,:) = 0.0
307 tmp_v(:,:) = 0.0
308 tmp_t(:,:) = 0.0
309 tmp_s(:,:) = 0.0
310
311 t_scale = us%C_to_degC ; s_scale = us%S_to_ppt
312
313 ! First collect local stats
314 do j=js,je ; do i=is,ie
315 tmp_a(i,j) = tmp_a(i,j) + g%areaT(i,j)
316 enddo ; enddo
317 t%minimum = 1.e34 ; t%maximum = -1.e34 ; t%average = 0.
318 s%minimum = 1.e34 ; s%maximum = -1.e34 ; s%average = 0.
319 h_minimum = 1.e34*gv%m_to_H
320 do k=1,nz ; do j=js,je ; do i=is,ie
321 if (g%mask2dT(i,j)>0.) then
322 dv = g%areaT(i,j)*h(i,j,k)
323 tmp_v(i,j) = tmp_v(i,j) + dv
324 if (do_ts .and. h(i,j,k)>0.) then
325 t%minimum = min( t%minimum, t_scale*temp(i,j,k) ) ; t%maximum = max( t%maximum, t_scale*temp(i,j,k) )
326 s%minimum = min( s%minimum, s_scale*salt(i,j,k) ) ; s%maximum = max( s%maximum, s_scale*salt(i,j,k) )
327 tmp_t(i,j) = tmp_t(i,j) + dv*temp(i,j,k)
328 tmp_s(i,j) = tmp_s(i,j) + dv*salt(i,j,k)
329 endif
330 if (h_minimum > h(i,j,k)) h_minimum = h(i,j,k)
331 endif
332 enddo ; enddo ; enddo
333 area = reproducing_sum( tmp_a, unscale=us%L_to_m**2 )
334 vol = reproducing_sum( tmp_v, unscale=us%L_to_m**2*gv%H_to_mks )
335 if (do_ts) then
336 call min_across_pes( t%minimum ) ; call max_across_pes( t%maximum )
337 call min_across_pes( s%minimum ) ; call max_across_pes( s%maximum )
338 t%average = t_scale*reproducing_sum( tmp_t, unscale=us%C_to_degC*us%L_to_m**2*gv%H_to_mks) / vol
339 s%average = s_scale*reproducing_sum( tmp_s, unscale=us%S_to_ppt*us%L_to_m**2*gv%H_to_mks) / vol
340 endif
341 if (is_root_pe()) then
342 if (.not.firstcall) then
343 dv = vol - oldvol
344 delt%minimum = t%minimum - oldt%minimum ; delt%maximum = t%maximum - oldt%maximum
345 delt%average = t%average - oldt%average
346 dels%minimum = s%minimum - olds%minimum ; dels%maximum = s%maximum - olds%maximum
347 dels%average = s%average - olds%average
348 write(lmsg(1:80),'(2(a,es12.4))') 'Mean thickness =', gv%H_to_mks*vol/area,' frac. delta=',dv/vol
349 call mom_mesg(lmsg//trim(mesg))
350 if (do_ts) then
351 write(lmsg(1:80),'(a,3es12.4)') 'Temp min/mean/max =',t%minimum,t%average,t%maximum
352 call mom_mesg(lmsg//trim(mesg))
353 write(lmsg(1:80),'(a,3es12.4)') 'delT min/mean/max =',delt%minimum,delt%average,delt%maximum
354 call mom_mesg(lmsg//trim(mesg))
355 write(lmsg(1:80),'(a,3es12.4)') 'Salt min/mean/max =',s%minimum,s%average,s%maximum
356 call mom_mesg(lmsg//trim(mesg))
357 write(lmsg(1:80),'(a,3es12.4)') 'delS min/mean/max =',dels%minimum,dels%average,dels%maximum
358 call mom_mesg(lmsg//trim(mesg))
359 endif
360 else
361 write(lmsg(1:80),'(a,es12.4)') 'Mean thickness =', gv%H_to_mks*vol/area
362 call mom_mesg(lmsg//trim(mesg))
363 if (do_ts) then
364 write(lmsg(1:80),'(a,3es12.4)') 'Temp min/mean/max =', t%minimum, t%average, t%maximum
365 call mom_mesg(lmsg//trim(mesg))
366 write(lmsg(1:80),'(a,3es12.4)') 'Salt min/mean/max =', s%minimum, s%average, s%maximum
367 call mom_mesg(lmsg//trim(mesg))
368 endif
369 endif
370 endif
371 firstcall = .false. ; oldvol = vol
372 oldt%minimum = t%minimum ; oldt%maximum = t%maximum ; oldt%average = t%average
373 olds%minimum = s%minimum ; olds%maximum = s%maximum ; olds%average = s%average
374
375 if (do_ts .and. t%minimum<-5.0) then
376 do j=js,je ; do i=is,ie
377 if (minval(t_scale*temp(i,j,:)) == t%minimum) then
378 write(0,'(a,2f12.5)') 'x,y=', g%geoLonT(i,j), g%geoLatT(i,j)
379 write(0,'(a3,3a12)') 'k','h','Temp','Salt'
380 do k = 1, nz
381 write(0,'(I0," ",3es12.4)') k, h(i,j,k), t_scale*temp(i,j,k), s_scale*salt(i,j,k)
382 enddo
383 stop 'Extremum detected'
384 endif
385 enddo ; enddo
386 endif
387
388 if (h_minimum<0.0) then
389 do j=js,je ; do i=is,ie
390 if (minval(h(i,j,:)) == h_minimum) then
391 write(0,'(a,2f12.5)') 'x,y=',g%geoLonT(i,j),g%geoLatT(i,j)
392 write(0,'(a3,3a12)') 'k','h','Temp','Salt'
393 do k = 1, nz
394 write(0,'(I0," ",3es12.4)') k, h(i,j,k), t_scale*temp(i,j,k), s_scale*salt(i,j,k)
395 enddo
396 stop 'Negative thickness detected'
397 endif
398 enddo ; enddo
399 endif
400
401end subroutine mom_state_stats
402
403end module mom_checksum_packages