6 use iso_fortran_env,
only : int64
7 use mom_coms,
only : sum_across_pes, pe_here, root_pe, num_pes, max_across_pes
17 use mom_io,
only : append_file, ascii_file, single_file, writeonly_file
20 use mom_time_manager,
only : time_type, get_time, get_date, set_time,
operator(>)
21 use mom_time_manager,
only :
operator(+),
operator(-),
operator(*),
operator(/)
22 use mom_time_manager,
only :
operator(/=),
operator(<=),
operator(>=),
operator(<)
23 use mom_time_manager,
only : get_calendar_type, time_type_to_real, no_calendar
28 use mpp_mod,
only : mpp_chksum
32 implicit none ;
private
34 #include <MOM_memory.h>
65 integer,
allocatable,
dimension(:) :: lh
68 logical :: do_ape_calc
73 character(len=200) :: depth_list_file
74 real :: d_list_min_inc
76 logical :: require_depth_list_chksum
79 logical :: update_depth_list_chksum
82 logical :: use_temperature
83 real :: fresh_water_input
89 real :: net_salt_input
93 real :: net_heat_input
103 type(time_type) :: energysavedays
105 type(time_type) :: energysavedays_geometric
109 logical :: energysave_geometric
111 type(time_type) :: write_energy_time
112 type(time_type) :: geometric_end_time
117 logical :: date_stamped_output
118 type(time_type) :: start_time
120 integer,
pointer :: ntrunc => null()
126 logical :: write_stocks
128 integer :: previous_calls = 0
129 integer :: prev_n = 0
130 integer :: fileenergy_nc
131 integer :: fileenergy_ascii
132 type(fieldtype),
dimension(NUM_FIELDS+MAX_FIELDS_) :: &
134 character(len=200) :: energyfile
141 Input_start_time, CS)
146 character(len=*),
intent(in) :: directory
147 integer,
target,
intent(inout) :: ntrnc
150 type(time_type),
intent(in) :: input_start_time
158 #include "version_variable.h"
159 character(len=40) :: mdl =
"MOM_sum_output"
160 character(len=200) :: energyfile
161 character(len=32) :: filename_appendix =
''
163 if (
associated(cs))
then
164 call mom_error(warning,
"MOM_sum_output_init called with associated control structure.")
171 call get_param(param_file, mdl,
"CALCULATE_APE", cs%do_APE_calc, &
172 "If true, calculate the available potential energy of "//&
173 "the interfaces. Setting this to false reduces the "//&
174 "memory footprint of high-PE-count models dramatically.", &
176 call get_param(param_file, mdl,
"WRITE_STOCKS", cs%write_stocks, &
177 "If true, write the integrated tracer amounts to stdout "//&
178 "when the energy files are written.", default=.true.)
179 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", cs%use_temperature, &
180 "If true, Temperature and salinity are used as state "//&
181 "variables.", default=.true.)
182 call get_param(param_file, mdl,
"DT", cs%dt_in_T, &
183 "The (baroclinic) dynamics time step.", &
184 units=
"s", scale=us%s_to_T, fail_if_missing=.true.)
185 call get_param(param_file, mdl,
"MAXTRUNC", cs%maxtrunc, &
186 "The run will be stopped, and the day set to a very "//&
187 "large value if the velocity is truncated more than "//&
188 "MAXTRUNC times between energy saves. Set MAXTRUNC to 0 "//&
189 "to stop if there is any truncation of velocities.", &
190 units=
"truncations save_interval-1", default=0)
192 call get_param(param_file, mdl,
"MAX_ENERGY", cs%max_Energy, &
193 "The maximum permitted average energy per unit mass; the "//&
194 "model will be stopped if there is more energy than "//&
195 "this. If zero or negative, this is set to 10*MAXVEL^2.", &
196 units=
"m2 s-2", default=0.0)
197 if (cs%max_Energy <= 0.0)
then
198 call get_param(param_file, mdl,
"MAXVEL", maxvel, &
199 "The maximum velocity allowed before the velocity "//&
200 "components are truncated.", units=
"m s-1", default=3.0e8)
201 cs%max_Energy = 10.0 * maxvel**2
202 call log_param(param_file, mdl,
"MAX_ENERGY as used", cs%max_Energy)
205 call get_param(param_file, mdl,
"ENERGYFILE", energyfile, &
206 "The file to use to write the energies and globally "//&
207 "summed diagnostics.", default=
"ocean.stats")
210 call get_filename_appendix(filename_appendix)
211 if (len_trim(filename_appendix) > 0)
then
212 energyfile = trim(energyfile) //
'.'//trim(filename_appendix)
215 cs%energyfile = trim(slasher(directory))//trim(energyfile)
216 call log_param(param_file, mdl,
"output_path/ENERGYFILE", cs%energyfile)
218 cs%energyfile = trim(cs%energyfile)//
"."//trim(adjustl(statslabel))
221 call get_param(param_file, mdl,
"DATE_STAMPED_STDOUT", cs%date_stamped_output, &
222 "If true, use dates (not times) in messages to stdout", &
224 call get_param(param_file, mdl,
"TIMEUNIT", cs%Timeunit, &
225 "The time unit in seconds a number of input fields", &
226 units=
"s", default=86400.0)
227 if (cs%Timeunit < 0.0) cs%Timeunit = 86400.0
231 if (cs%do_APE_calc)
then
232 call get_param(param_file, mdl,
"READ_DEPTH_LIST", cs%read_depth_list, &
233 "Read the depth list from a file if it exists or "//&
234 "create that file otherwise.", default=.false.)
235 call get_param(param_file, mdl,
"DEPTH_LIST_MIN_INC", cs%D_list_min_inc, &
236 "The minimum increment between the depths of the "//&
237 "entries in the depth-list file.", &
238 units=
"m", default=1.0e-10, scale=us%m_to_Z)
239 if (cs%read_depth_list)
then
240 call get_param(param_file, mdl,
"DEPTH_LIST_FILE", cs%depth_list_file, &
241 "The name of the depth list file.", default=
"Depth_list.nc")
242 if (scan(cs%depth_list_file,
'/') == 0) &
243 cs%depth_list_file = trim(slasher(directory))//trim(cs%depth_list_file)
245 call get_param(param_file, mdl,
"REQUIRE_DEPTH_LIST_CHECKSUMS", &
246 cs%require_depth_list_chksum, &
247 "Require that matching checksums be in Depth_list.nc "//&
248 "when reading the file.", default=.true.)
249 if (.not. cs%require_depth_list_chksum) &
250 call get_param(param_file, mdl,
"UPDATE_DEPTH_LIST_CHECKSUMS", &
251 cs%update_depth_list_chksum, &
252 "Automatically update the Depth_list.nc file if the "//&
253 "checksums are missing or do not match current values.", &
257 allocate(cs%lH(g%ke))
263 call get_param(param_file, mdl,
"TIMEUNIT", time_unit, &
264 "The time unit for ENERGYSAVEDAYS.", &
265 units=
"s", default=86400.0)
266 call get_param(param_file, mdl,
"ENERGYSAVEDAYS",cs%energysavedays, &
267 "The interval in units of TIMEUNIT between saves of the "//&
268 "energies of the run and other globally summed diagnostics.",&
269 default=set_time(0,days=1), timeunit=time_unit)
270 call get_param(param_file, mdl,
"ENERGYSAVEDAYS_GEOMETRIC",cs%energysavedays_geometric, &
271 "The starting interval in units of TIMEUNIT for the first call "//&
272 "to save the energies of the run and other globally summed diagnostics. "//&
273 "The interval increases by a factor of 2. after each call to write_energy.",&
274 default=set_time(seconds=0), timeunit=time_unit)
276 if ((time_type_to_real(cs%energysavedays_geometric) > 0.) .and. &
277 (cs%energysavedays_geometric < cs%energysavedays))
then
278 cs%energysave_geometric = .true.
280 cs%energysave_geometric = .false.
283 cs%Start_time = input_start_time
292 if (
associated(cs))
then
293 if (cs%do_APE_calc)
then
294 deallocate(cs%lH, cs%DL)
303 subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_forcing)
307 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
309 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
311 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
315 type(time_type),
intent(in) :: day
316 integer,
intent(in) :: n
321 optional,
pointer :: tracer_csp
323 optional,
pointer :: obc
324 type(time_type),
optional,
intent(in) :: dt_forcing
326 real :: eta(szi_(g),szj_(g),szk_(g)+1)
327 real :: areatm(szi_(g),szj_(g))
329 real :: pe(szk_(g)+1)
332 real :: z_0ape(szk_(g)+1)
334 real :: h_0ape(szk_(g)+1)
339 real :: vol_lay(szk_(g))
341 real :: mass_lay(szk_(g))
358 real :: salin_mass_in
377 salt_efp, heat_efp, salt_chg_efp, heat_chg_efp, mass_chg_efp, &
378 mass_anom_efp, salt_anom_efp, heat_anom_efp
383 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
385 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
387 real,
dimension(SZI_(G),SZJ_(G)) :: &
390 real :: ke_scale_factor
392 integer :: num_nc_fields
394 integer :: i, j, k, is, ie, js, je, ns, nz, m, isq, ieq, jsq, jeq
395 integer :: l, lbelow, labove
398 integer :: start_of_day, num_days
400 character(len=240) :: energypath_nc
401 character(len=200) :: mesg
402 character(len=32) :: mesg_intro, time_units, day_str, n_str, date_str
403 logical :: date_stamped
404 type(time_type) :: dt_force
405 real :: tr_stocks(max_fields_)
406 real :: tr_min(max_fields_), tr_max(max_fields_)
407 real :: tr_min_x(max_fields_), tr_min_y(max_fields_), tr_min_z(max_fields_)
408 real :: tr_max_x(max_fields_), tr_max_y(max_fields_), tr_max_z(max_fields_)
409 logical :: tr_minmax_got(max_fields_) = .false.
410 character(len=40),
dimension(MAX_FIELDS_) :: &
412 integer :: ntr_stocks
413 integer :: iyear, imonth, iday, ihour, iminute, isecond, itick
414 logical :: local_open_bc
421 dt_force = set_time(seconds=2) ;
if (
present(dt_forcing)) dt_force = dt_forcing
422 if (cs%previous_calls == 0)
then
423 if (cs%energysave_geometric)
then
424 if (cs%energysavedays_geometric < cs%energysavedays)
then
425 cs%write_energy_time = day + cs%energysavedays_geometric
426 cs%geometric_end_time = cs%Start_time + cs%energysavedays * &
427 (1 + (day - cs%Start_time) / cs%energysavedays)
429 cs%write_energy_time = cs%Start_time + cs%energysavedays * &
430 (1 + (day - cs%Start_time) / cs%energysavedays)
433 cs%write_energy_time = cs%Start_time + cs%energysavedays * &
434 (1 + (day - cs%Start_time) / cs%energysavedays)
436 elseif (day + (dt_force/2) <= cs%write_energy_time)
then
439 if (cs%energysave_geometric)
then
440 if (cs%write_energy_time + cs%energysavedays_geometric >= &
441 cs%geometric_end_time)
then
442 cs%write_energy_time = cs%geometric_end_time
443 cs%energysave_geometric = .false.
445 cs%write_energy_time = cs%write_energy_time + cs%energysavedays_geometric
447 cs%energysavedays_geometric = cs%energysavedays_geometric*2
449 cs%write_energy_time = cs%write_energy_time + cs%energysavedays
454 if (.not.cs%use_temperature) num_nc_fields = 11
455 vars(1) = var_desc(
"Ntrunc",
"Nondim",
"Number of Velocity Truncations",
'1',
'1')
456 vars(2) = var_desc(
"En",
"Joules",
"Total Energy",
'1',
'1')
457 vars(3) = var_desc(
"APE",
"Joules",
"Total Interface APE",
'1',
'i')
458 vars(4) = var_desc(
"KE",
"Joules",
"Total Layer KE",
'1',
'L')
459 vars(5) = var_desc(
"H0",
"meter",
"Zero APE Depth of Interface",
'1',
'i')
460 vars(6) = var_desc(
"Mass_lay",
"kg",
"Total Layer Mass",
'1',
'L')
461 vars(7) = var_desc(
"Mass",
"kg",
"Total Mass",
'1',
'1')
462 vars(8) = var_desc(
"Mass_chg",
"kg",
"Total Mass Change between Entries",
'1',
'1')
463 vars(9) = var_desc(
"Mass_anom",
"kg",
"Anomalous Total Mass Change",
'1',
'1')
464 vars(10) = var_desc(
"max_CFL_trans",
"Nondim",
"Maximum finite-volume CFL",
'1',
'1')
465 vars(11) = var_desc(
"max_CFL_lin",
"Nondim",
"Maximum finite-difference CFL",
'1',
'1')
466 if (cs%use_temperature)
then
467 vars(12) = var_desc(
"Salt",
"kg",
"Total Salt",
'1',
'1')
468 vars(13) = var_desc(
"Salt_chg",
"kg",
"Total Salt Change between Entries",
'1',
'1')
469 vars(14) = var_desc(
"Salt_anom",
"kg",
"Anomalous Total Salt Change",
'1',
'1')
470 vars(15) = var_desc(
"Heat",
"Joules",
"Total Heat",
'1',
'1')
471 vars(16) = var_desc(
"Heat_chg",
"Joules",
"Total Heat Change between Entries",
'1',
'1')
472 vars(17) = var_desc(
"Heat_anom",
"Joules",
"Anomalous Total Heat Change",
'1',
'1')
475 local_open_bc = .false.
476 if (
present(obc))
then ;
if (
associated(obc))
then
477 local_open_bc = (obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)
480 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
481 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
482 h_to_kg_m2 = gv%H_to_kg_m2
484 if (.not.
associated(cs))
call mom_error(fatal, &
485 "write_energy: Module must be initialized before it is used.")
487 do j=js,je ;
do i=is,ie
488 areatm(i,j) = g%mask2dT(i,j)*us%L_to_m**2*g%areaT(i,j)
491 if (gv%Boussinesq)
then
493 do k=1,nz ;
do j=js,je ;
do i=is,ie
494 tmp1(i,j,k) = h(i,j,k) * (h_to_kg_m2*areatm(i,j))
495 enddo ;
enddo ;
enddo
500 if (local_open_bc)
then
501 do ns=1, obc%number_of_segments
502 segment => obc%segment(ns)
503 if (.not. segment%on_pe .or. segment%specified) cycle
504 i=segment%HI%IsdB ; j=segment%HI%JsdB
506 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
510 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
514 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
518 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
526 do k=1,nz ; vol_lay(k) = (gv%H_to_Z/h_to_kg_m2)*mass_lay(k) ;
enddo
529 if (cs%do_APE_calc)
then
530 do k=1,nz ;
do j=js,je ;
do i=is,ie
531 tmp1(i,j,k) = h_to_kg_m2 * h(i,j,k) * areatm(i,j)
532 enddo ;
enddo ;
enddo
535 call find_eta(h, tv, g, gv, us, eta)
536 do k=1,nz ;
do j=js,je ;
do i=is,ie
537 tmp1(i,j,k) = us%Z_to_m*(eta(i,j,k)-eta(i,j,k+1)) * areatm(i,j)
538 enddo ;
enddo ;
enddo
540 do k=1,nz ; vol_lay(k) = us%m_to_Z * vol_lay(k) ;
enddo
542 do k=1,nz ;
do j=js,je ;
do i=is,ie
543 tmp1(i,j,k) = h_to_kg_m2 * h(i,j,k) * areatm(i,j)
544 enddo ;
enddo ;
enddo
546 do k=1,nz ; vol_lay(k) = us%m_to_Z * (mass_lay(k) / (us%R_to_kg_m3*gv%Rho0)) ;
enddo
551 if (
present(tracer_csp))
then
553 stock_units=tr_units, num_stocks=ntr_stocks,&
554 got_min_max=tr_minmax_got, global_min=tr_min, global_max=tr_max, &
555 xgmin=tr_min_x, ygmin=tr_min_y, zgmin=tr_min_z,&
556 xgmax=tr_max_x, ygmax=tr_max_y, zgmax=tr_max_z)
557 if (ntr_stocks > 0)
then
559 vars(num_nc_fields+m) = var_desc(tr_names(m), units=tr_units(m), &
560 longname=tr_names(m), hor_grid=
'1', z_grid=
'1')
562 num_nc_fields = num_nc_fields + ntr_stocks
566 if (cs%previous_calls == 0)
then
567 cs%mass_prev = mass_tot ; cs%fresh_water_input = 0.0
569 cs%mass_prev_EFP = mass_efp
570 cs%fresh_water_in_EFP = real_to_efp(0.0)
574 if (day > cs%Start_time)
then
575 call open_file(cs%fileenergy_ascii, trim(cs%energyfile), &
576 action=append_file, form=ascii_file, nohdrs=.true.)
578 call open_file(cs%fileenergy_ascii, trim(cs%energyfile), &
579 action=writeonly_file, form=ascii_file, nohdrs=.true.)
580 if (abs(cs%timeunit - 86400.0) < 1.0)
then
581 if (cs%use_temperature)
then
582 write(cs%fileenergy_ascii,
'(" Step,",7x,"Day, Truncs, &
583 &Energy/Mass, Maximum CFL, Mean Sea Level, Total Mass, Mean Salin, &
584 &Mean Temp, Frac Mass Err, Salin Err, Temp Err")')
585 write(cs%fileenergy_ascii,
'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,&
586 &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",8x,"[degC]")')
588 write(cs%fileenergy_ascii,
'(" Step,",7x,"Day, Truncs, &
589 &Energy/Mass, Maximum CFL, Mean sea level, Total Mass, Frac Mass Err")')
590 write(cs%fileenergy_ascii,
'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,&
591 &"[kg]",11x,"[Nondim]")')
594 if ((cs%timeunit >= 0.99) .and. (cs%timeunit < 1.01))
then
595 time_units =
" [seconds] "
596 elseif ((cs%timeunit >= 3599.0) .and. (cs%timeunit < 3601.0))
then
597 time_units =
" [hours] "
598 elseif ((cs%timeunit >= 86399.0) .and. (cs%timeunit < 86401.0))
then
599 time_units =
" [days] "
600 elseif ((cs%timeunit >= 3.0e7) .and. (cs%timeunit < 3.2e7))
then
601 time_units =
" [years] "
603 write(time_units,
'(9x,"[",es8.2," s] ")') cs%timeunit
606 if (cs%use_temperature)
then
607 write(cs%fileenergy_ascii,
'(" Step,",7x,"Time, Truncs, &
608 &Energy/Mass, Maximum CFL, Mean Sea Level, Total Mass, Mean Salin, &
609 &Mean Temp, Frac Mass Err, Salin Err, Temp Err")')
610 write(cs%fileenergy_ascii,
'(A25,10x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,&
611 &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",6x,&
612 &"[degC]")') time_units
614 write(cs%fileenergy_ascii,
'(" Step,",7x,"Time, Truncs, &
615 &Energy/Mass, Maximum CFL, Mean sea level, Total Mass, Frac Mass Err")')
616 write(cs%fileenergy_ascii,
'(A25,10x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,&
617 &"[kg]",11x,"[Nondim]")') time_units
623 energypath_nc = trim(cs%energyfile) //
".nc"
624 if (day > cs%Start_time)
then
625 call reopen_file(cs%fileenergy_nc, trim(energypath_nc), vars, &
626 num_nc_fields, cs%fields, single_file, cs%timeunit, &
629 call create_file(cs%fileenergy_nc, trim(energypath_nc), vars, &
630 num_nc_fields, cs%fields, single_file, cs%timeunit, &
635 if (cs%do_APE_calc)
then
636 lbelow = 1 ; volbelow = 0.0
638 volbelow = volbelow + vol_lay(k)
639 if ((volbelow >= cs%DL(cs%lH(k))%vol_below) .and. &
640 (volbelow < cs%DL(cs%lH(k)+1)%vol_below))
then
643 labove=cs%list_size+1
644 l = (labove + lbelow) / 2
645 do while (l > lbelow)
646 if (volbelow < cs%DL(l)%vol_below)
then ; labove = l
647 else ; lbelow = l ;
endif
648 l = (labove + lbelow) / 2
653 z_0ape(k) = cs%DL(l)%depth - (volbelow - cs%DL(l)%vol_below) / cs%DL(l)%area
655 z_0ape(nz+1) = cs%DL(2)%depth
662 if (gv%Boussinesq)
then
663 do j=js,je ;
do i=is,ie
666 hbelow = hbelow + h(i,j,k) * gv%H_to_Z
667 hint = z_0ape(k) + (hbelow - g%bathyT(i,j))
668 hbot = z_0ape(k) - g%bathyT(i,j)
669 hbot = (hbot + abs(hbot)) * 0.5
670 pe_pt(i,j,k) = 0.5 * areatm(i,j) * us%Z_to_m*us%L_T_to_m_s**2*(us%R_to_kg_m3*gv%Rho0*gv%g_prime(k)) * &
671 (hint * hint - hbot * hbot)
675 do j=js,je ;
do i=is,ie
677 hint = z_0ape(k) + eta(i,j,k)
678 hbot = max(z_0ape(k) - g%bathyT(i,j), 0.0)
679 pe_pt(i,j,k) = 0.5 * (areatm(i,j) * us%Z_to_m*us%L_T_to_m_s**2*(us%R_to_kg_m3*gv%Rho0*gv%g_prime(k))) * &
680 (hint * hint - hbot * hbot)
686 do k=1,nz+1 ; h_0ape(k) = us%Z_to_m*z_0ape(k) ;
enddo
689 do k=1,nz+1 ; pe(k) = 0.0 ; h_0ape(k) = 0.0 ;
enddo
693 ke_scale_factor = gv%H_to_kg_m2*us%L_T_to_m_s**2
695 do k=1,nz ;
do j=js,je ;
do i=is,ie
696 tmp1(i,j,k) = (0.25 * ke_scale_factor * (areatm(i,j) * h(i,j,k))) * &
697 (u(i-1,j,k)**2 + u(i,j,k)**2 + v(i,j-1,k)**2 + v(i,j,k)**2)
698 enddo ;
enddo ;
enddo
701 toten = ke_tot + pe_tot
703 salt = 0.0 ; heat = 0.0
704 if (cs%use_temperature)
then
705 temp_int(:,:) = 0.0 ; salt_int(:,:) = 0.0
706 do k=1,nz ;
do j=js,je ;
do i=is,ie
707 salt_int(i,j) = salt_int(i,j) + tv%S(i,j,k) * &
708 (h(i,j,k)*(h_to_kg_m2 * areatm(i,j)))
709 temp_int(i,j) = temp_int(i,j) + (tv%C_p * tv%T(i,j,k)) * &
710 (h(i,j,k)*(h_to_kg_m2 * areatm(i,j)))
711 enddo ;
enddo ;
enddo
718 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
719 if (u(i,j,k) < 0.0)
then
720 cfl_trans = (-u(i,j,k) * cs%dt_in_T) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
722 cfl_trans = (u(i,j,k) * cs%dt_in_T) * (g%dy_Cu(i,j) * g%IareaT(i,j))
724 cfl_lin = abs(u(i,j,k) * cs%dt_in_T) * g%IdxCu(i,j)
725 max_cfl(1) = max(max_cfl(1), cfl_trans)
726 max_cfl(2) = max(max_cfl(2), cfl_lin)
727 enddo ;
enddo ;
enddo
728 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
729 if (v(i,j,k) < 0.0)
then
730 cfl_trans = (-v(i,j,k) * cs%dt_in_T) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
732 cfl_trans = (v(i,j,k) * cs%dt_in_T) * (g%dx_Cv(i,j) * g%IareaT(i,j))
734 cfl_lin = abs(v(i,j,k) * cs%dt_in_T) * g%IdyCv(i,j)
735 max_cfl(1) = max(max_cfl(1), cfl_trans)
736 max_cfl(2) = max(max_cfl(2), cfl_lin)
737 enddo ;
enddo ;
enddo
739 call sum_across_pes(cs%ntrunc)
743 if (ntr_stocks > 0)
call sum_across_pes(tr_stocks,ntr_stocks)
745 call max_across_pes(max_cfl(1))
746 call max_across_pes(max_cfl(2))
747 if (cs%use_temperature .and. cs%previous_calls == 0)
then
748 cs%salt_prev = salt ; cs%net_salt_input = 0.0
749 cs%heat_prev = heat ; cs%net_heat_input = 0.0
751 cs%salt_prev_EFP = salt_efp ; cs%net_salt_in_EFP = real_to_efp(0.0)
752 cs%heat_prev_EFP = heat_efp ; cs%net_heat_in_EFP = real_to_efp(0.0)
754 irho0 = 1.0 / (us%R_to_kg_m3*gv%Rho0)
756 if (cs%use_temperature)
then
757 salt_chg_efp = salt_efp - cs%salt_prev_EFP
758 salt_anom_efp = salt_chg_efp - cs%net_salt_in_EFP
759 salt_chg = efp_to_real(salt_chg_efp) ; salt_anom = efp_to_real(salt_anom_efp)
760 heat_chg_efp = heat_efp - cs%heat_prev_EFP
761 heat_anom_efp = heat_chg_efp - cs%net_heat_in_EFP
762 heat_chg = efp_to_real(heat_chg_efp) ; heat_anom = efp_to_real(heat_anom_efp)
765 mass_chg_efp = mass_efp - cs%mass_prev_EFP
767 if (gv%Boussinesq)
then
768 mass_anom_efp = mass_chg_efp - cs%fresh_water_in_EFP
771 mass_anom_efp = mass_chg_efp - cs%fresh_water_in_EFP
772 if (cs%use_temperature) &
773 salin_mass_in = 0.001*efp_to_real(cs%net_salt_in_EFP)
775 mass_chg = efp_to_real(mass_chg_efp)
776 mass_anom = efp_to_real(mass_anom_efp) - salin_mass_in
778 if (cs%use_temperature)
then
779 salin = salt / mass_tot ; salin_anom = salt_anom / mass_tot
781 temp = heat / (mass_tot*tv%C_p) ; temp_anom = heat_anom / (mass_tot*tv%C_p)
783 en_mass = toten / mass_tot
785 call get_time(day, start_of_day, num_days)
786 date_stamped = (cs%date_stamped_output .and. (get_calendar_type() /= no_calendar))
788 call get_date(day, iyear, imonth, iday, ihour, iminute, isecond, itick)
789 if (abs(cs%timeunit - 86400.0) < 1.0)
then
790 reday = real(num_days)+ (real(start_of_day)/86400.0)
791 mesg_intro =
"MOM Day "
793 reday = real(num_days)*(86400.0/cs%timeunit) + &
794 REAL(start_of_day)/abs(cs%timeunit)
795 mesg_intro =
"MOM Time "
797 if (reday < 1.0e8)
then ;
write(day_str,
'(F12.3)') reday
798 elseif (reday < 1.0e11)
then ;
write(day_str,
'(F15.3)') reday
799 else ;
write(day_str,
'(ES15.9)') reday ;
endif
801 if (n < 1000000)
then ;
write(n_str,
'(I6)') n
802 elseif (n < 10000000)
then ;
write(n_str,
'(I7)') n
803 elseif (n < 100000000)
then ;
write(n_str,
'(I8)') n
804 else ;
write(n_str,
'(I10)') n ;
endif
806 if (date_stamped)
then
807 write(date_str,
'("MOM Date",i7,2("/",i2.2)," ",i2.2,2(":",i2.2))') &
808 iyear, imonth, iday, ihour, iminute, isecond
810 date_str = trim(mesg_intro)//trim(day_str)
814 if (cs%use_temperature)
then
815 write(*,
'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", &
816 & ES18.12, ", Salt ", F15.11,", Temp ", F15.11)') &
817 trim(date_str), trim(n_str), en_mass, max_cfl(1), mass_tot, salin, temp
819 write(*,
'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", &
821 trim(date_str), trim(n_str), en_mass, max_cfl(1), mass_tot
824 if (cs%use_temperature)
then
825 write(cs%fileenergy_ascii,
'(A,",",A,",", I6,", En ",ES22.16, &
826 &", CFL ", F8.5, ", SL ",&
827 &es11.4,", M ",ES11.5,", S",f8.4,", T",f8.4,&
828 &", Me ",ES9.2,", Se ",ES9.2,", Te ",ES9.2)') &
829 trim(n_str), trim(day_str), cs%ntrunc, en_mass, max_cfl(1), &
830 -h_0ape(1), mass_tot, salin, temp, mass_anom/mass_tot, salin_anom, &
833 write(cs%fileenergy_ascii,
'(A,",",A,",", I6,", En ",ES22.16, &
834 &", CFL ", F8.5, ", SL ",&
835 &ES11.4,", Mass ",ES11.5,", Me ",ES9.2)') &
836 trim(n_str), trim(day_str), cs%ntrunc, en_mass, max_cfl(1), &
837 -h_0ape(1), mass_tot, mass_anom/mass_tot
840 if (cs%ntrunc > 0)
then
841 write(*,
'(A," Energy/Mass:",ES12.5," Truncations ",I0)') &
842 trim(date_str), en_mass, cs%ntrunc
845 if (cs%write_stocks)
then
846 write(*,
'(" Total Energy: ",Z16.16,ES24.16)') toten, toten
847 write(*,
'(" Total Mass: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
848 mass_tot, mass_chg, mass_anom, mass_anom/mass_tot
849 if (cs%use_temperature)
then
851 write(*,
'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') &
852 salt*0.001, salt_chg*0.001, salt_anom*0.001
854 write(*,
'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
855 salt*0.001, salt_chg*0.001, salt_anom*0.001, salt_anom/salt
858 write(*,
'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') &
859 heat, heat_chg, heat_anom
861 write(*,
'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
862 heat, heat_chg, heat_anom, heat_anom/heat
867 write(*,
'(" Total ",a,": ",ES24.16,X,a)') &
868 trim(tr_names(m)), tr_stocks(m), trim(tr_units(m))
870 if (tr_minmax_got(m))
then
871 write(*,
'(64X,"Global Min:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') &
872 tr_min(m),tr_min_x(m),tr_min_y(m),tr_min_z(m)
873 write(*,
'(64X,"Global Max:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') &
874 tr_max(m),tr_max_x(m),tr_max_y(m),tr_max_z(m)
881 var = real(cs%ntrunc)
882 call write_field(cs%fileenergy_nc, cs%fields(1), var, reday)
883 call write_field(cs%fileenergy_nc, cs%fields(2), toten, reday)
884 call write_field(cs%fileenergy_nc, cs%fields(3), pe, reday)
885 call write_field(cs%fileenergy_nc, cs%fields(4), ke, reday)
886 call write_field(cs%fileenergy_nc, cs%fields(5), h_0ape, reday)
887 call write_field(cs%fileenergy_nc, cs%fields(6), mass_lay, reday)
889 call write_field(cs%fileenergy_nc, cs%fields(7), mass_tot, reday)
890 call write_field(cs%fileenergy_nc, cs%fields(8), mass_chg, reday)
891 call write_field(cs%fileenergy_nc, cs%fields(9), mass_anom, reday)
892 call write_field(cs%fileenergy_nc, cs%fields(10), max_cfl(1), reday)
893 call write_field(cs%fileenergy_nc, cs%fields(11), max_cfl(1), reday)
894 if (cs%use_temperature)
then
895 call write_field(cs%fileenergy_nc, cs%fields(12), 0.001*salt, reday)
896 call write_field(cs%fileenergy_nc, cs%fields(13), 0.001*salt_chg, reday)
897 call write_field(cs%fileenergy_nc, cs%fields(14), 0.001*salt_anom, reday)
898 call write_field(cs%fileenergy_nc, cs%fields(15), heat, reday)
899 call write_field(cs%fileenergy_nc, cs%fields(16), heat_chg, reday)
900 call write_field(cs%fileenergy_nc, cs%fields(17), heat_anom, reday)
902 call write_field(cs%fileenergy_nc, cs%fields(17+m), tr_stocks(m), reday)
906 call write_field(cs%fileenergy_nc, cs%fields(11+m), tr_stocks(m), reday)
910 call flush_file(cs%fileenergy_nc)
913 if ((en_mass>cs%max_Energy) .or. &
914 ((en_mass>cs%max_Energy) .and. (en_mass<cs%max_Energy)))
then
915 write(mesg,
'("Energy per unit mass of ",ES11.4," exceeds ",ES11.4)') &
916 en_mass, cs%max_Energy
918 "write_energy : Excessive energy per unit mass or NaNs forced model termination.")
920 if (cs%ntrunc>cs%maxtrunc)
then
921 call mom_error(fatal,
"write_energy : Ocean velocity has been truncated too many times.")
924 cs%previous_calls = cs%previous_calls + 1
925 cs%mass_prev = mass_tot ; cs%fresh_water_input = 0.0
926 if (cs%use_temperature)
then
927 cs%salt_prev = salt ; cs%net_salt_input = 0.0
928 cs%heat_prev = heat ; cs%net_heat_input = 0.0
931 cs%mass_prev_EFP = mass_efp ; cs%fresh_water_in_EFP = real_to_efp(0.0)
932 if (cs%use_temperature)
then
933 cs%salt_prev_EFP = salt_efp ; cs%net_salt_in_EFP = real_to_efp(0.0)
934 cs%heat_prev_EFP = heat_efp ; cs%net_heat_in_EFP = real_to_efp(0.0)
943 type(
surface),
intent(in) :: sfc_state
947 real,
intent(in) :: dt
953 real,
dimension(SZI_(G),SZJ_(G)) :: &
974 integer :: i, j, is, ie, js, je
976 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
978 rzl2_to_kg = us%L_to_m**2*us%R_to_kg_m3*us%Z_to_m
980 fw_in(:,:) = 0.0 ; fw_input = 0.0
981 if (
associated(fluxes%evap))
then
982 if (
associated(fluxes%lprec) .and.
associated(fluxes%fprec))
then
983 do j=js,je ;
do i=is,ie
984 fw_in(i,j) = rzl2_to_kg * dt*g%areaT(i,j)*(fluxes%evap(i,j) + &
985 (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + fluxes%lrunoff(i,j)) + &
986 (fluxes%fprec(i,j) + fluxes%frunoff(i,j))))
990 "accumulate_net_input called with associated evap field, but no precip field.")
994 if (
associated(fluxes%seaice_melt))
then ;
do j=js,je ;
do i=is,ie
995 fw_in(i,j) = fw_in(i,j) + rzl2_to_kg*dt * &
996 g%areaT(i,j) * fluxes%seaice_melt(i,j)
997 enddo ;
enddo ;
endif
999 salt_in(:,:) = 0.0 ; heat_in(:,:) = 0.0
1000 if (cs%use_temperature)
then
1002 if (
associated(fluxes%sw))
then ;
do j=js,je ;
do i=is,ie
1003 heat_in(i,j) = heat_in(i,j) + dt*us%T_to_s*us%L_to_m**2*g%areaT(i,j) * (fluxes%sw(i,j) + &
1004 (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(i,j))))
1005 enddo ;
enddo ;
endif
1007 if (
associated(fluxes%seaice_melt_heat))
then ;
do j=js,je ;
do i=is,ie
1008 heat_in(i,j) = heat_in(i,j) + dt*us%T_to_s*us%L_to_m**2*g%areaT(i,j) * fluxes%seaice_melt_heat(i,j)
1009 enddo ;
enddo ;
endif
1022 if (
associated(tv%TempxPmE))
then
1023 do j=js,je ;
do i=is,ie
1024 heat_in(i,j) = heat_in(i,j) + (c_p * rzl2_to_kg*g%areaT(i,j)) * tv%TempxPmE(i,j)
1026 elseif (
associated(fluxes%evap))
then
1027 do j=js,je ;
do i=is,ie
1028 heat_in(i,j) = heat_in(i,j) + (c_p * sfc_state%SST(i,j)) * fw_in(i,j)
1034 if (
associated(tv%internal_heat))
then
1035 do j=js,je ;
do i=is,ie
1036 heat_in(i,j) = heat_in(i,j) + (c_p * us%L_to_m**2*g%areaT(i,j)) * &
1037 tv%internal_heat(i,j)
1040 if (
associated(tv%frazil))
then ;
do j=js,je ;
do i=is,ie
1041 heat_in(i,j) = heat_in(i,j) + us%L_to_m**2*g%areaT(i,j) * tv%frazil(i,j)
1042 enddo ;
enddo ;
endif
1043 if (
associated(fluxes%heat_added))
then ;
do j=js,je ;
do i=is,ie
1044 heat_in(i,j) = heat_in(i,j) + dt*us%T_to_s*us%L_to_m**2*g%areaT(i,j)*fluxes%heat_added(i,j)
1045 enddo ;
enddo ;
endif
1050 if (
associated(fluxes%salt_flux))
then ;
do j=js,je ;
do i=is,ie
1052 salt_in(i,j) = rzl2_to_kg * dt * &
1053 g%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j))
1054 enddo ;
enddo ;
endif
1057 if ((cs%use_temperature) .or.
associated(fluxes%lprec) .or. &
1058 associated(fluxes%evap))
then
1063 cs%fresh_water_input = cs%fresh_water_input + fw_input
1064 cs%net_salt_input = cs%net_salt_input + salt_input
1065 cs%net_heat_input = cs%net_heat_input + heat_input
1067 cs%fresh_water_in_EFP = cs%fresh_water_in_EFP + fw_in_efp
1068 cs%net_salt_in_EFP = cs%net_salt_in_EFP + salt_in_efp
1069 cs%net_heat_in_EFP = cs%net_heat_in_EFP + heat_in_efp
1086 if (cs%read_depth_list)
then
1091 trim(cs%depth_list_file)//
" does not exist. Creating a new file.")
1101 cs%lH(k) = cs%list_size
1113 real,
dimension(G%Domain%niglobal*G%Domain%njglobal + 1) :: &
1114 Dlist, & !< The global list of bottom depths [Z ~> m].
1116 integer,
dimension(G%Domain%niglobal*G%Domain%njglobal+1) :: &
1123 logical :: add_to_list
1125 integer :: ir, indxt
1126 integer :: mls, list_size
1127 integer :: list_pos, i_global, j_global
1128 integer :: i, j, k, kl
1130 mls = g%Domain%niglobal*g%Domain%njglobal
1135 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1137 j_global = j + g%jdg_offset - (g%jsg-1)
1138 i_global = i + g%idg_offset - (g%isg-1)
1140 list_pos = (j_global-1)*g%Domain%niglobal + i_global
1141 dlist(list_pos) = g%bathyT(i,j)
1142 arealist(list_pos) = g%mask2dT(i,j) * g%US%L_to_m**2*g%areaT(i,j)
1146 call sum_across_pes(dlist, mls+1)
1147 call sum_across_pes(arealist, mls+1)
1149 do j=1,mls+1 ; indx2(j) = j ;
enddo
1150 k = mls / 2 + 1 ; ir = mls
1159 indx2(ir) = indx2(1)
1161 if (ir == 1)
then ; indx2(1) = indxt ;
exit ;
endif
1164 do ;
if (j > ir)
exit
1165 if (j < ir .AND. dlist(indx2(j)) < dlist(indx2(j+1))) j = j + 1
1166 if (dnow < dlist(indx2(j)))
then ; indx2(i) = indx2(j) ; i = j ; j = j + i
1167 else ; j = ir+1 ;
endif
1179 d_list_prev = dlist(indx2(mls))
1182 if (dlist(indx2(k)) < d_list_prev-cs%D_list_min_inc)
then
1183 list_size = list_size + 1
1184 d_list_prev = dlist(indx2(k))
1188 cs%list_size = list_size
1189 allocate(cs%DL(cs%list_size+1))
1191 vol = 0.0 ; area = 0.0
1192 dprev = dlist(indx2(mls))
1198 vol = vol + area * (dprev - dlist(i))
1199 area = area + arealist(i)
1201 add_to_list = .false.
1202 if ((kl == 0) .or. (k==1))
then
1203 add_to_list = .true.
1204 elseif (dlist(indx2(k-1)) < d_list_prev-cs%D_list_min_inc)
then
1205 add_to_list = .true.
1206 d_list_prev = dlist(indx2(k-1))
1209 if (add_to_list)
then
1211 cs%DL(kl)%depth = dlist(i)
1212 cs%DL(kl)%area = area
1213 cs%DL(kl)%vol_below = vol
1218 do while (kl < list_size)
1221 cs%DL(kl)%vol_below = cs%DL(kl-1)%vol_below * 1.000001
1222 cs%DL(kl)%area = cs%DL(kl-1)%area
1223 cs%DL(kl)%depth = cs%DL(kl-1)%depth
1226 cs%DL(cs%list_size+1)%vol_below = cs%DL(cs%list_size)%vol_below * 1000.0
1227 cs%DL(cs%list_size+1)%area = cs%DL(cs%list_size)%area
1228 cs%DL(cs%list_size+1)%depth = cs%DL(cs%list_size)%depth
1238 character(len=*),
intent(in) :: filename
1239 integer,
intent(in) :: list_size
1241 real,
allocatable :: tmp(:)
1242 integer :: ncid, dimid(1), Did, Aid, Vid, status, k
1243 character(len=16) :: depth_chksum, area_chksum
1250 allocate(tmp(list_size)) ; tmp(:) = 0.0
1252 status = nf90_create(filename, 0, ncid)
1253 if (status /= nf90_noerr)
then
1254 call mom_error(warning, filename//trim(nf90_strerror(status)))
1258 status = nf90_def_dim(ncid,
"list", list_size, dimid(1))
1259 if (status /= nf90_noerr)
call mom_error(warning, &
1260 filename//trim(nf90_strerror(status)))
1262 status = nf90_def_var(ncid,
"depth", nf90_double, dimid, did)
1263 if (status /= nf90_noerr)
call mom_error(warning, &
1264 filename//
" depth "//trim(nf90_strerror(status)))
1265 status = nf90_put_att(ncid, did,
"long_name",
"Sorted depth")
1266 if (status /= nf90_noerr)
call mom_error(warning, &
1267 filename//
" depth "//trim(nf90_strerror(status)))
1268 status = nf90_put_att(ncid, did,
"units",
"m")
1269 if (status /= nf90_noerr)
call mom_error(warning, &
1270 filename//
" depth "//trim(nf90_strerror(status)))
1272 status = nf90_def_var(ncid,
"area", nf90_double, dimid, aid)
1273 if (status /= nf90_noerr)
call mom_error(warning, &
1274 filename//
" area "//trim(nf90_strerror(status)))
1275 status = nf90_put_att(ncid, aid,
"long_name",
"Open area at depth")
1276 if (status /= nf90_noerr)
call mom_error(warning, &
1277 filename//
" area "//trim(nf90_strerror(status)))
1278 status = nf90_put_att(ncid, aid,
"units",
"m2")
1279 if (status /= nf90_noerr)
call mom_error(warning, &
1280 filename//
" area "//trim(nf90_strerror(status)))
1282 status = nf90_def_var(ncid,
"vol_below", nf90_double, dimid, vid)
1283 if (status /= nf90_noerr)
call mom_error(warning, &
1284 filename//
" vol_below "//trim(nf90_strerror(status)))
1285 status = nf90_put_att(ncid, vid,
"long_name",
"Open volume below depth")
1286 if (status /= nf90_noerr)
call mom_error(warning, &
1287 filename//
" vol_below "//trim(nf90_strerror(status)))
1288 status = nf90_put_att(ncid, vid,
"units",
"m3")
1289 if (status /= nf90_noerr)
call mom_error(warning, &
1290 filename//
" vol_below "//trim(nf90_strerror(status)))
1294 if (status /= nf90_noerr)
call mom_error(warning, &
1298 if (status /= nf90_noerr)
call mom_error(warning, &
1301 status = nf90_enddef(ncid)
1302 if (status /= nf90_noerr)
call mom_error(warning, &
1303 filename//trim(nf90_strerror(status)))
1305 do k=1,list_size ; tmp(k) = us%Z_to_m*cs%DL(k)%depth ;
enddo
1306 status = nf90_put_var(ncid, did, tmp)
1307 if (status /= nf90_noerr)
call mom_error(warning, &
1308 filename//
" depth "//trim(nf90_strerror(status)))
1310 do k=1,list_size ; tmp(k) = cs%DL(k)%area ;
enddo
1311 status = nf90_put_var(ncid, aid, tmp)
1312 if (status /= nf90_noerr)
call mom_error(warning, &
1313 filename//
" area "//trim(nf90_strerror(status)))
1315 do k=1,list_size ; tmp(k) = us%Z_to_m*cs%DL(k)%vol_below ;
enddo
1316 status = nf90_put_var(ncid, vid, tmp)
1317 if (status /= nf90_noerr)
call mom_error(warning, &
1318 filename//
" vol_below "//trim(nf90_strerror(status)))
1320 status = nf90_close(ncid)
1321 if (status /= nf90_noerr)
call mom_error(warning, &
1322 filename//trim(nf90_strerror(status)))
1333 character(len=*),
intent(in) :: filename
1335 character(len=32) :: mdl
1336 character(len=240) :: var_name, var_msg
1337 real,
allocatable :: tmp(:)
1338 integer :: ncid, status, varid, list_size, k
1339 integer :: ndim, len, var_dim_ids(NF90_MAX_VAR_DIMS)
1340 character(len=16) :: depth_file_chksum, depth_grid_chksum
1341 character(len=16) :: area_file_chksum, area_grid_chksum
1342 integer :: depth_attr_status, area_attr_status
1344 mdl =
"MOM_sum_output read_depth_list:"
1346 status = nf90_open(filename, nf90_nowrite, ncid)
1347 if (status /= nf90_noerr)
then
1348 call mom_error(fatal,mdl//
" Difficulties opening "//trim(filename)// &
1349 " - "//trim(nf90_strerror(status)))
1358 if (any([depth_attr_status, area_attr_status] == nf90_enotatt))
then
1359 var_msg = trim(cs%depth_list_file) //
" checksums are missing;"
1360 if (cs%require_depth_list_chksum)
then
1361 call mom_error(fatal, trim(var_msg) //
" aborting.")
1362 elseif (cs%update_depth_list_chksum)
then
1363 call mom_error(warning, trim(var_msg) //
" updating file.")
1369 trim(var_msg) //
" some diagnostics may not be reproducible.")
1373 if (depth_attr_status /= nf90_noerr)
then
1374 var_msg = mdl //
"Failed to read " // trim(filename) //
":" &
1377 trim(var_msg) //
" - " // nf90_strerror(depth_attr_status))
1380 if (area_attr_status /= nf90_noerr)
then
1381 var_msg = mdl //
"Failed to read " // trim(filename) //
":" &
1384 trim(var_msg) //
" - " // nf90_strerror(area_attr_status))
1389 if (depth_grid_chksum /= depth_file_chksum &
1390 .or. area_grid_chksum /= area_file_chksum)
then
1391 var_msg = trim(cs%depth_list_file) //
" checksums do not match;"
1392 if (cs%require_depth_list_chksum)
then
1393 call mom_error(fatal, trim(var_msg) //
" aborting.")
1394 elseif (cs%update_depth_list_chksum)
then
1395 call mom_error(warning, trim(var_msg) //
" updating file.")
1401 trim(var_msg) //
" some diagnostics may not be reproducible.")
1407 var_msg = trim(var_name)//
" in "//trim(filename)//
" - "
1408 status = nf90_inq_varid(ncid, var_name, varid)
1409 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1410 " Difficulties finding variable "//trim(var_msg)//&
1411 trim(nf90_strerror(status)))
1413 status = nf90_inquire_variable(ncid, varid, ndims=ndim, dimids=var_dim_ids)
1414 if (status /= nf90_noerr)
then
1415 call mom_error(fatal,mdl//
" cannot inquire about "//trim(var_msg)//&
1416 trim(nf90_strerror(status)))
1417 elseif (ndim > 1)
then
1418 call mom_error(fatal,mdl//
" "//trim(var_msg)//&
1419 " has too many or too few dimensions.")
1423 status = nf90_inquire_dimension(ncid, var_dim_ids(1), len=list_size)
1424 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1425 " cannot inquire about dimension(1) of "//trim(var_msg)//&
1426 trim(nf90_strerror(status)))
1428 cs%list_size = list_size-1
1429 allocate(cs%DL(list_size))
1430 allocate(tmp(list_size))
1432 status = nf90_get_var(ncid, varid, tmp)
1433 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1434 " Difficulties reading variable "//trim(var_msg)//&
1435 trim(nf90_strerror(status)))
1437 do k=1,list_size ; cs%DL(k)%depth = us%m_to_Z*tmp(k) ;
enddo
1440 var_msg = trim(var_name)//
" in "//trim(filename)//
" - "
1441 status = nf90_inq_varid(ncid, var_name, varid)
1442 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1443 " Difficulties finding variable "//trim(var_msg)//&
1444 trim(nf90_strerror(status)))
1445 status = nf90_get_var(ncid, varid, tmp)
1446 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1447 " Difficulties reading variable "//trim(var_msg)//&
1448 trim(nf90_strerror(status)))
1450 do k=1,list_size ; cs%DL(k)%area = tmp(k) ;
enddo
1452 var_name =
"vol_below"
1453 var_msg = trim(var_name)//
" in "//trim(filename)
1454 status = nf90_inq_varid(ncid, var_name, varid)
1455 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1456 " Difficulties finding variable "//trim(var_msg)//&
1457 trim(nf90_strerror(status)))
1458 status = nf90_get_var(ncid, varid, tmp)
1459 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1460 " Difficulties reading variable "//trim(var_msg)//&
1461 trim(nf90_strerror(status)))
1463 do k=1,list_size ; cs%DL(k)%vol_below = us%m_to_Z*tmp(k) ;
enddo
1465 status = nf90_close(ncid)
1466 if (status /= nf90_noerr)
call mom_error(warning, mdl// &
1467 " Difficulties closing "//trim(filename)//
" - "//trim(nf90_strerror(status)))
1486 character(len=16),
intent(out) :: depth_chksum
1487 character(len=16),
intent(out) :: area_chksum
1490 real,
allocatable :: field(:,:)
1492 allocate(field(g%isc:g%iec, g%jsc:g%jec))
1495 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1496 field(i,j) = g%bathyT(i,j)
1498 write(depth_chksum,
'(Z16)') mpp_chksum(field(:,:))
1501 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1502 field(i,j) = g%mask2dT(i,j) * g%US%L_to_m**2*g%areaT(i,j)
1504 write(area_chksum,
'(Z16)') mpp_chksum(field(:,:))