Line data Source code
1 : !=======================================================================
2 :
3 : ! Diagnostic information output during run
4 : !
5 : ! authors: Elizabeth C. Hunke, LANL
6 : ! Bruce P. Briegleb, NCAR
7 : !
8 : ! 2004: Block structure added by William Lipscomb
9 : ! 2006: Converted to free source form (F90) by Elizabeth Hunke
10 :
11 : module ice_diagnostics
12 :
13 : use ice_kinds_mod
14 : use ice_blocks, only: nx_block, ny_block
15 : use ice_communicate, only: my_task, master_task
16 : use ice_constants, only: c0, c1
17 : use ice_calendar, only: istep1
18 : use ice_domain_size, only: nslyr
19 : use ice_fileunits, only: nu_diag
20 : use ice_fileunits, only: flush_fileunit
21 : use ice_exit, only: abort_ice
22 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
23 : use icepack_intfc, only: icepack_max_aero, icepack_max_iso
24 : use icepack_intfc, only: icepack_query_parameters
25 : use icepack_intfc, only: icepack_query_tracer_flags
26 : use icepack_intfc, only: icepack_query_tracer_indices
27 :
28 : implicit none
29 : private
30 : public :: runtime_diags, init_mass_diags, init_diags, debug_ice, &
31 : print_state, diagnostic_abort
32 :
33 : ! diagnostic output file
34 : character (len=char_len), public :: diag_file
35 :
36 : ! point print data
37 :
38 : logical (kind=log_kind), public :: &
39 : debug_model , & ! if true, debug model at high level ! LCOV_EXCL_LINE
40 : print_points , & ! if true, print point data ! LCOV_EXCL_LINE
41 : print_global ! if true, print global data
42 :
43 : integer (kind=int_kind), public :: &
44 : debug_model_step = 0 ! begin printing at istep1=debug_model_step
45 :
46 : integer (kind=int_kind), parameter, public :: &
47 : npnt = 2 ! total number of points to be printed
48 :
49 : ! Set to true to identify unstable fast-moving ice.
50 : logical (kind=log_kind), parameter :: &
51 : check_umax = .false. ! if true, check for speed > umax_stab
52 :
53 : real (kind=dbl_kind), parameter :: &
54 : umax_stab = 1.0_dbl_kind , & ! ice speed threshold for instability (m/s) ! LCOV_EXCL_LINE
55 : aice_extmin = 0.15_dbl_kind ! min aice value for ice extent calc
56 :
57 : real (kind=dbl_kind), dimension(npnt), public :: &
58 : latpnt , & ! latitude of diagnostic points ! LCOV_EXCL_LINE
59 : lonpnt ! longitude of diagnostic points
60 :
61 : integer (kind=int_kind) :: &
62 : iindx , & ! i index for points ! LCOV_EXCL_LINE
63 : jindx , & ! j index for points ! LCOV_EXCL_LINE
64 : bindx ! block index for points
65 :
66 : ! for water and heat budgets
67 : real (kind=dbl_kind), dimension(npnt) :: &
68 : pdhi , & ! change in mean ice thickness (m) ! LCOV_EXCL_LINE
69 : pdhs , & ! change in mean snow thickness (m) ! LCOV_EXCL_LINE
70 : pde ! change in ice and snow energy (W m-2)
71 :
72 : real (kind=dbl_kind), dimension(npnt), public :: &
73 : plat, plon ! latitude, longitude of points
74 :
75 : integer (kind=int_kind), dimension(npnt), public :: &
76 : piloc, pjloc, pbloc, pmloc ! location of diagnostic points
77 :
78 : integer (kind=int_kind), public :: &
79 : debug_model_i = -1, & ! location of debug_model point, local i index ! LCOV_EXCL_LINE
80 : debug_model_j = -1, & ! location of debug_model point, local j index ! LCOV_EXCL_LINE
81 : debug_model_iblk = -1, & ! location of debug_model point, local block number ! LCOV_EXCL_LINE
82 : debug_model_task = -1 ! location of debug_model point, local task number
83 :
84 : ! for hemispheric water and heat budgets
85 : real (kind=dbl_kind) :: &
86 : totmn , & ! total ice/snow water mass (nh) ! LCOV_EXCL_LINE
87 : totms , & ! total ice/snow water mass (sh) ! LCOV_EXCL_LINE
88 : totmin , & ! total ice water mass (nh) ! LCOV_EXCL_LINE
89 : totmis , & ! total ice water mass (sh) ! LCOV_EXCL_LINE
90 : totsn , & ! total salt mass (nh) ! LCOV_EXCL_LINE
91 : totss , & ! total salt mass (sh) ! LCOV_EXCL_LINE
92 : toten , & ! total ice/snow energy (J) ! LCOV_EXCL_LINE
93 : totes ! total ice/snow energy (J)
94 :
95 : real (kind=dbl_kind), dimension(icepack_max_iso) :: &
96 : totison , & ! total isotope mass ! LCOV_EXCL_LINE
97 : totisos ! total isotope mass
98 :
99 : real (kind=dbl_kind), dimension(icepack_max_aero) :: &
100 : totaeron , & ! total aerosol mass ! LCOV_EXCL_LINE
101 : totaeros ! total aerosol mass
102 :
103 : !=======================================================================
104 :
105 : contains
106 :
107 : !=======================================================================
108 :
109 : ! Writes diagnostic info (max, min, global sums, etc) to standard out
110 : !
111 : ! authors: Elizabeth C. Hunke, LANL
112 : ! Bruce P. Briegleb, NCAR
113 : ! Cecilia M. Bitz, UW
114 :
115 5784 : subroutine runtime_diags (dt)
116 :
117 : use ice_arrays_column, only: floe_rad_c
118 : use ice_broadcast, only: broadcast_scalar
119 : use ice_constants, only: c1, c1000, c2, p001, p5, &
120 : field_loc_center, m2_to_km2
121 : use ice_domain, only: distrb_info, nblocks
122 : use ice_domain_size, only: ncat, n_iso, n_aero, max_blocks, nfsd
123 : use ice_flux, only: alvdr, alidr, alvdf, alidf, evap, fsnow, frazil, &
124 : fswabs, fswthru, flw, flwout, fsens, fsurf, flat, frzmlt_init, frain, fpond, & ! LCOV_EXCL_LINE
125 : fhocn_ai, fsalt_ai, fresh_ai, frazil_diag, & ! LCOV_EXCL_LINE
126 : update_ocn_f, cpl_frazil, Tair, Qa, fsw, fcondtop, meltt, meltb, meltl, snoice, & ! LCOV_EXCL_LINE
127 : dsnow, congel, sst, sss, Tf, fhocn, & ! LCOV_EXCL_LINE
128 : swvdr, swvdf, swidr, swidf, & ! LCOV_EXCL_LINE
129 : alvdr_init, alvdf_init, alidr_init, alidf_init
130 : use ice_flux_bgc, only: faero_atm, faero_ocn, fiso_atm, fiso_ocn
131 : use ice_global_reductions, only: global_sum, global_sum_prod, global_maxval
132 : use ice_grid, only: lmask_n, lmask_s, tarean, tareas, grid_ice, grid_average_X2Y
133 : use ice_state ! everything
134 : ! tcraig, this is likely to cause circular dependency because ice_prescribed_mod is high level routine
135 : #ifdef CESMCOUPLED
136 : use ice_prescribed_mod, only: prescribed_ice
137 : #endif
138 :
139 : real (kind=dbl_kind), intent(in) :: &
140 : dt ! time step
141 :
142 : ! local variables
143 :
144 : integer (kind=int_kind) :: &
145 : i, j, k, n, iblk, nc, & ! LCOV_EXCL_LINE
146 : ktherm, & ! LCOV_EXCL_LINE
147 : nt_tsfc, nt_aero, nt_fbri, nt_apnd, nt_hpnd, nt_fsd, & ! LCOV_EXCL_LINE
148 : nt_isosno, nt_isoice, nt_rsnw, nt_rhos, nt_smice, nt_smliq
149 :
150 : logical (kind=log_kind) :: &
151 : tr_pond_topo, tr_brine, tr_iso, tr_aero, calc_Tsfc, tr_fsd, & ! LCOV_EXCL_LINE
152 : tr_snow, snwgrain
153 :
154 : real (kind=dbl_kind) :: &
155 : rhow, rhos, rhoi, puny, awtvdr, awtidr, awtvdf, awtidf, & ! LCOV_EXCL_LINE
156 4320 : rhofresh, lfresh, lvap, ice_ref_salinity, Tffresh
157 :
158 : character (len=char_len) :: &
159 : snwredist, saltflux_option
160 :
161 : ! hemispheric state quantities
162 : real (kind=dbl_kind) :: &
163 : umaxn, hmaxn, shmaxn, arean, snwmxn, extentn, shmaxnt, & ! LCOV_EXCL_LINE
164 : umaxs, hmaxs, shmaxs, areas, snwmxs, extents, shmaxst, & ! LCOV_EXCL_LINE
165 : etotn, mtotn, micen, msnwn, pmaxn, ketotn, & ! LCOV_EXCL_LINE
166 : etots, mtots, mices, msnws, pmaxs, ketots, & ! LCOV_EXCL_LINE
167 : stotn, & ! LCOV_EXCL_LINE
168 : stots, & ! LCOV_EXCL_LINE
169 : urmsn, albtotn, arean_alb, mpndn, ptotn, spondn, & ! LCOV_EXCL_LINE
170 5760 : urmss, albtots, areas_alb, mpnds, ptots, sponds
171 :
172 : ! hemispheric flux quantities
173 : real (kind=dbl_kind) :: &
174 : rnn, snn, frzn, hnetn, fhocnn, fhatmn, fhfrzn, & ! LCOV_EXCL_LINE
175 : rns, sns, frzs, hnets, fhocns, fhatms, fhfrzs, & ! LCOV_EXCL_LINE
176 : fswnetn, fswnets, fswdnn, fswdns, swerrn, swerrs, & ! LCOV_EXCL_LINE
177 : sfsaltn, sfreshn, evpn, fluxn , delmxn, delmin, & ! LCOV_EXCL_LINE
178 : sfsalts, sfreshs, evps, fluxs , delmxs, delmis, & ! LCOV_EXCL_LINE
179 : delein, werrn, herrn, msltn, delmsltn, serrn, & ! LCOV_EXCL_LINE
180 7200 : deleis, werrs, herrs, mslts, delmslts, serrs
181 :
182 : ! isotope diagnostics
183 : real (kind=dbl_kind), dimension(icepack_max_aero) :: &
184 : fisoan, fisoon, isorn, & ! LCOV_EXCL_LINE
185 : fisoas, fisoos, isors, & ! LCOV_EXCL_LINE
186 : isomx1n, isomx1s, & ! LCOV_EXCL_LINE
187 18720 : isototn, isotots
188 :
189 : ! aerosol diagnostics
190 : real (kind=dbl_kind), dimension(icepack_max_aero) :: &
191 : faeran, faeron, aerrn, & ! LCOV_EXCL_LINE
192 : faeras, faeros, aerrs, & ! LCOV_EXCL_LINE
193 : aeromx1n, aeromx1s, & ! LCOV_EXCL_LINE
194 18720 : aerototn, aerotots
195 :
196 : ! fields at diagnostic points
197 : real (kind=dbl_kind), dimension(npnt) :: &
198 : paice, pTair, pQa, pfsnow, pfrain, pfsw, pflw, & ! LCOV_EXCL_LINE
199 : pTsfc, pevap, pfswabs, pflwout, pflat, pfsens, & ! LCOV_EXCL_LINE
200 : pfsurf, pfcondtop, psst, psss, pTf, hiavg, hsavg, hbravg, & ! LCOV_EXCL_LINE
201 : pfhocn, psalt, fsdavg, & ! LCOV_EXCL_LINE
202 : pmeltt, pmeltb, pmeltl, psnoice, pdsnow, pfrazil, pcongel, & ! LCOV_EXCL_LINE
203 17280 : prsnwavg, prhosavg, psmicetot, psmliqtot, psmtot
204 :
205 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
206 : uvelT, vvelT, & ! u,v on T points ! LCOV_EXCL_LINE
207 10021008 : work1, work2 ! temporary
208 :
209 : real (kind=dbl_kind), parameter :: &
210 : maxval_spval = -0.9_dbl_kind*HUGE(0.0_dbl_kind) ! spval to detect
211 : ! undefined values returned from global_maxval. if global_maxval
212 : ! is applied to a region that does not exist (for instance
213 : ! southern hemisphere in box cases), global_maxval
214 : ! returns -HUGE which we want to avoid writing. The
215 : ! return value is checked against maxval_spval before writing.
216 :
217 : character(len=*), parameter :: subname = '(runtime_diags)'
218 :
219 5784 : call icepack_query_parameters(ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc)
220 : call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_aero_out=tr_aero, &
221 : tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd, tr_iso_out=tr_iso, & ! LCOV_EXCL_LINE
222 5784 : tr_snow_out=tr_snow)
223 : call icepack_query_tracer_indices(nt_fbri_out=nt_fbri, nt_Tsfc_out=nt_Tsfc, &
224 : nt_aero_out=nt_aero, nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, & ! LCOV_EXCL_LINE
225 : nt_fsd_out=nt_fsd,nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, & ! LCOV_EXCL_LINE
226 : nt_rsnw_out=nt_rsnw, nt_rhos_out=nt_rhos, & ! LCOV_EXCL_LINE
227 5784 : nt_smice_out=nt_smice, nt_smliq_out=nt_smliq)
228 : call icepack_query_parameters(Tffresh_out=Tffresh, rhos_out=rhos, &
229 : rhow_out=rhow, rhoi_out=rhoi, puny_out=puny, & ! LCOV_EXCL_LINE
230 : awtvdr_out=awtvdr, awtidr_out=awtidr, awtvdf_out=awtvdf, awtidf_out=awtidf, & ! LCOV_EXCL_LINE
231 : rhofresh_out=rhofresh, lfresh_out=lfresh, lvap_out=lvap, & ! LCOV_EXCL_LINE
232 : ice_ref_salinity_out=ice_ref_salinity,snwredist_out=snwredist, & ! LCOV_EXCL_LINE
233 5784 : snwgrain_out=snwgrain, saltflux_option_out=saltflux_option)
234 5784 : call icepack_warnings_flush(nu_diag)
235 5784 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
236 0 : file=__FILE__, line=__LINE__)
237 :
238 : !-----------------------------------------------------------------
239 : ! state of the ice
240 : !-----------------------------------------------------------------
241 : ! hemispheric quantities
242 :
243 : ! total ice area
244 5784 : arean = c0
245 5784 : areas = c0
246 5784 : arean = global_sum(aice, distrb_info, field_loc_center, tarean)
247 5784 : areas = global_sum(aice, distrb_info, field_loc_center, tareas)
248 5784 : arean = arean * m2_to_km2
249 5784 : areas = areas * m2_to_km2
250 :
251 : ! ice extent (= area of grid cells with aice > aice_extmin)
252 16221984 : work1(:,:,:) = c0
253 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
254 8688 : do iblk = 1, nblocks
255 210240 : do j = 1, ny_block
256 7151880 : do i = 1, nx_block
257 7146096 : if (aice(i,j,iblk) >= aice_extmin) work1(i,j,iblk) = c1
258 : enddo
259 : enddo
260 : enddo
261 : !$OMP END PARALLEL DO
262 5784 : extentn = c0
263 5784 : extents = c0
264 5784 : extentn = global_sum(work1, distrb_info, field_loc_center, tarean)
265 5784 : extents = global_sum(work1, distrb_info, field_loc_center, tareas)
266 5784 : extentn = extentn * m2_to_km2
267 5784 : extents = extents * m2_to_km2
268 :
269 : ! total ice volume
270 5784 : shmaxn = c0
271 5784 : shmaxs = c0
272 5784 : shmaxn = global_sum(vice, distrb_info, field_loc_center, tarean)
273 5784 : shmaxs = global_sum(vice, distrb_info, field_loc_center, tareas)
274 :
275 : ! total snow volume
276 5784 : snwmxn = c0
277 5784 : snwmxs = c0
278 5784 : snwmxn = global_sum(vsno, distrb_info, field_loc_center, tarean)
279 5784 : snwmxs = global_sum(vsno, distrb_info, field_loc_center, tareas)
280 :
281 : ! total pond volume
282 5784 : ptotn = c0
283 5784 : ptots = c0
284 5784 : if (tr_pond_topo) then
285 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,n)
286 0 : do iblk = 1, nblocks
287 0 : do j = 1, ny_block
288 0 : do i = 1, nx_block
289 0 : work1(i,j,iblk) = c0
290 0 : do n = 1, ncat
291 : work1(i,j,iblk) = work1(i,j,iblk) &
292 : + aicen(i,j,n,iblk) & ! LCOV_EXCL_LINE
293 : * trcrn(i,j,nt_apnd,n,iblk) & ! LCOV_EXCL_LINE
294 0 : * trcrn(i,j,nt_hpnd,n,iblk)
295 : enddo
296 : enddo
297 : enddo
298 : enddo
299 : !$OMP END PARALLEL DO
300 0 : ptotn = global_sum(work1, distrb_info, field_loc_center, tarean)
301 0 : ptots = global_sum(work1, distrb_info, field_loc_center, tareas)
302 : endif
303 :
304 : ! total ice-snow kinetic energy, on T points.
305 5784 : if (grid_ice == 'B') then
306 5784 : call grid_average_X2Y('A',uvel ,'U',uvelT,'T')
307 5784 : call grid_average_X2Y('A',vvel ,'U',vvelT,'T')
308 0 : elseif (grid_ice == 'C') then
309 0 : call grid_average_X2Y('A',uvelE,'E',uvelT,'T')
310 0 : call grid_average_X2Y('A',vvelN,'N',vvelT,'T')
311 0 : elseif (grid_ice == 'CD') then
312 0 : call grid_average_X2Y('A',uvelE,'E',uvelN,'N',uvelT,'T')
313 0 : call grid_average_X2Y('A',vvelE,'E',vvelN,'N',vvelT,'T')
314 : endif
315 :
316 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
317 8688 : do iblk = 1, nblocks
318 210240 : do j = 1, ny_block
319 7151880 : do i = 1, nx_block
320 : work1(i,j,iblk) = p5 * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) &
321 7146096 : * (uvelT(i,j,iblk)**2 + vvelT(i,j,iblk)**2)
322 : enddo
323 : enddo
324 : enddo
325 : !$OMP END PARALLEL DO
326 5784 : ketotn = c0
327 5784 : ketots = c0
328 5784 : ketotn = global_sum(work1, distrb_info, field_loc_center, tarean)
329 5784 : ketots = global_sum(work1, distrb_info, field_loc_center, tareas)
330 :
331 : ! rms ice speed
332 5784 : urmsn = c2*ketotn/(rhoi*shmaxn + rhos*snwmxn + puny)
333 5784 : if (urmsn > puny) then
334 5784 : urmsn = sqrt(urmsn)
335 : else
336 0 : urmsn = c0
337 : endif
338 :
339 5784 : urmss = c2*ketots/(rhoi*shmaxs + rhos*snwmxs + puny)
340 5784 : if (urmss > puny) then
341 2904 : urmss = sqrt(urmss)
342 : else
343 2880 : urmss = c0
344 : endif
345 :
346 : ! average ice albedo
347 : ! mask out cells where sun is below horizon (for delta-Eddington)
348 :
349 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
350 8688 : do iblk = 1, nblocks
351 210240 : do j = 1, ny_block
352 7151880 : do i = 1, nx_block
353 : work1(i,j,iblk) = alvdr(i,j,iblk)*awtvdr &
354 : + alidr(i,j,iblk)*awtidr & ! LCOV_EXCL_LINE
355 : + alvdf(i,j,iblk)*awtvdf & ! LCOV_EXCL_LINE
356 6947424 : + alidf(i,j,iblk)*awtidf
357 7146096 : if (work1(i,j,iblk) > puny) then
358 5765985 : work2(i,j,iblk) = tarean(i,j,iblk)
359 : else
360 1181439 : work2(i,j,iblk) = c0
361 : endif
362 : enddo
363 : enddo
364 : enddo
365 : !$OMP END PARALLEL DO
366 :
367 5784 : arean_alb = global_sum(aice, distrb_info, field_loc_center, work2)
368 :
369 0 : albtotn = global_sum_prod(aice, work1, distrb_info, &
370 5784 : field_loc_center, work2)
371 :
372 5784 : if (arean_alb > c0) then
373 5784 : albtotn = albtotn / arean_alb
374 : else
375 0 : albtotn = c0
376 : endif
377 :
378 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
379 8688 : do iblk = 1, nblocks
380 210240 : do j = 1, ny_block
381 7151880 : do i = 1, nx_block
382 7146096 : if (work1(i,j,iblk) > puny) then
383 5765985 : work2(i,j,iblk) = tareas(i,j,iblk)
384 : else
385 1181439 : work2(i,j,iblk) = c0
386 : endif
387 : enddo
388 : enddo
389 : enddo
390 : !$OMP END PARALLEL DO
391 :
392 5784 : areas_alb = global_sum(aice, distrb_info, field_loc_center, work2)
393 :
394 0 : albtots = global_sum_prod(aice, work1, distrb_info, &
395 5784 : field_loc_center, work2)
396 :
397 5784 : if (areas_alb > c0) then
398 2904 : albtots = albtots / areas_alb
399 : else
400 2880 : albtots = c0
401 : endif
402 :
403 : ! maximum ice volume (= mean thickness including open water)
404 5784 : hmaxn = c0
405 5784 : hmaxs = c0
406 5784 : hmaxn = global_maxval(vice, distrb_info, lmask_n)
407 5784 : hmaxs = global_maxval(vice, distrb_info, lmask_s)
408 5784 : if (hmaxn < maxval_spval) hmaxn = c0
409 5784 : if (hmaxs < maxval_spval) hmaxs = c0
410 :
411 : ! maximum ice speed
412 5784 : if (grid_ice == 'CD') then
413 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
414 0 : do iblk = 1, nblocks
415 0 : do j = 1, ny_block
416 0 : do i = 1, nx_block
417 : work1(i,j,iblk) = max(sqrt(uvelE(i,j,iblk)**2 + vvelE(i,j,iblk)**2), &
418 0 : sqrt(uvelN(i,j,iblk)**2 + vvelN(i,j,iblk)**2))
419 : enddo
420 : enddo
421 : enddo
422 : !$OMP END PARALLEL DO
423 5784 : elseif (grid_ice == 'C') then
424 : ! map uvelE to N and vvelN to E then compute max on E and N
425 0 : call grid_average_X2Y('A',uvelE,'E',work1,'N') ! work1 =~ uvelN
426 0 : call grid_average_X2Y('A',vvelN,'N',work2,'E') ! work2 =~ vvelE
427 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
428 0 : do iblk = 1, nblocks
429 0 : do j = 1, ny_block
430 0 : do i = 1, nx_block
431 : work1(i,j,iblk) = max(sqrt(uvelE(i,j,iblk)**2 + work2(i,j,iblk)**2), &
432 0 : sqrt(work1(i,j,iblk)**2 + vvelN(i,j,iblk)**2))
433 : enddo
434 : enddo
435 : enddo
436 : !$OMP END PARALLEL DO
437 : else
438 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
439 8688 : do iblk = 1, nblocks
440 210240 : do j = 1, ny_block
441 7151880 : do i = 1, nx_block
442 7146096 : work1(i,j,iblk) = sqrt(uvel(i,j,iblk)**2 + vvel(i,j,iblk)**2)
443 : enddo
444 : enddo
445 : enddo
446 : !$OMP END PARALLEL DO
447 : endif
448 :
449 5784 : umaxn = c0
450 5784 : umaxs = c0
451 5784 : umaxn = global_maxval(work1, distrb_info, lmask_n)
452 5784 : umaxs = global_maxval(work1, distrb_info, lmask_s)
453 5784 : if (umaxn < maxval_spval) umaxn = c0
454 5784 : if (umaxs < maxval_spval) umaxs = c0
455 :
456 : ! Write warning message if ice speed is too big
457 : ! (Ice speeds of ~1 m/s or more usually indicate instability)
458 :
459 : if (check_umax) then
460 : if (umaxn > umax_stab) then
461 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
462 : do iblk = 1, nblocks
463 : do j = 1, ny_block
464 : do i = 1, nx_block
465 : if (abs(work1(i,j,iblk) - umaxn) < puny) then
466 : write(nu_diag,*) ' '
467 : write(nu_diag,*) 'Warning, large ice speed'
468 : write(nu_diag,*) 'my_task, iblk, i, j, umaxn:', &
469 : my_task, iblk, i, j, umaxn
470 : endif
471 : enddo
472 : enddo
473 : enddo
474 : !$OMP END PARALLEL DO
475 : elseif (umaxs > umax_stab) then
476 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
477 : do iblk = 1, nblocks
478 : do j = 1, ny_block
479 : do i = 1, nx_block
480 : if (abs(work1(i,j,iblk) - umaxs) < puny) then
481 : write(nu_diag,*) ' '
482 : write(nu_diag,*) 'Warning, large ice speed'
483 : write(nu_diag,*) 'my_task, iblk, i, j, umaxs:', &
484 : my_task, iblk, i, j, umaxs
485 : endif
486 : enddo
487 : enddo
488 : enddo
489 : !$OMP END PARALLEL DO
490 : endif ! umax
491 : endif ! check_umax
492 :
493 : ! maximum ice strength
494 :
495 5784 : pmaxn = c0
496 5784 : pmaxs = c0
497 5784 : pmaxn = global_maxval(strength, distrb_info, lmask_n)
498 5784 : pmaxs = global_maxval(strength, distrb_info, lmask_s)
499 5784 : if (pmaxn < maxval_spval) pmaxn = c0
500 5784 : if (pmaxs < maxval_spval) pmaxs = c0
501 :
502 5784 : pmaxn = pmaxn / c1000 ! convert to kN/m
503 5784 : pmaxs = pmaxs / c1000
504 :
505 5784 : if (print_global) then
506 :
507 : ! total ice/snow internal energy
508 5784 : call total_energy (work1)
509 :
510 5784 : etotn = c0
511 5784 : etots = c0
512 0 : etotn = global_sum(work1, distrb_info, &
513 5784 : field_loc_center, tarean)
514 0 : etots = global_sum(work1, distrb_info, &
515 5784 : field_loc_center, tareas)
516 :
517 : ! total salt volume
518 5784 : call total_salt (work2)
519 :
520 0 : stotn = global_sum(work2, distrb_info, &
521 5784 : field_loc_center, tarean)
522 0 : stots = global_sum(work2, distrb_info, &
523 5784 : field_loc_center, tareas)
524 :
525 :
526 : !-----------------------------------------------------------------
527 : ! various fluxes
528 : !-----------------------------------------------------------------
529 : ! evap, fsens, and flwout need to be multiplied by aice because
530 : ! regrettably they have been divided by aice for the coupler
531 : !-----------------------------------------------------------------
532 :
533 : ! evaporation
534 :
535 5784 : evpn = c0
536 5784 : evps = c0
537 : evpn = global_sum_prod(evap, aice, distrb_info, &
538 5784 : field_loc_center, tarean)
539 : evps = global_sum_prod(evap, aice, distrb_info, &
540 5784 : field_loc_center, tareas)
541 5784 : evpn = evpn*dt
542 5784 : evps = evps*dt
543 :
544 : ! total brine tracer
545 5784 : shmaxnt = c0
546 5784 : shmaxst = c0
547 5784 : if (tr_brine) then
548 0 : shmaxnt = global_sum(vice(:,:,:)*trcr(:,:,nt_fbri,:), distrb_info, &
549 0 : field_loc_center, tarean)
550 0 : shmaxst = global_sum(vice(:,:,:)*trcr(:,:,nt_fbri,:), distrb_info, &
551 0 : field_loc_center, tareas)
552 : endif
553 :
554 : ! salt flux
555 5784 : sfsaltn = c0
556 5784 : sfsalts = c0
557 : sfsaltn = global_sum(fsalt_ai, distrb_info, &
558 5784 : field_loc_center, tarean)
559 : sfsalts = global_sum(fsalt_ai, distrb_info, &
560 5784 : field_loc_center, tareas)
561 5784 : sfsaltn = sfsaltn*dt
562 5784 : sfsalts = sfsalts*dt
563 :
564 : ! fresh water flux
565 5784 : sfreshn = c0
566 5784 : sfreshs = c0
567 : sfreshn = global_sum(fresh_ai, distrb_info, &
568 5784 : field_loc_center, tarean)
569 : sfreshs = global_sum(fresh_ai, distrb_info, &
570 5784 : field_loc_center, tareas)
571 5784 : sfreshn = sfreshn*dt
572 5784 : sfreshs = sfreshs*dt
573 :
574 : ! pond water flux
575 5784 : spondn = c0
576 5784 : sponds = c0
577 5784 : if (tr_pond_topo) then
578 : spondn = global_sum(fpond, distrb_info, &
579 0 : field_loc_center, tarean)
580 : sponds = global_sum(fpond, distrb_info, &
581 0 : field_loc_center, tareas)
582 0 : spondn = spondn*dt
583 0 : sponds = sponds*dt
584 : endif
585 :
586 : ! ocean heat
587 : ! Note: fswthru not included because it does not heat ice
588 5784 : fhocnn = c0
589 5784 : fhocns = c0
590 : fhocnn = global_sum(fhocn_ai, distrb_info, &
591 5784 : field_loc_center, tarean)
592 : fhocns = global_sum(fhocn_ai, distrb_info, &
593 5784 : field_loc_center, tareas)
594 :
595 : ! latent heat
596 : ! You may be wondering, where is the latent heat flux?
597 : ! It is not included here because it cancels with
598 : ! the evaporative flux times the enthalpy of the
599 : ! ice/snow that evaporated.
600 :
601 : ! atmo heat flux
602 : ! Note: flwout includes the reflected longwave down, needed by the
603 : ! atmosphere as an upwards radiative boundary condition.
604 : ! Also note: fswabs includes solar radiation absorbed in ocean,
605 : ! which must be subtracted here.
606 :
607 5784 : if (calc_Tsfc) then
608 :
609 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
610 8688 : do iblk = 1, nblocks
611 210240 : do j = 1, ny_block
612 7151880 : do i = 1, nx_block
613 : work1(i,j,iblk) = &
614 : (fswabs(i,j,iblk) - fswthru(i,j,iblk) & ! LCOV_EXCL_LINE
615 : + fsens (i,j,iblk) + flwout (i,j,iblk)) & ! LCOV_EXCL_LINE
616 : * aice (i,j,iblk) & ! LCOV_EXCL_LINE
617 7146096 : + flw (i,j,iblk) * aice_init (i,j,iblk)
618 : enddo
619 : enddo
620 : enddo
621 : !$OMP END PARALLEL DO
622 :
623 : else ! fsurf is computed by atmosphere model
624 :
625 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
626 0 : do iblk = 1, nblocks
627 0 : do j = 1, ny_block
628 0 : do i = 1, nx_block
629 : work1(i,j,iblk) = &
630 : (fsurf(i,j,iblk) - flat(i,j,iblk)) & ! LCOV_EXCL_LINE
631 0 : * aice(i,j,iblk)
632 : enddo
633 : enddo
634 : enddo
635 : !$OMP END PARALLEL DO
636 :
637 : endif ! calc_Tsfc
638 :
639 5784 : fhatmn = c0
640 5784 : fhatms = c0
641 0 : fhatmn = global_sum(work1, distrb_info, &
642 5784 : field_loc_center, tarean)
643 0 : fhatms = global_sum(work1, distrb_info, &
644 5784 : field_loc_center, tareas)
645 :
646 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
647 8688 : do iblk = 1, nblocks
648 210240 : do j = 1, ny_block
649 7151880 : do i = 1, nx_block
650 : work1(i,j,iblk) = &
651 7146096 : fswabs(i,j,iblk) * aice(i,j,iblk)
652 : enddo
653 : enddo
654 : enddo
655 : !$OMP END PARALLEL DO
656 :
657 5784 : fswnetn = c0
658 5784 : fswnets = c0
659 0 : fswnetn = global_sum(work1, distrb_info, &
660 5784 : field_loc_center, tarean)
661 0 : fswnets = global_sum(work1, distrb_info, &
662 5784 : field_loc_center, tareas)
663 :
664 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
665 8688 : do iblk = 1, nblocks
666 210240 : do j = 1, ny_block
667 7151880 : do i = 1, nx_block
668 : work1(i,j,iblk) = (aice_init(i,j,iblk)-alvdr_init(i,j,iblk))*swvdr(i,j,iblk) &
669 : + (aice_init(i,j,iblk)-alidr_init(i,j,iblk))*swidr(i,j,iblk) & ! LCOV_EXCL_LINE
670 : + (aice_init(i,j,iblk)-alvdf_init(i,j,iblk))*swvdf(i,j,iblk) & ! LCOV_EXCL_LINE
671 7146096 : + (aice_init(i,j,iblk)-alidf_init(i,j,iblk))*swidf(i,j,iblk)
672 : enddo
673 : enddo
674 : enddo
675 : !$OMP END PARALLEL DO
676 :
677 5784 : fswdnn = c0
678 5784 : fswdns = c0
679 0 : fswdnn = global_sum(work1, distrb_info, &
680 5784 : field_loc_center, tarean)
681 0 : fswdns = global_sum(work1, distrb_info, &
682 5784 : field_loc_center, tareas)
683 :
684 : ! freezing potential
685 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
686 8688 : do iblk = 1, nblocks
687 210240 : do j = 1, ny_block
688 7151880 : do i = 1, nx_block
689 7146096 : work1(i,j,iblk) = max(c0,frzmlt_init(i,j,iblk))
690 : enddo
691 : enddo
692 : enddo
693 : !$OMP END PARALLEL DO
694 :
695 5784 : fhfrzn = c0
696 5784 : fhfrzs = c0
697 0 : fhfrzn = global_sum(work1, distrb_info, &
698 5784 : field_loc_center, tarean)
699 0 : fhfrzs = global_sum(work1, distrb_info, &
700 5784 : field_loc_center, tareas)
701 :
702 : ! rain
703 5784 : rnn = c0
704 5784 : rns = c0
705 : rnn = global_sum_prod(frain, aice_init, distrb_info, &
706 5784 : field_loc_center, tarean)
707 : rns = global_sum_prod(frain, aice_init, distrb_info, &
708 5784 : field_loc_center, tareas)
709 5784 : rnn = rnn*dt
710 5784 : rns = rns*dt
711 :
712 : ! snow
713 5784 : snn = c0
714 5784 : sns = c0
715 : snn = global_sum_prod(fsnow, aice_init, distrb_info, &
716 5784 : field_loc_center, tarean)
717 : sns = global_sum_prod(fsnow, aice_init, distrb_info, &
718 5784 : field_loc_center, tareas)
719 5784 : snn = snn*dt
720 5784 : sns = sns*dt
721 :
722 : ! frazil ice growth !! should not be multiplied by aice
723 : ! m/step->kg/m^2/s
724 16221984 : work1(:,:,:) = frazil(:,:,:)*rhoi/dt
725 5784 : if (.not. update_ocn_f .and. ktherm == 2 .and. cpl_frazil == 'fresh_ice_correction') then
726 16221984 : work1(:,:,:) = (frazil(:,:,:)-frazil_diag(:,:,:))*rhoi/dt
727 : endif
728 5784 : frzn = c0
729 5784 : frzs = c0
730 0 : frzn = global_sum(work1, distrb_info, &
731 5784 : field_loc_center, tarean)
732 0 : frzs = global_sum(work1, distrb_info, &
733 5784 : field_loc_center, tareas)
734 5784 : frzn = frzn*dt
735 5784 : frzs = frzs*dt
736 :
737 : ! ice, snow, pond mass
738 5784 : micen = rhoi*shmaxn
739 5784 : msnwn = rhos*snwmxn
740 5784 : mices = rhoi*shmaxs
741 5784 : msnws = rhos*snwmxs
742 5784 : mpndn = rhofresh*ptotn
743 5784 : mpnds = rhofresh*ptots
744 :
745 : ! total ice, snow and pond mass
746 5784 : mtotn = micen + msnwn + mpndn
747 5784 : mtots = mices + msnws + mpnds
748 :
749 : ! mass change since beginning of time step
750 5784 : delmin = mtotn - totmn
751 5784 : delmis = mtots - totms
752 :
753 : ! ice mass change including frazil ice formation
754 5784 : delmxn = micen - totmin
755 5784 : delmxs = mices - totmis
756 5784 : if (.not. update_ocn_f) then
757 : ! ice mass change excluding frazil ice formation
758 5784 : delmxn = delmxn - frzn
759 5784 : delmxs = delmxs - frzs
760 : endif
761 :
762 : ! total water flux
763 5784 : fluxn = c0
764 5784 : fluxs = c0
765 5784 : if( arean > c0) then
766 : ! water associated with frazil ice included in fresh
767 5784 : fluxn = rnn + snn + evpn - sfreshn
768 5784 : if (.not. update_ocn_f) then
769 5784 : fluxn = fluxn + frzn
770 : endif
771 : endif
772 5784 : if( areas > c0) then
773 : ! water associated with frazil ice included in fresh
774 2904 : fluxs = rns + sns + evps - sfreshs
775 2904 : if (.not. update_ocn_f) then
776 2904 : fluxs = fluxs + frzs
777 : endif
778 : endif
779 :
780 5784 : werrn = (fluxn-delmin)/(mtotn + c1)
781 5784 : werrs = (fluxs-delmis)/(mtots + c1)
782 :
783 : ! energy change
784 5784 : delein = etotn - toten
785 5784 : deleis = etots - totes
786 :
787 5784 : fhatmn = fhatmn + ( - snn * Lfresh + evpn * Lvap ) / dt
788 5784 : fhatms = fhatms + ( - sns * Lfresh + evps * Lvap ) / dt
789 :
790 5784 : hnetn = (fhatmn - fhocnn - fhfrzn) * dt
791 5784 : hnets = (fhatms - fhocns - fhfrzs) * dt
792 :
793 5784 : herrn = (hnetn - delein) / (etotn - c1)
794 5784 : herrs = (hnets - deleis) / (etots - c1)
795 :
796 5784 : swerrn = (fswnetn - fswdnn) / (fswnetn - c1)
797 5784 : swerrs = (fswnets - fswdns) / (fswnets - c1)
798 :
799 : ! salt mass
800 5784 : if (saltflux_option == 'prognostic') then
801 : ! compute the total salt mass
802 0 : msltn = stotn*rhoi*p001
803 0 : mslts = stots*rhoi*p001
804 :
805 : ! change in salt mass
806 0 : delmsltn = rhoi*(stotn-totsn)*p001
807 0 : delmslts = rhoi*(stots-totss)*p001
808 : else
809 5784 : msltn = micen*ice_ref_salinity*p001
810 5784 : mslts = mices*ice_ref_salinity*p001
811 :
812 : ! change in salt mass
813 5784 : delmsltn = delmxn*ice_ref_salinity*p001
814 5784 : delmslts = delmxs*ice_ref_salinity*p001
815 : endif
816 :
817 : ! salt error
818 5784 : serrn = (sfsaltn + delmsltn) / (msltn + c1)
819 5784 : serrs = (sfsalts + delmslts) / (mslts + c1)
820 :
821 : ! isotopes
822 5784 : if (tr_iso) then
823 0 : fisoan = c0
824 0 : fisoas = c0
825 0 : fisoon = c0
826 0 : fisoos = c0
827 0 : isototn = c0
828 0 : isotots = c0
829 0 : isomx1n = c0
830 0 : isomx1s = c0
831 0 : isorn = c0
832 0 : isors = c0
833 0 : do n = 1, n_iso
834 0 : fisoan(n) = global_sum_prod(fiso_atm(:,:,n,:), aice_init, &
835 0 : distrb_info, field_loc_center, tarean)
836 0 : fisoas(n) = global_sum_prod(fiso_atm(:,:,n,:), aice_init, &
837 0 : distrb_info, field_loc_center, tareas)
838 0 : fisoan(n) = fisoan(n)*dt
839 0 : fisoas(n) = fisoas(n)*dt
840 0 : fisoon(n) = global_sum_prod(fiso_ocn(:,:,n,:), aice, &
841 0 : distrb_info, field_loc_center, tarean)
842 0 : fisoos(n) = global_sum_prod(fiso_ocn(:,:,n,:), aice, &
843 0 : distrb_info, field_loc_center, tareas)
844 0 : fisoon(n) = fisoon(n)*dt
845 0 : fisoos(n) = fisoos(n)*dt
846 :
847 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,k)
848 0 : do iblk = 1, nblocks
849 0 : do j = 1, ny_block
850 0 : do i = 1, nx_block
851 0 : work1(i,j,iblk) = c0
852 0 : do k = 1, n_iso
853 : work1(i,j,iblk) = work1(i,j,iblk) &
854 : + vsno(i,j,iblk)*trcr(i,j,nt_isosno+k-1,iblk) & ! LCOV_EXCL_LINE
855 0 : + vice(i,j,iblk)*trcr(i,j,nt_isoice+k-1,iblk)
856 : enddo
857 : enddo
858 : enddo
859 : enddo
860 : !$OMP END PARALLEL DO
861 0 : isototn(n) = global_sum(work1, distrb_info, field_loc_center, tarean)
862 0 : isotots(n) = global_sum(work1, distrb_info, field_loc_center, tareas)
863 0 : isomx1n(n) = global_maxval(work1, distrb_info, lmask_n)
864 0 : isomx1s(n) = global_maxval(work1, distrb_info, lmask_s)
865 0 : if (isomx1n(n) < maxval_spval) isomx1n(n) = c0
866 0 : if (isomx1s(n) < maxval_spval) isomx1s(n) = c0
867 0 : isorn(n) = (totison(n)-isototn(n)+fisoan(n)-fisoon(n))/(isototn(n)+c1)
868 0 : isors(n) = (totisos(n)-isotots(n)+fisoas(n)-fisoos(n))/(isotots(n)+c1)
869 : enddo ! n_iso
870 : endif ! tr_iso
871 :
872 : ! aerosols
873 5784 : if (tr_aero) then
874 0 : faeran = c0
875 0 : faeras = c0
876 0 : faeron = c0
877 0 : faeros = c0
878 0 : aerototn = c0
879 0 : aerotots = c0
880 0 : aeromx1n = c0
881 0 : aeromx1s = c0
882 0 : aerrn = c0
883 0 : aerrs = c0
884 0 : do n = 1, n_aero
885 0 : faeran(n) = global_sum_prod(faero_atm(:,:,n,:), aice_init, &
886 0 : distrb_info, field_loc_center, tarean)
887 0 : faeras(n) = global_sum_prod(faero_atm(:,:,n,:), aice_init, &
888 0 : distrb_info, field_loc_center, tareas)
889 0 : faeran(n) = faeran(n)*dt
890 0 : faeras(n) = faeras(n)*dt
891 0 : faeron(n) = global_sum_prod(faero_ocn(:,:,n,:), aice, &
892 0 : distrb_info, field_loc_center, tarean)
893 0 : faeros(n) = global_sum_prod(faero_ocn(:,:,n,:), aice, &
894 0 : distrb_info, field_loc_center, tareas)
895 0 : faeron(n) = faeron(n)*dt
896 0 : faeros(n) = faeros(n)*dt
897 :
898 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
899 0 : do iblk = 1, nblocks
900 0 : do j = 1, ny_block
901 0 : do i = 1, nx_block
902 : work1(i,j,iblk) = &
903 : trcr(i,j,nt_aero +4*(n-1),iblk)*vsno(i,j,iblk) & ! LCOV_EXCL_LINE
904 : + trcr(i,j,nt_aero+1+4*(n-1),iblk)*vsno(i,j,iblk) & ! LCOV_EXCL_LINE
905 : + trcr(i,j,nt_aero+2+4*(n-1),iblk)*vice(i,j,iblk) & ! LCOV_EXCL_LINE
906 0 : + trcr(i,j,nt_aero+3+4*(n-1),iblk)*vice(i,j,iblk)
907 : enddo
908 : enddo
909 : enddo
910 : !$OMP END PARALLEL DO
911 0 : aerototn(n) = global_sum(work1, distrb_info, field_loc_center, tarean)
912 0 : aerotots(n) = global_sum(work1, distrb_info, field_loc_center, tareas)
913 0 : aeromx1n(n) = global_maxval(work1, distrb_info, lmask_n)
914 0 : aeromx1s(n) = global_maxval(work1, distrb_info, lmask_s)
915 0 : if (aeromx1n(n) < maxval_spval) aeromx1n(n) = c0
916 0 : if (aeromx1s(n) < maxval_spval) aeromx1s(n) = c0
917 :
918 0 : aerrn(n) = (totaeron(n)-aerototn(n)+faeran(n)-faeron(n)) &
919 0 : / (aerototn(n) + c1)
920 0 : aerrs(n) = (totaeros(n)-aerotots(n)+faeras(n)-faeros(n)) &
921 0 : / (aerotots(n) + c1)
922 : enddo ! n_aero
923 : endif ! tr_aero
924 :
925 : endif ! print_global
926 :
927 5784 : if (print_points) then
928 :
929 : !-----------------------------------------------------------------
930 : ! state of the ice and associated fluxes for 2 defined points
931 : ! NOTE these are computed for the last timestep only (not avg)
932 : !-----------------------------------------------------------------
933 :
934 5784 : call total_energy (work1)
935 5784 : call total_salt (work2)
936 :
937 17352 : do n = 1, npnt
938 11568 : if (my_task == pmloc(n)) then
939 1968 : i = piloc(n)
940 1968 : j = pjloc(n)
941 1968 : iblk = pbloc(n)
942 :
943 1968 : pTair(n) = Tair(i,j,iblk) - Tffresh ! air temperature
944 1968 : pQa(n) = Qa(i,j,iblk) ! specific humidity
945 1968 : pfsnow(n) = fsnow(i,j,iblk)*dt/rhos ! snowfall
946 1968 : pfrain(n) = frain(i,j,iblk)*dt/rhow ! rainfall
947 1968 : pfsw(n) = fsw(i,j,iblk) ! shortwave radiation
948 1968 : pflw(n) = flw(i,j,iblk) ! longwave radiation
949 1968 : paice(n) = aice(i,j,iblk) ! ice area
950 :
951 1968 : fsdavg(n) = c0 ! avg floe effective radius
952 1968 : hiavg(n) = c0 ! avg snow/ice thickness
953 1968 : hsavg(n) = c0
954 1968 : hbravg(n) = c0 ! avg brine thickness
955 1968 : if (paice(n) /= c0) then
956 1968 : hiavg(n) = vice(i,j,iblk)/paice(n)
957 1968 : hsavg(n) = vsno(i,j,iblk)/paice(n)
958 1968 : if (tr_brine) hbravg(n) = trcr(i,j,nt_fbri,iblk)* hiavg(n)
959 1968 : if (tr_fsd) then
960 : ! not sure why this does not work
961 : ! do k = 1, nfsd
962 : ! fsdavg(n) = fsdavg(n) &
963 : ! + trcr(i,j,nt_fsd+k-1,iblk) * floe_rad_c(k) & ! LCOV_EXCL_LINE
964 : ! / trcr(i,j,nt_fsd+k-1,iblk)
965 : ! enddo
966 : ! this works
967 0 : do nc = 1, ncat
968 0 : do k = 1, nfsd
969 0 : fsdavg(n) = fsdavg(n) &
970 : + trcrn(i,j,nt_fsd+k-1,nc,iblk) * floe_rad_c(k) & ! LCOV_EXCL_LINE
971 0 : * aicen(i,j,nc,iblk) / paice(n)
972 : enddo
973 : enddo
974 : endif
975 : endif
976 1968 : if (tr_snow) then ! snow tracer quantities
977 0 : prsnwavg (n) = c0 ! avg snow grain radius
978 0 : prhosavg (n) = c0 ! avg snow density
979 0 : psmicetot(n) = c0 ! total mass of ice in snow (kg/m2)
980 0 : psmliqtot(n) = c0 ! total mass of liquid in snow (kg/m2)
981 0 : psmtot (n) = c0 ! total mass of snow volume (kg/m2)
982 0 : if (vsno(i,j,iblk) > c0) then
983 0 : do k = 1, nslyr
984 0 : prsnwavg (n) = prsnwavg (n) + trcr(i,j,nt_rsnw +k-1,iblk) ! snow grain radius
985 0 : prhosavg (n) = prhosavg (n) + trcr(i,j,nt_rhos +k-1,iblk) ! compacted snow density
986 0 : psmicetot(n) = psmicetot(n) + trcr(i,j,nt_smice+k-1,iblk) * vsno(i,j,iblk)
987 0 : psmliqtot(n) = psmliqtot(n) + trcr(i,j,nt_smliq+k-1,iblk) * vsno(i,j,iblk)
988 : end do
989 : endif
990 0 : psmtot (n) = rhos * vsno(i,j,iblk) ! mass of ice in standard density snow
991 0 : prsnwavg (n) = prsnwavg (n) / real(nslyr,kind=dbl_kind) ! snow grain radius
992 0 : prhosavg (n) = prhosavg (n) / real(nslyr,kind=dbl_kind) ! compacted snow density
993 0 : psmicetot(n) = psmicetot(n) / real(nslyr,kind=dbl_kind) ! mass of ice in snow
994 0 : psmliqtot(n) = psmliqtot(n) / real(nslyr,kind=dbl_kind) ! mass of liquid in snow
995 : end if
996 1968 : psalt(n) = c0
997 1968 : if (vice(i,j,iblk) /= c0) psalt(n) = work2(i,j,iblk)/vice(i,j,iblk)
998 1968 : pTsfc(n) = trcr(i,j,nt_Tsfc,iblk) ! ice/snow sfc temperature
999 1968 : pevap(n) = evap(i,j,iblk)*dt/rhoi ! sublimation/condensation
1000 1968 : pfswabs(n) = fswabs(i,j,iblk) ! absorbed solar flux
1001 1968 : pflwout(n) = flwout(i,j,iblk) ! outward longwave flux
1002 1968 : pflat(n) = flat(i,j,iblk) ! latent heat flux
1003 1968 : pfsens(n) = fsens(i,j,iblk) ! sensible heat flux
1004 1968 : pfsurf(n) = fsurf(i,j,iblk) ! total sfc heat flux
1005 1968 : pfcondtop(n) = fcondtop(i,j,iblk) ! top sfc cond flux
1006 1968 : pmeltt(n) = meltt(i,j,iblk) ! top melt
1007 1968 : pmeltb(n) = meltb(i,j,iblk) ! bottom melt
1008 1968 : pmeltl(n) = meltl(i,j,iblk) ! lateral melt
1009 1968 : psnoice(n) = snoice(i,j,iblk) ! snow ice
1010 1968 : pdsnow(n) = dsnow(i,j,iblk) ! snow change
1011 1968 : pfrazil(n) = frazil(i,j,iblk) ! frazil ice
1012 1968 : pcongel(n) = congel(i,j,iblk) ! congelation ice
1013 1968 : pdhi(n) = vice(i,j,iblk) - pdhi(n) ! ice thickness change
1014 1968 : pdhs(n) = vsno(i,j,iblk) - pdhs(n) ! snow thickness change
1015 1968 : pde(n) =-(work1(i,j,iblk)- pde(n))/dt ! ice/snow energy change
1016 1968 : psst(n) = sst(i,j,iblk) ! sea surface temperature
1017 1968 : psss(n) = sss(i,j,iblk) ! sea surface salinity
1018 1968 : pTf(n) = Tf(i,j,iblk) ! freezing temperature
1019 1968 : pfhocn(n) = -fhocn(i,j,iblk) ! ocean heat used by ice
1020 :
1021 : endif ! my_task = pmloc
1022 :
1023 11568 : call broadcast_scalar(pTair (n), pmloc(n))
1024 11568 : call broadcast_scalar(pQa (n), pmloc(n))
1025 11568 : call broadcast_scalar(pfsnow (n), pmloc(n))
1026 11568 : call broadcast_scalar(pfrain (n), pmloc(n))
1027 11568 : call broadcast_scalar(pfsw (n), pmloc(n))
1028 11568 : call broadcast_scalar(pflw (n), pmloc(n))
1029 11568 : call broadcast_scalar(paice (n), pmloc(n))
1030 11568 : call broadcast_scalar(hsavg (n), pmloc(n))
1031 11568 : call broadcast_scalar(hiavg (n), pmloc(n))
1032 11568 : call broadcast_scalar(fsdavg (n), pmloc(n))
1033 11568 : call broadcast_scalar(psalt (n), pmloc(n))
1034 11568 : call broadcast_scalar(hbravg (n), pmloc(n))
1035 11568 : call broadcast_scalar(pTsfc (n), pmloc(n))
1036 11568 : call broadcast_scalar(pevap (n), pmloc(n))
1037 11568 : call broadcast_scalar(pfswabs (n), pmloc(n))
1038 11568 : call broadcast_scalar(pflwout (n), pmloc(n))
1039 11568 : call broadcast_scalar(pflat (n), pmloc(n))
1040 11568 : call broadcast_scalar(pfsens (n), pmloc(n))
1041 11568 : call broadcast_scalar(pfsurf (n), pmloc(n))
1042 11568 : call broadcast_scalar(pfcondtop(n), pmloc(n))
1043 11568 : call broadcast_scalar(pmeltt (n), pmloc(n))
1044 11568 : call broadcast_scalar(pmeltb (n), pmloc(n))
1045 11568 : call broadcast_scalar(pmeltl (n), pmloc(n))
1046 11568 : call broadcast_scalar(psnoice (n), pmloc(n))
1047 11568 : call broadcast_scalar(pdsnow (n), pmloc(n))
1048 11568 : call broadcast_scalar(psmtot (n), pmloc(n))
1049 11568 : call broadcast_scalar(prsnwavg (n), pmloc(n))
1050 11568 : call broadcast_scalar(prhosavg (n), pmloc(n))
1051 11568 : call broadcast_scalar(psmicetot(n), pmloc(n))
1052 11568 : call broadcast_scalar(psmliqtot(n), pmloc(n))
1053 11568 : call broadcast_scalar(pfrazil (n), pmloc(n))
1054 11568 : call broadcast_scalar(pcongel (n), pmloc(n))
1055 11568 : call broadcast_scalar(pdhi (n), pmloc(n))
1056 11568 : call broadcast_scalar(pdhs (n), pmloc(n))
1057 11568 : call broadcast_scalar(pde (n), pmloc(n))
1058 11568 : call broadcast_scalar(psst (n), pmloc(n))
1059 11568 : call broadcast_scalar(psss (n), pmloc(n))
1060 11568 : call broadcast_scalar(pTf (n), pmloc(n))
1061 123912 : call broadcast_scalar(pfhocn (n), pmloc(n))
1062 :
1063 : enddo ! npnt
1064 : endif ! print_points
1065 :
1066 : !-----------------------------------------------------------------
1067 : ! start spewing
1068 : !-----------------------------------------------------------------
1069 :
1070 5784 : if (my_task == master_task) then
1071 :
1072 984 : write(nu_diag,899) 'Arctic','Antarctic'
1073 :
1074 984 : write(nu_diag,901) 'total ice area (km^2) = ',arean, areas
1075 984 : write(nu_diag,901) 'total ice extent(km^2) = ',extentn,extents
1076 984 : write(nu_diag,901) 'total ice volume (m^3) = ',shmaxn, shmaxs
1077 984 : write(nu_diag,901) 'total snw volume (m^3) = ',snwmxn, snwmxs
1078 984 : write(nu_diag,901) 'tot kinetic energy (J) = ',ketotn, ketots
1079 984 : write(nu_diag,900) 'rms ice speed (m/s) = ',urmsn, urmss
1080 984 : write(nu_diag,900) 'average albedo = ',albtotn,albtots
1081 984 : write(nu_diag,900) 'max ice volume (m) = ',hmaxn, hmaxs
1082 984 : write(nu_diag,900) 'max ice speed (m/s) = ',umaxn, umaxs
1083 984 : write(nu_diag,900) 'max strength (kN/m) = ',pmaxn, pmaxs
1084 :
1085 984 : if (print_global) then ! global diags for conservations checks
1086 :
1087 : ! tcraig, this is likely to cause circular dependency because ice_prescribed_mod is high level routine
1088 : #ifdef CESMCOUPLED
1089 : if (prescribed_ice) then
1090 : write (nu_diag,*) '----------------------------'
1091 : write (nu_diag,*) 'This is the prescribed ice option.'
1092 : write (nu_diag,*) 'Heat and water will not be conserved.'
1093 : write (nu_diag,*) '----------------------------'
1094 : endif
1095 : #endif
1096 :
1097 984 : write(nu_diag,*) '----------------------------'
1098 984 : write(nu_diag,901) 'arwt rain h2o kg in dt = ',rnn,rns
1099 984 : write(nu_diag,901) 'arwt snow h2o kg in dt = ',snn,sns
1100 984 : write(nu_diag,901) 'arwt evap h2o kg in dt = ',evpn,evps
1101 984 : write(nu_diag,901) 'arwt frzl h2o kg in dt = ',frzn,frzs
1102 984 : if (tr_pond_topo) &
1103 0 : write(nu_diag,901) 'arwt fpnd h2o kg in dt = ',spondn,sponds
1104 984 : write(nu_diag,901) 'arwt frsh h2o kg in dt = ',sfreshn,sfreshs
1105 :
1106 984 : write(nu_diag,901) 'arwt ice mass (kg) = ',micen,mices
1107 984 : write(nu_diag,901) 'arwt snw mass (kg) = ',msnwn,msnws
1108 984 : if (tr_pond_topo) &
1109 0 : write(nu_diag,901) 'arwt pnd mass (kg) = ',mpndn,mpnds
1110 :
1111 984 : write(nu_diag,901) 'arwt tot mass (kg) = ',mtotn,mtots
1112 984 : write(nu_diag,901) 'arwt tot mass chng(kg) = ',delmin,delmis
1113 984 : write(nu_diag,901) 'arwt water flux = ',fluxn,fluxs
1114 984 : if (update_ocn_f) then
1115 0 : write (nu_diag,*) '(=rain+snow+evap-fresh) '
1116 : else
1117 984 : write (nu_diag,*) '(=rain+snow+evap+frzl-fresh) '
1118 : endif
1119 984 : write(nu_diag,901) 'water flux error = ',werrn,werrs
1120 :
1121 984 : write(nu_diag,*) '----------------------------'
1122 984 : write(nu_diag,901) 'arwt atm heat flux (W) = ',fhatmn,fhatms
1123 984 : write(nu_diag,901) 'arwt ocn heat flux (W) = ',fhocnn,fhocns
1124 984 : write(nu_diag,901) 'arwt frzl heat flux(W) = ',fhfrzn,fhfrzs
1125 984 : write(nu_diag,901) 'arwt tot energy (J) = ',etotn,etots
1126 984 : write(nu_diag,901) 'arwt net heat (J) = ',hnetn,hnets
1127 984 : write(nu_diag,901) 'arwt tot energy chng(J)= ',delein,deleis
1128 984 : write(nu_diag,901) 'arwt heat error = ',herrn,herrs
1129 984 : write(nu_diag,*) '----------------------------'
1130 984 : write(nu_diag,901) 'arwt incoming sw (W) = ',fswdnn,fswdns
1131 984 : write(nu_diag,901) 'arwt absorbed sw (W) = ',fswnetn,fswnets
1132 984 : write(nu_diag,901) 'arwt swdn error = ',swerrn,swerrs
1133 :
1134 984 : write(nu_diag,*) '----------------------------'
1135 984 : write(nu_diag,901) 'total brine tr (m^3) = ',shmaxnt, shmaxst
1136 984 : write(nu_diag,901) 'arwt salt mass (kg) = ',msltn,mslts
1137 984 : write(nu_diag,901) 'arwt salt mass chng(kg)= ',delmsltn, &
1138 1968 : delmslts
1139 984 : write(nu_diag,901) 'arwt salt flx in dt(kg)= ',sfsaltn, &
1140 1968 : sfsalts
1141 984 : write(nu_diag,901) 'arwt salt flx error = ',serrn,serrs
1142 :
1143 984 : write(nu_diag,*) '----------------------------'
1144 984 : if (tr_iso) then
1145 0 : do n = 1, n_iso
1146 0 : write(nu_diag,*) ' isotope ',n
1147 0 : write(nu_diag,901) 'fiso_atm (kg/m2) = ', fisoan(n), fisoas(n)
1148 0 : write(nu_diag,901) 'fiso_ocn (kg/m2) = ', fisoon(n), fisoos(n)
1149 0 : write(nu_diag,901) 'total iso (kg/m2) = ', isototn(n), isotots(n)
1150 0 : write(nu_diag,901) 'iso error = ', isorn(n), isors(n)
1151 0 : write(nu_diag,901) 'maximum iso (kg/m2) = ', isomx1n(n),isomx1s(n)
1152 : enddo
1153 0 : write(nu_diag,*) '----------------------------'
1154 : endif ! tr_iso
1155 984 : if (tr_aero) then
1156 0 : do n = 1, n_aero
1157 0 : write(nu_diag,*) ' aerosol ',n
1158 0 : write(nu_diag,901) 'faero_atm (kg/m2) = ', faeran(n), faeras(n)
1159 0 : write(nu_diag,901) 'faero_ocn (kg/m2) = ', faeron(n), faeros(n)
1160 0 : write(nu_diag,901) 'total aero (kg/m2) = ', aerototn(n), aerotots(n)
1161 0 : write(nu_diag,901) 'aero error = ', aerrn(n), aerrs(n)
1162 0 : write(nu_diag,901) 'maximum aero (kg/m2) = ', aeromx1n(n),aeromx1s(n)
1163 : enddo
1164 0 : write(nu_diag,*) '----------------------------'
1165 : endif ! tr_aero
1166 :
1167 : endif ! print_global
1168 :
1169 984 : call flush_fileunit(nu_diag)
1170 :
1171 : !-----------------------------------------------------------------
1172 : ! diagnostics for Arctic and Antarctic points
1173 : !-----------------------------------------------------------------
1174 :
1175 984 : if (print_points) then
1176 :
1177 984 : write(nu_diag,*) ' '
1178 984 : write(nu_diag,902) ' Lat, Long ',plat(1),plon(1), &
1179 1968 : plat(2),plon(2)
1180 984 : write(nu_diag,903) ' my_task, iblk, i, j ', &
1181 : pmloc(1),pbloc(1),piloc(1),pjloc(1), & ! LCOV_EXCL_LINE
1182 1968 : pmloc(2),pbloc(2),piloc(2),pjloc(2)
1183 984 : write(nu_diag,*) '----------atm----------'
1184 984 : write(nu_diag,900) 'air temperature (C) = ',pTair(1),pTair(2)
1185 984 : write(nu_diag,900) 'specific humidity = ',pQa(1),pQa(2)
1186 984 : write(nu_diag,900) 'snowfall (m) = ',pfsnow(1), &
1187 1968 : pfsnow(2)
1188 984 : write(nu_diag,900) 'rainfall (m) = ',pfrain(1), &
1189 1968 : pfrain(2)
1190 984 : if (.not.calc_Tsfc) then
1191 0 : write(nu_diag,900) 'total surface heat flux= ',pfsurf(1),pfsurf(2)
1192 0 : write(nu_diag,900) 'top sfc conductive flux= ',pfcondtop(1), &
1193 0 : pfcondtop(2)
1194 0 : write(nu_diag,900) 'latent heat flx = ',pflat(1),pflat(2)
1195 : else
1196 984 : write(nu_diag,900) 'shortwave radiation sum= ',pfsw(1),pfsw(2)
1197 984 : write(nu_diag,900) 'longwave radiation = ',pflw(1),pflw(2)
1198 : endif
1199 984 : write(nu_diag,*) '----------ice----------'
1200 984 : write(nu_diag,900) 'area fraction = ',paice(1),paice(2)
1201 984 : write(nu_diag,900) 'avg ice thickness (m) = ',hiavg(1),hiavg(2)
1202 984 : write(nu_diag,900) 'avg snow depth (m) = ',hsavg(1),hsavg(2)
1203 984 : write(nu_diag,900) 'avg salinity (ppt) = ',psalt(1),psalt(2)
1204 984 : write(nu_diag,900) 'avg brine thickness (m)= ',hbravg(1),hbravg(2)
1205 984 : if (tr_fsd) &
1206 0 : write(nu_diag,900) 'avg fsd rep radius (m) = ',fsdavg(1),fsdavg(2)
1207 :
1208 984 : if (calc_Tsfc) then
1209 984 : write(nu_diag,900) 'surface temperature(C) = ',pTsfc(1),pTsfc(2)
1210 984 : write(nu_diag,900) 'absorbed shortwave flx = ',pfswabs(1), &
1211 1968 : pfswabs(2)
1212 984 : write(nu_diag,900) 'outward longwave flx = ',pflwout(1), &
1213 1968 : pflwout(2)
1214 984 : write(nu_diag,900) 'sensible heat flx = ',pfsens(1), &
1215 1968 : pfsens(2)
1216 984 : write(nu_diag,900) 'latent heat flx = ',pflat(1),pflat(2)
1217 : endif
1218 984 : write(nu_diag,900) 'subl/cond (m ice) = ',pevap(1),pevap(2)
1219 984 : write(nu_diag,900) 'top melt (m) = ',pmeltt(1) &
1220 1968 : ,pmeltt(2)
1221 984 : write(nu_diag,900) 'bottom melt (m) = ',pmeltb(1) &
1222 1968 : ,pmeltb(2)
1223 984 : write(nu_diag,900) 'lateral melt (m) = ',pmeltl(1) &
1224 1968 : ,pmeltl(2)
1225 984 : write(nu_diag,900) 'new ice (m) = ',pfrazil(1), &
1226 1968 : pfrazil(2)
1227 984 : write(nu_diag,900) 'congelation (m) = ',pcongel(1), &
1228 1968 : pcongel(2)
1229 984 : write(nu_diag,900) 'snow-ice (m) = ',psnoice(1), &
1230 1968 : psnoice(2)
1231 984 : write(nu_diag,900) 'snow change (m) = ',pdsnow(1), &
1232 1968 : pdsnow(2)
1233 984 : write(nu_diag,900) 'effective dhi (m) = ',pdhi(1),pdhi(2)
1234 984 : write(nu_diag,900) 'effective dhs (m) = ',pdhs(1),pdhs(2)
1235 984 : write(nu_diag,900) 'intnl enrgy chng(W/m^2)= ',pde (1),pde (2)
1236 :
1237 984 : if (tr_snow) then
1238 0 : if (trim(snwredist) /= 'none') then
1239 0 : write(nu_diag,900) 'avg snow density(kg/m3)= ',prhosavg(1) &
1240 0 : ,prhosavg(2)
1241 : endif
1242 0 : if (snwgrain) then
1243 0 : write(nu_diag,900) 'avg snow grain radius = ',prsnwavg(1) &
1244 0 : ,prsnwavg(2)
1245 0 : write(nu_diag,900) 'mass ice in snow(kg/m2)= ',psmicetot(1) &
1246 0 : ,psmicetot(2)
1247 0 : write(nu_diag,900) 'mass liq in snow(kg/m2)= ',psmliqtot(1) &
1248 0 : ,psmliqtot(2)
1249 0 : write(nu_diag,900) 'mass std snow (kg/m2)= ',psmtot(1) &
1250 0 : ,psmtot(2)
1251 0 : write(nu_diag,900) 'max ice+liq (kg/m2)= ',rhow * hsavg(1) &
1252 0 : ,rhow * hsavg(2)
1253 : endif
1254 : endif
1255 :
1256 984 : write(nu_diag,*) '----------ocn----------'
1257 984 : write(nu_diag,900) 'sst (C) = ',psst(1),psst(2)
1258 984 : write(nu_diag,900) 'sss (ppt) = ',psss(1),psss(2)
1259 984 : write(nu_diag,900) 'freezing temp (C) = ',pTf(1),pTf(2)
1260 984 : write(nu_diag,900) 'heat used (W/m^2) = ',pfhocn(1), &
1261 1968 : pfhocn(2)
1262 :
1263 : endif ! print_points
1264 : endif ! my_task = master_task
1265 :
1266 : 899 format (27x,a24,2x,a24)
1267 : 900 format (a25,2x,f24.17,2x,f24.17)
1268 : 901 format (a25,2x,1pe24.17,2x,1pe24.17)
1269 : 902 format (a25,10x,f6.1,1x,f6.1,9x,f6.1,1x,f6.1)
1270 : 903 format (a25,5x,i4,1x,i4,1x,i4,1x,i4,7x,i4,1x,i4,1x,i4,1x,i4)
1271 :
1272 5784 : end subroutine runtime_diags
1273 :
1274 : !=======================================================================
1275 :
1276 : ! Computes global combined ice and snow mass sum
1277 : !
1278 : ! author: Elizabeth C. Hunke, LANL
1279 :
1280 5784 : subroutine init_mass_diags
1281 :
1282 : use ice_constants, only: field_loc_center
1283 : use ice_domain, only: distrb_info, nblocks
1284 : use ice_domain_size, only: n_iso, n_aero, ncat, max_blocks
1285 : use ice_global_reductions, only: global_sum
1286 : use ice_grid, only: tareas, tarean
1287 : use ice_state, only: aicen, vice, vsno, trcrn, trcr
1288 :
1289 : integer (kind=int_kind) :: n, i, j, k, iblk, &
1290 : nt_hpnd, nt_apnd, nt_aero, nt_isosno, nt_isoice
1291 :
1292 : logical (kind=log_kind) :: &
1293 : tr_iso, tr_aero, tr_pond_topo
1294 :
1295 : real (kind=dbl_kind) :: &
1296 : shmaxn, snwmxn, shmaxs, snwmxs, totpn, totps, & ! LCOV_EXCL_LINE
1297 1440 : rhoi, rhos, rhofresh
1298 :
1299 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
1300 10021008 : work1, work2
1301 :
1302 : character(len=*), parameter :: subname = '(init_mass_diags)'
1303 :
1304 5784 : call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_pond_topo_out=tr_pond_topo)
1305 5784 : call icepack_query_tracer_flags(tr_iso_out=tr_iso)
1306 : call icepack_query_tracer_indices(nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
1307 5784 : nt_hpnd_out=nt_hpnd, nt_apnd_out=nt_apnd, nt_aero_out=nt_aero)
1308 : call icepack_query_parameters( &
1309 5784 : rhoi_out=rhoi, rhos_out=rhos, rhofresh_out=rhofresh)
1310 5784 : call icepack_warnings_flush(nu_diag)
1311 5784 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1312 0 : file=__FILE__, line=__LINE__)
1313 :
1314 : ! total ice volume
1315 5784 : shmaxn = global_sum(vice, distrb_info, field_loc_center, tarean)
1316 5784 : shmaxs = global_sum(vice, distrb_info, field_loc_center, tareas)
1317 :
1318 : ! total snow volume
1319 5784 : snwmxn = global_sum(vsno, distrb_info, field_loc_center, tarean)
1320 5784 : snwmxs = global_sum(vsno, distrb_info, field_loc_center, tareas)
1321 :
1322 : ! north/south ice mass
1323 5784 : totmin = rhoi*shmaxn
1324 5784 : totmis = rhoi*shmaxs
1325 :
1326 : ! north/south ice+snow mass
1327 5784 : totmn = totmin + rhos*snwmxn
1328 5784 : totms = totmis + rhos*snwmxs
1329 :
1330 : ! north/south ice+snow energy
1331 5784 : call total_energy (work1)
1332 5784 : toten = global_sum(work1, distrb_info, field_loc_center, tarean)
1333 5784 : totes = global_sum(work1, distrb_info, field_loc_center, tareas)
1334 :
1335 : ! north/south salt
1336 5784 : call total_salt (work2)
1337 5784 : totsn = global_sum(work2, distrb_info, field_loc_center, tarean)
1338 5784 : totss = global_sum(work2, distrb_info, field_loc_center, tareas)
1339 :
1340 :
1341 5784 : if (print_points) then
1342 17352 : do n = 1, npnt
1343 17352 : if (my_task == pmloc(n)) then
1344 1968 : i = piloc(n)
1345 1968 : j = pjloc(n)
1346 1968 : iblk = pbloc(n)
1347 :
1348 1968 : pdhi(n) = vice(i,j,iblk)
1349 1968 : pdhs(n) = vsno(i,j,iblk)
1350 1968 : pde(n) = work1(i,j,iblk)
1351 : endif
1352 : enddo ! npnt
1353 : endif ! print_points
1354 :
1355 5784 : if (tr_iso) then
1356 0 : do n=1,n_iso
1357 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,k)
1358 0 : do iblk = 1, nblocks
1359 0 : do j = 1, ny_block
1360 0 : do i = 1, nx_block
1361 0 : work1(i,j,iblk) = c0
1362 0 : do k = 1, n_iso
1363 : work1(i,j,iblk) = work1(i,j,iblk) &
1364 : + vsno(i,j,iblk)*trcr(i,j,nt_isosno+k-1,iblk) & ! LCOV_EXCL_LINE
1365 0 : + vice(i,j,iblk)*trcr(i,j,nt_isoice+k-1,iblk)
1366 : enddo
1367 : enddo
1368 : enddo
1369 : enddo
1370 : !$OMP END PARALLEL DO
1371 0 : totison(n)= global_sum(work1, distrb_info, field_loc_center, tarean)
1372 0 : totisos(n)= global_sum(work1, distrb_info, field_loc_center, tareas)
1373 : enddo
1374 : endif
1375 :
1376 5784 : if (tr_aero) then
1377 0 : do n=1,n_aero
1378 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
1379 0 : do iblk = 1, nblocks
1380 0 : do j = 1, ny_block
1381 0 : do i = 1, nx_block
1382 : work1(i,j,iblk) = trcr(i,j,nt_aero +4*(n-1),iblk)*vsno(i,j,iblk) &
1383 : + trcr(i,j,nt_aero+1+4*(n-1),iblk)*vsno(i,j,iblk) & ! LCOV_EXCL_LINE
1384 : + trcr(i,j,nt_aero+2+4*(n-1),iblk)*vice(i,j,iblk) & ! LCOV_EXCL_LINE
1385 0 : + trcr(i,j,nt_aero+3+4*(n-1),iblk)*vice(i,j,iblk)
1386 : enddo
1387 : enddo
1388 : enddo
1389 : !$OMP END PARALLEL DO
1390 0 : totaeron(n)= global_sum(work1, distrb_info, field_loc_center, tarean)
1391 0 : totaeros(n)= global_sum(work1, distrb_info, field_loc_center, tareas)
1392 : enddo
1393 : endif
1394 :
1395 5784 : if (tr_pond_topo) then
1396 0 : totpn = c0
1397 0 : totps = c0
1398 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,n)
1399 0 : do iblk = 1, nblocks
1400 0 : do j = 1, ny_block
1401 0 : do i = 1, nx_block
1402 0 : work1(i,j,iblk) = c0
1403 0 : do n = 1, ncat
1404 : work1(i,j,iblk) = work1(i,j,iblk) &
1405 : + aicen(i,j,n,iblk) & ! LCOV_EXCL_LINE
1406 : * trcrn(i,j,nt_apnd,n,iblk) & ! LCOV_EXCL_LINE
1407 0 : * trcrn(i,j,nt_hpnd,n,iblk)
1408 : enddo
1409 : enddo
1410 : enddo
1411 : enddo
1412 : !$OMP END PARALLEL DO
1413 0 : totpn = global_sum(work1, distrb_info, field_loc_center, tarean)
1414 0 : totps = global_sum(work1, distrb_info, field_loc_center, tareas)
1415 :
1416 : ! north/south ice+snow+pond mass
1417 0 : totmn = totmn + totpn*rhofresh
1418 0 : totms = totms + totps*rhofresh
1419 : endif
1420 :
1421 5784 : end subroutine init_mass_diags
1422 :
1423 : !=======================================================================
1424 :
1425 : ! Computes total energy of ice and snow in a grid cell.
1426 : !
1427 : ! authors: E. C. Hunke, LANL
1428 :
1429 17352 : subroutine total_energy (work)
1430 :
1431 : use ice_domain, only: nblocks
1432 : use ice_domain_size, only: ncat, nilyr, nslyr, max_blocks
1433 : use ice_grid, only: tmask
1434 : use ice_state, only: vicen, vsnon, trcrn
1435 :
1436 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(out) :: &
1437 : work ! total energy
1438 :
1439 : ! local variables
1440 :
1441 : integer (kind=int_kind) :: &
1442 : icells ! number of ocean/ice cells
1443 :
1444 : integer (kind=int_kind), dimension (nx_block*ny_block) :: &
1445 : indxi, & ! compressed indices in i/j directions ! LCOV_EXCL_LINE
1446 34704 : indxj
1447 :
1448 : integer (kind=int_kind) :: &
1449 : i, j, k, n, iblk, ij, & ! LCOV_EXCL_LINE
1450 : nt_qice, nt_qsno
1451 :
1452 : character(len=*), parameter :: subname = '(total_energy)'
1453 :
1454 17352 : call icepack_query_tracer_indices(nt_qice_out=nt_qice, nt_qsno_out=nt_qsno)
1455 17352 : call icepack_warnings_flush(nu_diag)
1456 17352 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1457 0 : file=__FILE__, line=__LINE__)
1458 :
1459 8640 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,n,k,ij,icells,indxi,indxj)
1460 26064 : do iblk = 1, nblocks
1461 :
1462 : !-----------------------------------------------------------------
1463 : ! Initialize
1464 : !-----------------------------------------------------------------
1465 :
1466 17352 : icells = 0
1467 613368 : do j = 1, ny_block
1468 21455640 : do i = 1, nx_block
1469 21438288 : if (tmask(i,j,iblk)) then
1470 19594800 : icells = icells + 1
1471 19594800 : indxi(icells) = i
1472 19594800 : indxj(icells) = j
1473 : endif ! tmask
1474 : enddo
1475 : enddo
1476 :
1477 21455640 : work(:,:,iblk) = c0
1478 :
1479 : !-----------------------------------------------------------------
1480 : ! Aggregate
1481 : !-----------------------------------------------------------------
1482 :
1483 121464 : do n = 1, ncat
1484 694080 : do k = 1, nilyr
1485 686512080 : do ij = 1, icells
1486 685818000 : i = indxi(ij)
1487 685818000 : j = indxj(ij)
1488 : work(i,j,iblk) = work(i,j,iblk) &
1489 : + trcrn(i,j,nt_qice+k-1,n,iblk) & ! LCOV_EXCL_LINE
1490 686425320 : * vicen(i,j,n,iblk) / real(nilyr,kind=dbl_kind)
1491 : enddo ! ij
1492 : enddo ! k
1493 :
1494 190872 : do k = 1, nslyr
1495 98147520 : do ij = 1, icells
1496 97974000 : i = indxi(ij)
1497 97974000 : j = indxj(ij)
1498 : work(i,j,iblk) = work(i,j,iblk) &
1499 : + trcrn(i,j,nt_qsno+k-1,n,iblk) & ! LCOV_EXCL_LINE
1500 98060760 : * vsnon(i,j,n,iblk) / real(nslyr,kind=dbl_kind)
1501 : enddo ! ij
1502 : enddo ! k
1503 : enddo ! n
1504 :
1505 : enddo ! iblk
1506 : !$OMP END PARALLEL DO
1507 :
1508 17352 : end subroutine total_energy
1509 :
1510 : !=======================================================================
1511 :
1512 : ! Computes bulk salinity of ice and snow in a grid cell.
1513 : ! author: E. C. Hunke, LANL
1514 :
1515 17352 : subroutine total_salt (work)
1516 :
1517 : use ice_domain, only: nblocks
1518 : use ice_domain_size, only: ncat, nilyr, max_blocks
1519 : use ice_grid, only: tmask
1520 : use ice_state, only: vicen, trcrn
1521 :
1522 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(out) :: &
1523 : work ! total salt
1524 :
1525 : ! local variables
1526 :
1527 : integer (kind=int_kind) :: &
1528 : icells ! number of ocean/ice cells
1529 :
1530 : integer (kind=int_kind), dimension (nx_block*ny_block) :: &
1531 : indxi, & ! compressed indices in i/j directions ! LCOV_EXCL_LINE
1532 34704 : indxj
1533 :
1534 : integer (kind=int_kind) :: &
1535 : i, j, k, n, iblk, ij, & ! LCOV_EXCL_LINE
1536 : nt_sice
1537 :
1538 : character(len=*), parameter :: subname = '(total_salt)'
1539 :
1540 17352 : call icepack_query_tracer_indices(nt_sice_out=nt_sice)
1541 17352 : call icepack_warnings_flush(nu_diag)
1542 17352 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1543 0 : file=__FILE__, line=__LINE__)
1544 :
1545 8640 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,n,k,ij,icells,indxi,indxj)
1546 26064 : do iblk = 1, nblocks
1547 :
1548 : !-----------------------------------------------------------------
1549 : ! Initialize
1550 : !-----------------------------------------------------------------
1551 :
1552 17352 : icells = 0
1553 613368 : do j = 1, ny_block
1554 21455640 : do i = 1, nx_block
1555 21438288 : if (tmask(i,j,iblk)) then
1556 19594800 : icells = icells + 1
1557 19594800 : indxi(icells) = i
1558 19594800 : indxj(icells) = j
1559 : endif ! tmask
1560 : enddo
1561 : enddo
1562 :
1563 21455640 : work(:,:,iblk) = c0
1564 :
1565 : !-----------------------------------------------------------------
1566 : ! Aggregate
1567 : !-----------------------------------------------------------------
1568 :
1569 121464 : do n = 1, ncat
1570 711432 : do k = 1, nilyr
1571 686512080 : do ij = 1, icells
1572 685818000 : i = indxi(ij)
1573 685818000 : j = indxj(ij)
1574 : work(i,j,iblk) = work(i,j,iblk) &
1575 : + trcrn(i,j,nt_sice+k-1,n,iblk) & ! LCOV_EXCL_LINE
1576 686425320 : * vicen(i,j,n,iblk) / real(nilyr,kind=dbl_kind)
1577 : enddo ! ij
1578 : enddo ! k
1579 : enddo ! n
1580 :
1581 : enddo ! iblk
1582 : !$OMP END PARALLEL DO
1583 :
1584 17352 : end subroutine total_salt
1585 :
1586 : !=======================================================================
1587 :
1588 : ! Find tasks for diagnostic points.
1589 : !
1590 : ! authors: Elizabeth C. Hunke and William H. Lipscomb, LANL
1591 :
1592 37 : subroutine init_diags
1593 :
1594 : use ice_grid, only: hm, TLAT, TLON
1595 : use ice_blocks, only: block, get_block
1596 : use ice_constants, only: c180, c360, p5
1597 : use ice_domain, only: blocks_ice, distrb_info, nblocks
1598 : use ice_global_reductions, only: global_minval, global_maxval
1599 :
1600 : real (kind=dbl_kind) :: &
1601 : rad_to_deg, & ! LCOV_EXCL_LINE
1602 : puny, & ! LCOV_EXCL_LINE
1603 : latdis , & ! latitude distance ! LCOV_EXCL_LINE
1604 : londis , & ! longitude distance ! LCOV_EXCL_LINE
1605 : totdis , & ! total distance ! LCOV_EXCL_LINE
1606 : mindis , & ! minimum distance from desired location ! LCOV_EXCL_LINE
1607 8 : mindis_g ! global minimum distance from desired location
1608 :
1609 : integer (kind=int_kind) :: &
1610 : n , & ! index for point search ! LCOV_EXCL_LINE
1611 : i,j , & ! grid indices ! LCOV_EXCL_LINE
1612 : iblk , & ! block index ! LCOV_EXCL_LINE
1613 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
1614 :
1615 : type (block) :: &
1616 : this_block ! block information for current block
1617 :
1618 : character(len=*), parameter :: subname = '(init_diags)'
1619 :
1620 : !tcraig, do this all the time now for print_points_state usage
1621 : ! if (print_points) then
1622 :
1623 37 : call icepack_query_parameters(puny_out=puny, rad_to_deg_out=rad_to_deg)
1624 37 : call icepack_warnings_flush(nu_diag)
1625 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1626 0 : file=__FILE__, line=__LINE__)
1627 :
1628 37 : if (my_task==master_task) then
1629 7 : write(nu_diag,*) ' '
1630 7 : write(nu_diag,*) ' Find indices of diagnostic points '
1631 : endif
1632 :
1633 111 : piloc(:) = -1
1634 111 : pjloc(:) = -1
1635 111 : pbloc(:) = -1
1636 111 : pmloc(:) = -999
1637 111 : plat(:) = -999._dbl_kind
1638 111 : plon(:) = -999._dbl_kind
1639 :
1640 : ! find minimum distance to diagnostic points on this processor
1641 111 : do n = 1, npnt
1642 74 : if (lonpnt(n) > c180) lonpnt(n) = lonpnt(n) - c360
1643 :
1644 74 : iindx = 0
1645 74 : jindx = 0
1646 74 : bindx = 0
1647 74 : mindis = 540.0_dbl_kind ! 360. + 180.
1648 :
1649 74 : if (abs(latpnt(n)) < c360 .and. abs(lonpnt(n)) < c360) then
1650 :
1651 : ! MDT, 09/2017: Comment out OpenMP directives since loop is not thread-safe
1652 : ! This is computing closest point, Could add a CRITICAL but it's just initialization
1653 : !!$XXXOMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,latdis,londis,totdis)
1654 394 : do iblk = 1, nblocks
1655 320 : this_block = get_block(blocks_ice(iblk),iblk)
1656 320 : ilo = this_block%ilo
1657 320 : ihi = this_block%ihi
1658 320 : jlo = this_block%jlo
1659 320 : jhi = this_block%jhi
1660 :
1661 10040 : do j = jlo, jhi
1662 191212 : do i = ilo, ihi
1663 190892 : if (hm(i,j,iblk) > p5) then
1664 143660 : latdis = abs(latpnt(n)-TLAT(i,j,iblk)*rad_to_deg)
1665 32024 : londis = abs(lonpnt(n)-TLON(i,j,iblk)*rad_to_deg) &
1666 143660 : * cos(TLAT(i,j,iblk))
1667 143660 : totdis = sqrt(latdis**2 + londis**2)
1668 143660 : if (totdis < mindis) then
1669 13957 : mindis = totdis
1670 13957 : jindx = j
1671 13957 : iindx = i
1672 13957 : bindx = iblk
1673 : endif ! totdis < mindis
1674 : endif ! hm > p5
1675 : enddo ! i
1676 : enddo ! j
1677 : enddo ! iblk
1678 : !!$XXXOMP END PARALLEL DO
1679 :
1680 : endif
1681 :
1682 : ! find global minimum distance to diagnostic points
1683 74 : mindis_g = global_minval(mindis, distrb_info)
1684 :
1685 : ! save indices of minimum-distance grid cell
1686 74 : if (mindis <= 180.0 .and. abs(mindis_g - mindis) < puny) then
1687 14 : piloc(n) = iindx
1688 14 : pjloc(n) = jindx
1689 14 : pbloc(n) = bindx
1690 14 : pmloc(n) = my_task
1691 14 : plat(n) = TLAT(iindx,jindx,bindx)*rad_to_deg
1692 14 : plon(n) = TLON(iindx,jindx,bindx)*rad_to_deg
1693 : endif
1694 :
1695 : ! communicate to all processors
1696 74 : piloc(n) = global_maxval(piloc(n), distrb_info)
1697 74 : pjloc(n) = global_maxval(pjloc(n), distrb_info)
1698 74 : pbloc(n) = global_maxval(pbloc(n), distrb_info)
1699 74 : pmloc(n) = global_maxval(pmloc(n), distrb_info)
1700 74 : plat(n) = global_maxval(plat(n), distrb_info)
1701 74 : plon(n) = global_maxval(plon(n), distrb_info)
1702 :
1703 : ! write to log file
1704 111 : if (my_task==master_task) then
1705 14 : write(nu_diag,*) ' '
1706 14 : write(nu_diag,100) n,latpnt(n),lonpnt(n),plat(n),plon(n), &
1707 28 : piloc(n), pjloc(n), pbloc(n), pmloc(n)
1708 : endif
1709 : 100 format(' found point',i4/ &
1710 : ' lat lon TLAT TLON i j block task'/ & ! LCOV_EXCL_LINE
1711 : 4(f6.1,1x),1x,4(i4,2x) )
1712 :
1713 : enddo ! npnt
1714 : ! endif ! print_points
1715 :
1716 37 : end subroutine init_diags
1717 :
1718 : !=======================================================================
1719 :
1720 : ! This routine is useful for debugging
1721 : ! author Elizabeth C. Hunke, LANL
1722 :
1723 0 : subroutine debug_ice(iblk, plabeld)
1724 :
1725 : character (char_len), intent(in) :: plabeld
1726 : integer (kind=int_kind), intent(in) :: iblk
1727 :
1728 : ! local
1729 : character(len=*), parameter :: subname='(debug_ice)'
1730 :
1731 0 : if (istep1 >= debug_model_step) then
1732 :
1733 : ! set debug point to 1st global point if not set as local values
1734 : if (debug_model_i < 0 .and. debug_model_j < 0 .and. &
1735 0 : debug_model_iblk < 0 .and. debug_model_task < 0) then
1736 0 : debug_model_i = piloc(1)
1737 0 : debug_model_j = pjloc(1)
1738 0 : debug_model_task = pmloc(1)
1739 0 : debug_model_iblk = pbloc(1)
1740 : endif
1741 :
1742 : ! if debug point is messed up, abort
1743 : if (debug_model_i < 0 .or. debug_model_j < 0 .or. &
1744 0 : debug_model_iblk < 0 .or. debug_model_task < 0) then
1745 0 : call abort_ice (subname//'ERROR: debug_model_[i,j,iblk,mytask] not set correctly')
1746 : endif
1747 :
1748 : ! write out debug info
1749 0 : if (debug_model_iblk == iblk .and. debug_model_task == my_task) then
1750 0 : call print_state(plabeld,debug_model_i,debug_model_j,debug_model_iblk)
1751 : endif
1752 :
1753 : endif
1754 :
1755 0 : end subroutine debug_ice
1756 :
1757 : !=======================================================================
1758 :
1759 : ! This routine is useful for debugging.
1760 : ! author: Elizabeth C. Hunke, LANL
1761 :
1762 0 : subroutine print_state(plabel,i,j,iblk)
1763 :
1764 : use ice_grid, only: grid_ice
1765 : use ice_blocks, only: block, get_block
1766 : use ice_domain, only: blocks_ice
1767 : use ice_domain_size, only: ncat, nilyr, nslyr, nfsd
1768 : use ice_grid, only: TLAT, TLON
1769 : use ice_state, only: aice, aice0, aicen, vicen, vsnon, uvel, vvel, &
1770 : uvelE, vvelE, uvelN, vvelN, trcrn
1771 : use ice_flux, only: uatm, vatm, potT, Tair, Qa, flw, frain, fsnow, &
1772 : fsens, flat, evap, flwout, swvdr, swvdf, swidr, swidf, rhoa, & ! LCOV_EXCL_LINE
1773 : frzmlt, sst, sss, Tf, Tref, Qref, Uref, uocn, vocn, strtltxU, strtltyU
1774 :
1775 : character (len=20), intent(in) :: plabel
1776 :
1777 : integer (kind=int_kind), intent(in) :: &
1778 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
1779 : iblk ! block index
1780 :
1781 : ! local variables
1782 :
1783 : real (kind=dbl_kind) :: &
1784 : eidebug, esdebug, & ! LCOV_EXCL_LINE
1785 : qi, qs, Tsnow, si, & ! LCOV_EXCL_LINE
1786 0 : rad_to_deg, puny, rhoi, lfresh, rhos, cp_ice
1787 :
1788 : integer (kind=int_kind) :: n, k, nt_Tsfc, nt_qice, nt_qsno, nt_fsd, &
1789 : nt_isosno, nt_isoice, nt_sice, nt_smice, nt_smliq
1790 :
1791 : logical (kind=log_kind) :: tr_fsd, tr_iso, tr_snow
1792 :
1793 : type (block) :: &
1794 : this_block ! block information for current block
1795 :
1796 : character(len=*), parameter :: subname = '(print_state)'
1797 :
1798 : call icepack_query_tracer_flags(tr_fsd_out=tr_fsd, tr_iso_out=tr_iso, &
1799 0 : tr_snow_out=tr_snow)
1800 : call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, &
1801 : nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, & ! LCOV_EXCL_LINE
1802 : nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, & ! LCOV_EXCL_LINE
1803 0 : nt_smice_out=nt_smice, nt_smliq_out=nt_smliq)
1804 : call icepack_query_parameters( &
1805 : rad_to_deg_out=rad_to_deg, puny_out=puny, rhoi_out=rhoi, lfresh_out=lfresh, & ! LCOV_EXCL_LINE
1806 0 : rhos_out=rhos, cp_ice_out=cp_ice)
1807 0 : call icepack_warnings_flush(nu_diag)
1808 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1809 0 : file=__FILE__, line=__LINE__)
1810 :
1811 0 : this_block = get_block(blocks_ice(iblk),iblk)
1812 :
1813 0 : write(nu_diag,*) subname,' ',trim(plabel)
1814 0 : write(nu_diag,*) subname,' istep1, my_task, i, j, iblk:', &
1815 0 : istep1, my_task, i, j, iblk
1816 0 : write(nu_diag,*) subname,' Global block:', this_block%block_id
1817 0 : write(nu_diag,*) subname,' Global i and j:', &
1818 : this_block%i_glob(i), & ! LCOV_EXCL_LINE
1819 0 : this_block%j_glob(j)
1820 0 : write (nu_diag,*) subname,' Lat, Lon (degrees):', &
1821 : TLAT(i,j,iblk)*rad_to_deg, & ! LCOV_EXCL_LINE
1822 0 : TLON(i,j,iblk)*rad_to_deg
1823 0 : write(nu_diag,*) ' '
1824 0 : write(nu_diag,*) 'aice ', aice(i,j,iblk)
1825 0 : write(nu_diag,*) 'aice0', aice0(i,j,iblk)
1826 0 : do n = 1, ncat
1827 0 : write(nu_diag,*) ' '
1828 0 : write(nu_diag,*) 'n =',n
1829 0 : write(nu_diag,*) 'aicen', aicen(i,j,n,iblk)
1830 0 : write(nu_diag,*) 'vicen', vicen(i,j,n,iblk)
1831 0 : write(nu_diag,*) 'vsnon', vsnon(i,j,n,iblk)
1832 0 : if (aicen(i,j,n,iblk) > puny) then
1833 0 : write(nu_diag,*) 'hin', vicen(i,j,n,iblk)/aicen(i,j,n,iblk)
1834 0 : write(nu_diag,*) 'hsn', vsnon(i,j,n,iblk)/aicen(i,j,n,iblk)
1835 : endif
1836 0 : write(nu_diag,*) 'Tsfcn',trcrn(i,j,nt_Tsfc,n,iblk)
1837 0 : if (tr_fsd) write(nu_diag,*) 'afsdn',trcrn(i,j,nt_fsd,n,iblk) ! fsd cat 1
1838 : ! layer 1 diagnostics
1839 : ! if (tr_iso) write(nu_diag,*) 'isosno',trcrn(i,j,nt_isosno,n,iblk) ! isotopes in snow
1840 : ! if (tr_iso) write(nu_diag,*) 'isoice',trcrn(i,j,nt_isoice,n,iblk) ! isotopes in ice
1841 : ! if (tr_snow) write(nu_diag,*) 'smice', trcrn(i,j,nt_smice, n,iblk) ! ice mass in snow
1842 : ! if (tr_snow) write(nu_diag,*) 'smliq', trcrn(i,j,nt_smliq, n,iblk) ! liquid mass in snow
1843 0 : write(nu_diag,*) ' '
1844 :
1845 : ! dynamics (transport and/or ridging) causes the floe size distribution to become non-normal
1846 : ! if (tr_fsd) then
1847 : ! if (abs(sum(trcrn(i,j,nt_fsd:nt_fsd+nfsd-1,n,iblk))-c1) > puny) &
1848 : ! write(nu_diag,*) 'afsdn not normal', & ! LCOV_EXCL_LINE
1849 : ! sum(trcrn(i,j,nt_fsd:nt_fsd+nfsd-1,n,iblk)), & ! LCOV_EXCL_LINE
1850 : ! trcrn(i,j,nt_fsd:nt_fsd+nfsd-1,n,iblk)
1851 : ! endif
1852 :
1853 : enddo ! n
1854 :
1855 0 : eidebug = c0
1856 0 : do n = 1,ncat
1857 0 : do k = 1,nilyr
1858 0 : qi = trcrn(i,j,nt_qice+k-1,n,iblk)
1859 0 : write(nu_diag,*) 'qice, cat ',n,' layer ',k, qi
1860 0 : eidebug = eidebug + qi
1861 0 : if (aicen(i,j,n,iblk) > puny) then
1862 0 : write(nu_diag,*) 'qi/rhoi', qi/rhoi
1863 : endif
1864 : enddo
1865 0 : write(nu_diag,*) ' '
1866 : enddo
1867 0 : write(nu_diag,*) 'qice(i,j)',eidebug
1868 0 : write(nu_diag,*) ' '
1869 :
1870 0 : esdebug = c0
1871 0 : do n = 1,ncat
1872 0 : if (vsnon(i,j,n,iblk) > puny) then
1873 0 : do k = 1,nslyr
1874 0 : qs = trcrn(i,j,nt_qsno+k-1,n,iblk)
1875 0 : write(nu_diag,*) 'qsnow, cat ',n,' layer ',k, qs
1876 0 : esdebug = esdebug + qs
1877 0 : Tsnow = (Lfresh + qs/rhos) / cp_ice
1878 0 : write(nu_diag,*) 'qs/rhos', qs/rhos
1879 0 : write(nu_diag,*) 'Tsnow', Tsnow
1880 : enddo
1881 0 : write(nu_diag,*) ' '
1882 : endif
1883 : enddo
1884 0 : write(nu_diag,*) 'qsnow(i,j)',esdebug
1885 0 : write(nu_diag,*) ' '
1886 :
1887 0 : do n = 1,ncat
1888 0 : do k = 1,nilyr
1889 0 : si = trcrn(i,j,nt_sice+k-1,n,iblk)
1890 0 : write(nu_diag,*) 'sice, cat ',n,' layer ',k, si
1891 : enddo
1892 : enddo
1893 0 : write(nu_diag,*) ' '
1894 :
1895 0 : write(nu_diag,*) 'uvel(i,j)',uvel(i,j,iblk)
1896 0 : write(nu_diag,*) 'vvel(i,j)',vvel(i,j,iblk)
1897 0 : if (grid_ice == 'C') then
1898 0 : write(nu_diag,*) 'uvelE(i,j)',uvelE(i,j,iblk)
1899 0 : write(nu_diag,*) 'uvelN(i,j)',uvelN(i,j,iblk)
1900 0 : elseif (grid_ice == 'CD') then
1901 0 : write(nu_diag,*) 'uvelE(i,j)',uvelE(i,j,iblk)
1902 0 : write(nu_diag,*) 'vvelE(i,j)',vvelE(i,j,iblk)
1903 0 : write(nu_diag,*) 'uvelN(i,j)',uvelN(i,j,iblk)
1904 0 : write(nu_diag,*) 'vvelN(i,j)',vvelN(i,j,iblk)
1905 : endif
1906 :
1907 0 : write(nu_diag,*) ' '
1908 0 : write(nu_diag,*) 'atm states and fluxes'
1909 0 : write(nu_diag,*) ' uatm = ',uatm (i,j,iblk)
1910 0 : write(nu_diag,*) ' vatm = ',vatm (i,j,iblk)
1911 0 : write(nu_diag,*) ' potT = ',potT (i,j,iblk)
1912 0 : write(nu_diag,*) ' Tair = ',Tair (i,j,iblk)
1913 0 : write(nu_diag,*) ' Qa = ',Qa (i,j,iblk)
1914 0 : write(nu_diag,*) ' rhoa = ',rhoa (i,j,iblk)
1915 0 : write(nu_diag,*) ' swvdr = ',swvdr(i,j,iblk)
1916 0 : write(nu_diag,*) ' swvdf = ',swvdf(i,j,iblk)
1917 0 : write(nu_diag,*) ' swidr = ',swidr(i,j,iblk)
1918 0 : write(nu_diag,*) ' swidf = ',swidf(i,j,iblk)
1919 0 : write(nu_diag,*) ' flw = ',flw (i,j,iblk)
1920 0 : write(nu_diag,*) ' frain = ',frain(i,j,iblk)
1921 0 : write(nu_diag,*) ' fsnow = ',fsnow(i,j,iblk)
1922 0 : write(nu_diag,*) ' '
1923 0 : write(nu_diag,*) 'ocn states and fluxes'
1924 0 : write(nu_diag,*) ' frzmlt = ',frzmlt (i,j,iblk)
1925 0 : write(nu_diag,*) ' sst = ',sst (i,j,iblk)
1926 0 : write(nu_diag,*) ' sss = ',sss (i,j,iblk)
1927 0 : write(nu_diag,*) ' Tf = ',Tf (i,j,iblk)
1928 0 : write(nu_diag,*) ' uocn = ',uocn (i,j,iblk)
1929 0 : write(nu_diag,*) ' vocn = ',vocn (i,j,iblk)
1930 0 : write(nu_diag,*) ' strtltxU= ',strtltxU(i,j,iblk)
1931 0 : write(nu_diag,*) ' strtltyU= ',strtltyU(i,j,iblk)
1932 0 : write(nu_diag,*) ' '
1933 0 : write(nu_diag,*) 'srf states and fluxes'
1934 0 : write(nu_diag,*) ' Tref = ',Tref (i,j,iblk)
1935 0 : write(nu_diag,*) ' Qref = ',Qref (i,j,iblk)
1936 0 : write(nu_diag,*) ' Uref = ',Uref (i,j,iblk)
1937 0 : write(nu_diag,*) ' fsens = ',fsens (i,j,iblk)
1938 0 : write(nu_diag,*) ' flat = ',flat (i,j,iblk)
1939 0 : write(nu_diag,*) ' evap = ',evap (i,j,iblk)
1940 0 : write(nu_diag,*) ' flwout = ',flwout(i,j,iblk)
1941 0 : write(nu_diag,*) ' '
1942 0 : call flush_fileunit(nu_diag)
1943 :
1944 0 : end subroutine print_state
1945 :
1946 : !=======================================================================
1947 :
1948 : ! prints error information prior to aborting
1949 :
1950 0 : subroutine diagnostic_abort(istop, jstop, iblk, stop_label)
1951 :
1952 : use ice_blocks, only: block, get_block
1953 : use ice_domain, only: blocks_ice
1954 : use ice_grid, only: TLAT, TLON
1955 : use ice_state, only: aice
1956 :
1957 : integer (kind=int_kind), intent(in) :: &
1958 : istop, jstop, & ! indices of grid cell where model aborts ! LCOV_EXCL_LINE
1959 : iblk ! block index
1960 :
1961 : character (len=*), intent(in) :: stop_label
1962 :
1963 : ! local variables
1964 :
1965 0 : real (kind=dbl_kind) :: rad_to_deg
1966 :
1967 : type (block) :: &
1968 : this_block ! block information for current block
1969 :
1970 : character(len=*), parameter :: subname = '(diagnostic_abort)'
1971 :
1972 0 : call icepack_query_parameters(rad_to_deg_out=rad_to_deg)
1973 0 : call icepack_warnings_flush(nu_diag)
1974 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1975 0 : file=__FILE__, line=__LINE__)
1976 :
1977 0 : this_block = get_block(blocks_ice(iblk),iblk)
1978 :
1979 0 : call flush_fileunit(nu_diag)
1980 0 : if (istop > 0 .and. jstop > 0) then
1981 0 : call print_state(trim(stop_label),istop,jstop,iblk)
1982 : else
1983 0 : write (nu_diag,*) subname,' istep1, my_task, iblk =', &
1984 0 : istep1, my_task, iblk
1985 0 : write (nu_diag,*) subname,' Global block:', this_block%block_id
1986 : endif
1987 0 : call flush_fileunit(nu_diag)
1988 0 : call abort_ice (subname//'ERROR: '//trim(stop_label))
1989 :
1990 0 : end subroutine diagnostic_abort
1991 :
1992 : !=======================================================================
1993 :
1994 : end module ice_diagnostics
1995 :
1996 : !=======================================================================
|
