6 use iso_fortran_env,
only : int64
7 use mom_coms,
only : sum_across_pes, pe_here, root_pe, num_pes, max_across_pes
15 use mom_io,
only : create_file, fieldtype, flush_file, open_file, reopen_file
17 use mom_io,
only : append_file, ascii_file, single_file, writeonly_file
19 use mom_open_boundary,
only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s
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>
36 public write_energy, accumulate_net_input, mom_sum_output_init
43 integer,
parameter :: num_fields = 17
44 character (*),
parameter :: depth_chksum_attr =
"bathyT_checksum"
47 character (*),
parameter :: area_chksum_attr =
"mask2dT_areaT_checksum"
65 integer,
allocatable,
dimension(:) :: lh
68 logical :: do_ape_calc
71 logical :: read_depth_list
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
140 subroutine mom_sum_output_init(G, US, param_file, directory, ntrnc, &
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, &
183 "The (baroclinic) dynamics time step.", units=
"s", &
184 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))
258 call depth_list_setup(g, us, cs)
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
286 end subroutine mom_sum_output_init
289 subroutine mom_sum_output_end(CS)
292 if (
associated(cs))
then
293 if (cs%do_APE_calc)
then
294 deallocate(cs%lH, cs%DL)
299 end subroutine mom_sum_output_end
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 integer :: num_nc_fields
392 integer :: i, j, k, is, ie, js, je, ns, nz, m, isq, ieq, jsq, jeq
393 integer :: l, lbelow, labove
396 integer :: start_of_day, num_days
398 character(len=240) :: energypath_nc
399 character(len=200) :: mesg
400 character(len=32) :: mesg_intro, time_units, day_str, n_str, date_str
401 logical :: date_stamped
402 type(time_type) :: dt_force
403 real :: tr_stocks(max_fields_)
404 real :: tr_min(max_fields_), tr_max(max_fields_)
405 real :: tr_min_x(max_fields_), tr_min_y(max_fields_), tr_min_z(max_fields_)
406 real :: tr_max_x(max_fields_), tr_max_y(max_fields_), tr_max_z(max_fields_)
407 logical :: tr_minmax_got(max_fields_) = .false.
408 character(len=40),
dimension(MAX_FIELDS_) :: &
410 integer :: ntr_stocks
411 integer :: iyear, imonth, iday, ihour, iminute, isecond, itick
412 logical :: local_open_bc
416 type(
vardesc) :: vars(num_fields+max_fields_)
419 dt_force = set_time(seconds=2) ;
if (
present(dt_forcing)) dt_force = dt_forcing
420 if (cs%previous_calls == 0)
then
421 if (cs%energysave_geometric)
then
422 if (cs%energysavedays_geometric < cs%energysavedays)
then
423 cs%write_energy_time = day + cs%energysavedays_geometric
424 cs%geometric_end_time = cs%Start_time + cs%energysavedays * &
425 (1 + (day - cs%Start_time) / cs%energysavedays)
427 cs%write_energy_time = cs%Start_time + cs%energysavedays * &
428 (1 + (day - cs%Start_time) / cs%energysavedays)
431 cs%write_energy_time = cs%Start_time + cs%energysavedays * &
432 (1 + (day - cs%Start_time) / cs%energysavedays)
434 elseif (day + (dt_force/2) <= cs%write_energy_time)
then
437 if (cs%energysave_geometric)
then
438 if (cs%write_energy_time + cs%energysavedays_geometric >= &
439 cs%geometric_end_time)
then
440 cs%write_energy_time = cs%geometric_end_time
441 cs%energysave_geometric = .false.
443 cs%write_energy_time = cs%write_energy_time + cs%energysavedays_geometric
445 cs%energysavedays_geometric = cs%energysavedays_geometric*2
447 cs%write_energy_time = cs%write_energy_time + cs%energysavedays
452 if (.not.cs%use_temperature) num_nc_fields = 11
453 vars(1) = var_desc(
"Ntrunc",
"Nondim",
"Number of Velocity Truncations",
'1',
'1')
454 vars(2) = var_desc(
"En",
"Joules",
"Total Energy",
'1',
'1')
455 vars(3) = var_desc(
"APE",
"Joules",
"Total Interface APE",
'1',
'i')
456 vars(4) = var_desc(
"KE",
"Joules",
"Total Layer KE",
'1',
'L')
457 vars(5) = var_desc(
"H0",
"meter",
"Zero APE Depth of Interface",
'1',
'i')
458 vars(6) = var_desc(
"Mass_lay",
"kg",
"Total Layer Mass",
'1',
'L')
459 vars(7) = var_desc(
"Mass",
"kg",
"Total Mass",
'1',
'1')
460 vars(8) = var_desc(
"Mass_chg",
"kg",
"Total Mass Change between Entries",
'1',
'1')
461 vars(9) = var_desc(
"Mass_anom",
"kg",
"Anomalous Total Mass Change",
'1',
'1')
462 vars(10) = var_desc(
"max_CFL_trans",
"Nondim",
"Maximum finite-volume CFL",
'1',
'1')
463 vars(11) = var_desc(
"max_CFL_lin",
"Nondim",
"Maximum finite-difference CFL",
'1',
'1')
464 if (cs%use_temperature)
then
465 vars(12) = var_desc(
"Salt",
"kg",
"Total Salt",
'1',
'1')
466 vars(13) = var_desc(
"Salt_chg",
"kg",
"Total Salt Change between Entries",
'1',
'1')
467 vars(14) = var_desc(
"Salt_anom",
"kg",
"Anomalous Total Salt Change",
'1',
'1')
468 vars(15) = var_desc(
"Heat",
"Joules",
"Total Heat",
'1',
'1')
469 vars(16) = var_desc(
"Heat_chg",
"Joules",
"Total Heat Change between Entries",
'1',
'1')
470 vars(17) = var_desc(
"Heat_anom",
"Joules",
"Anomalous Total Heat Change",
'1',
'1')
473 local_open_bc = .false.
474 if (
present(obc))
then ;
if (
associated(obc))
then
475 local_open_bc = (obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)
478 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
479 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
480 h_to_kg_m2 = gv%H_to_kg_m2
482 if (.not.
associated(cs))
call mom_error(fatal, &
483 "write_energy: Module must be initialized before it is used.")
485 do j=js,je ;
do i=is,ie
486 areatm(i,j) = g%mask2dT(i,j)*g%areaT(i,j)
489 if (gv%Boussinesq)
then
491 do k=1,nz ;
do j=js,je ;
do i=is,ie
492 tmp1(i,j,k) = h(i,j,k) * (h_to_kg_m2*areatm(i,j))
493 enddo ;
enddo ;
enddo
498 if (local_open_bc)
then
499 do ns=1, obc%number_of_segments
500 segment => obc%segment(ns)
501 if (.not. segment%on_pe .or. segment%specified) cycle
502 i=segment%HI%IsdB ; j=segment%HI%JsdB
503 if (segment%direction == obc_direction_e)
then
504 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
507 elseif (segment%direction == obc_direction_w)
then
508 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
511 elseif (segment%direction == obc_direction_n)
then
512 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
515 elseif (segment%direction == obc_direction_s)
then
516 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
524 do k=1,nz ; vol_lay(k) = (gv%H_to_Z/h_to_kg_m2)*mass_lay(k) ;
enddo
527 if (cs%do_APE_calc)
then
528 do k=1,nz ;
do j=js,je ;
do i=is,ie
529 tmp1(i,j,k) = h_to_kg_m2 * h(i,j,k) * areatm(i,j)
530 enddo ;
enddo ;
enddo
533 call find_eta(h, tv, g, gv, us, eta)
534 do k=1,nz ;
do j=js,je ;
do i=is,ie
535 tmp1(i,j,k) = (eta(i,j,k)-eta(i,j,k+1)) * areatm(i,j)
536 enddo ;
enddo ;
enddo
539 do k=1,nz ;
do j=js,je ;
do i=is,ie
540 tmp1(i,j,k) = h_to_kg_m2 * h(i,j,k) * areatm(i,j)
541 enddo ;
enddo ;
enddo
543 do k=1,nz ; vol_lay(k) = us%m_to_Z * (mass_lay(k) / gv%Rho0) ;
enddo
548 if (
present(tracer_csp))
then
549 call call_tracer_stocks(h, tr_stocks, g, gv, tracer_csp, stock_names=tr_names, &
550 stock_units=tr_units, num_stocks=ntr_stocks,&
551 got_min_max=tr_minmax_got, global_min=tr_min, global_max=tr_max, &
552 xgmin=tr_min_x, ygmin=tr_min_y, zgmin=tr_min_z,&
553 xgmax=tr_max_x, ygmax=tr_max_y, zgmax=tr_max_z)
554 if (ntr_stocks > 0)
then
556 vars(num_nc_fields+m) = var_desc(tr_names(m), units=tr_units(m), &
557 longname=tr_names(m), hor_grid=
'1', z_grid=
'1')
559 num_nc_fields = num_nc_fields + ntr_stocks
563 if (cs%previous_calls == 0)
then
564 cs%mass_prev = mass_tot ; cs%fresh_water_input = 0.0
566 cs%mass_prev_EFP = mass_efp
567 cs%fresh_water_in_EFP = real_to_efp(0.0)
570 if (is_root_pe())
then
571 if (day > cs%Start_time)
then
572 call open_file(cs%fileenergy_ascii, trim(cs%energyfile), &
573 action=append_file, form=ascii_file, nohdrs=.true.)
575 call open_file(cs%fileenergy_ascii, trim(cs%energyfile), &
576 action=writeonly_file, form=ascii_file, nohdrs=.true.)
577 if (abs(cs%timeunit - 86400.0) < 1.0)
then
578 if (cs%use_temperature)
then
579 write(cs%fileenergy_ascii,
'(" Step,",7x,"Day, Truncs, &
580 &Energy/Mass, Maximum CFL, Mean Sea Level, Total Mass, Mean Salin, &
581 &Mean Temp, Frac Mass Err, Salin Err, Temp Err")')
582 write(cs%fileenergy_ascii,
'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,&
583 &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",8x,"[degC]")')
585 write(cs%fileenergy_ascii,
'(" Step,",7x,"Day, Truncs, &
586 &Energy/Mass, Maximum CFL, Mean sea level, Total Mass, Frac Mass Err")')
587 write(cs%fileenergy_ascii,
'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,&
588 &"[kg]",11x,"[Nondim]")')
591 if ((cs%timeunit >= 0.99) .and. (cs%timeunit < 1.01))
then
592 time_units =
" [seconds] "
593 elseif ((cs%timeunit >= 3599.0) .and. (cs%timeunit < 3601.0))
then
594 time_units =
" [hours] "
595 elseif ((cs%timeunit >= 86399.0) .and. (cs%timeunit < 86401.0))
then
596 time_units =
" [days] "
597 elseif ((cs%timeunit >= 3.0e7) .and. (cs%timeunit < 3.2e7))
then
598 time_units =
" [years] "
600 write(time_units,
'(9x,"[",es8.2," s] ")') cs%timeunit
603 if (cs%use_temperature)
then
604 write(cs%fileenergy_ascii,
'(" Step,",7x,"Time, Truncs, &
605 &Energy/Mass, Maximum CFL, Mean Sea Level, Total Mass, Mean Salin, &
606 &Mean Temp, Frac Mass Err, Salin Err, Temp Err")')
607 write(cs%fileenergy_ascii,
'(A25,10x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,&
608 &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",6x,&
609 &"[degC]")') time_units
611 write(cs%fileenergy_ascii,
'(" Step,",7x,"Time, Truncs, &
612 &Energy/Mass, Maximum CFL, Mean sea level, Total Mass, Frac Mass Err")')
613 write(cs%fileenergy_ascii,
'(A25,10x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,&
614 &"[kg]",11x,"[Nondim]")') time_units
620 energypath_nc = trim(cs%energyfile) //
".nc"
621 if (day > cs%Start_time)
then
622 call reopen_file(cs%fileenergy_nc, trim(energypath_nc), vars, &
623 num_nc_fields, cs%fields, single_file, cs%timeunit, &
626 call create_file(cs%fileenergy_nc, trim(energypath_nc), vars, &
627 num_nc_fields, cs%fields, single_file, cs%timeunit, &
632 if (cs%do_APE_calc)
then
633 lbelow = 1 ; volbelow = 0.0
635 volbelow = volbelow + vol_lay(k)
636 if ((volbelow >= cs%DL(cs%lH(k))%vol_below) .and. &
637 (volbelow < cs%DL(cs%lH(k)+1)%vol_below))
then
640 labove=cs%list_size+1
641 l = (labove + lbelow) / 2
642 do while (l > lbelow)
643 if (volbelow < cs%DL(l)%vol_below)
then ; labove = l
644 else ; lbelow = l ;
endif
645 l = (labove + lbelow) / 2
650 z_0ape(k) = cs%DL(l)%depth - (volbelow - cs%DL(l)%vol_below) / cs%DL(l)%area
652 z_0ape(nz+1) = cs%DL(2)%depth
659 if (gv%Boussinesq)
then
660 do j=js,je ;
do i=is,ie
663 hbelow = hbelow + h(i,j,k) * gv%H_to_Z
664 hint = z_0ape(k) + (hbelow - g%bathyT(i,j))
665 hbot = z_0ape(k) - g%bathyT(i,j)
666 hbot = (hbot + abs(hbot)) * 0.5
667 pe_pt(i,j,k) = 0.5 * areatm(i,j) * us%Z_to_m*(gv%Rho0*us%L_to_m**2*us%s_to_T**2*gv%g_prime(k)) * &
668 (hint * hint - hbot * hbot)
672 do j=js,je ;
do i=is,ie
674 hint = z_0ape(k) + eta(i,j,k)
675 hbot = max(z_0ape(k) - g%bathyT(i,j), 0.0)
676 pe_pt(i,j,k) = 0.5 * (areatm(i,j) * us%Z_to_m*(gv%Rho0*us%L_to_m**2*us%s_to_T**2*gv%g_prime(k))) * &
677 (hint * hint - hbot * hbot)
683 do k=1,nz+1 ; h_0ape(k) = us%Z_to_m*z_0ape(k) ;
enddo
686 do k=1,nz+1 ; pe(k) = 0.0 ; h_0ape(k) = 0.0 ;
enddo
691 do k=1,nz ;
do j=js,je ;
do i=is,ie
692 tmp1(i,j,k) = (0.25 * h_to_kg_m2 * (areatm(i,j) * h(i,j,k))) * &
693 (u(i-1,j,k)**2 + u(i,j,k)**2 + v(i,j-1,k)**2 + v(i,j,k)**2)
694 enddo ;
enddo ;
enddo
697 toten = ke_tot + pe_tot
699 salt = 0.0 ; heat = 0.0
700 if (cs%use_temperature)
then
701 temp_int(:,:) = 0.0 ; salt_int(:,:) = 0.0
702 do k=1,nz ;
do j=js,je ;
do i=is,ie
703 salt_int(i,j) = salt_int(i,j) + tv%S(i,j,k) * &
704 (h(i,j,k)*(h_to_kg_m2 * areatm(i,j)))
705 temp_int(i,j) = temp_int(i,j) + (tv%C_p * tv%T(i,j,k)) * &
706 (h(i,j,k)*(h_to_kg_m2 * areatm(i,j)))
707 enddo ;
enddo ;
enddo
714 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
715 if (u(i,j,k) < 0.0)
then
716 cfl_trans = (-u(i,j,k) * cs%dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
718 cfl_trans = (u(i,j,k) * cs%dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
720 cfl_lin = abs(u(i,j,k) * cs%dt) * g%IdxCu(i,j)
721 max_cfl(1) = max(max_cfl(1), cfl_trans)
722 max_cfl(2) = max(max_cfl(2), cfl_lin)
723 enddo ;
enddo ;
enddo
724 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
725 if (v(i,j,k) < 0.0)
then
726 cfl_trans = (-v(i,j,k) * cs%dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
728 cfl_trans = (v(i,j,k) * cs%dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
730 cfl_lin = abs(v(i,j,k) * cs%dt) * g%IdyCv(i,j)
731 max_cfl(1) = max(max_cfl(1), cfl_trans)
732 max_cfl(2) = max(max_cfl(2), cfl_lin)
733 enddo ;
enddo ;
enddo
735 call sum_across_pes(cs%ntrunc)
739 if (ntr_stocks > 0)
call sum_across_pes(tr_stocks,ntr_stocks)
741 call max_across_pes(max_cfl(1))
742 call max_across_pes(max_cfl(2))
743 if (cs%use_temperature .and. cs%previous_calls == 0)
then
744 cs%salt_prev = salt ; cs%net_salt_input = 0.0
745 cs%heat_prev = heat ; cs%net_heat_input = 0.0
747 cs%salt_prev_EFP = salt_efp ; cs%net_salt_in_EFP = real_to_efp(0.0)
748 cs%heat_prev_EFP = heat_efp ; cs%net_heat_in_EFP = real_to_efp(0.0)
752 if (cs%use_temperature)
then
753 salt_chg_efp = salt_efp - cs%salt_prev_EFP
754 salt_anom_efp = salt_chg_efp - cs%net_salt_in_EFP
755 salt_chg = efp_to_real(salt_chg_efp) ; salt_anom = efp_to_real(salt_anom_efp)
756 heat_chg_efp = heat_efp - cs%heat_prev_EFP
757 heat_anom_efp = heat_chg_efp - cs%net_heat_in_EFP
758 heat_chg = efp_to_real(heat_chg_efp) ; heat_anom = efp_to_real(heat_anom_efp)
761 mass_chg_efp = mass_efp - cs%mass_prev_EFP
763 if (gv%Boussinesq)
then
764 mass_anom_efp = mass_chg_efp - cs%fresh_water_in_EFP
767 mass_anom_efp = mass_chg_efp - cs%fresh_water_in_EFP
768 if (cs%use_temperature) &
769 salin_mass_in = 0.001*efp_to_real(cs%net_salt_in_EFP)
771 mass_chg = efp_to_real(mass_chg_efp)
772 mass_anom = efp_to_real(mass_anom_efp) - salin_mass_in
774 if (cs%use_temperature)
then
775 salin = salt / mass_tot ; salin_anom = salt_anom / mass_tot
777 temp = heat / (mass_tot*tv%C_p) ; temp_anom = heat_anom / (mass_tot*tv%C_p)
779 en_mass = toten / mass_tot
781 call get_time(day, start_of_day, num_days)
782 date_stamped = (cs%date_stamped_output .and. (get_calendar_type() /= no_calendar))
784 call get_date(day, iyear, imonth, iday, ihour, iminute, isecond, itick)
785 if (abs(cs%timeunit - 86400.0) < 1.0)
then
786 reday = real(num_days)+ (real(start_of_day)/86400.0)
787 mesg_intro =
"MOM Day "
789 reday = real(num_days)*(86400.0/cs%timeunit) + &
790 REAL(start_of_day)/abs(cs%timeunit)
791 mesg_intro =
"MOM Time "
793 if (reday < 1.0e8)
then ;
write(day_str,
'(F12.3)') reday
794 elseif (reday < 1.0e11)
then ;
write(day_str,
'(F15.3)') reday
795 else ;
write(day_str,
'(ES15.9)') reday ;
endif
797 if (n < 1000000)
then ;
write(n_str,
'(I6)') n
798 elseif (n < 10000000)
then ;
write(n_str,
'(I7)') n
799 elseif (n < 100000000)
then ;
write(n_str,
'(I8)') n
800 else ;
write(n_str,
'(I10)') n ;
endif
802 if (date_stamped)
then
803 write(date_str,
'("MOM Date",i7,2("/",i2.2)," ",i2.2,2(":",i2.2))') &
804 iyear, imonth, iday, ihour, iminute, isecond
806 date_str = trim(mesg_intro)//trim(day_str)
809 if (is_root_pe())
then
810 if (cs%use_temperature)
then
811 write(*,
'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", &
812 & ES18.12, ", Salt ", F15.11,", Temp ", F15.11)') &
813 trim(date_str), trim(n_str), en_mass, max_cfl(1), mass_tot, salin, temp
815 write(*,
'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", &
817 trim(date_str), trim(n_str), en_mass, max_cfl(1), mass_tot
820 if (cs%use_temperature)
then
821 write(cs%fileenergy_ascii,
'(A,",",A,",", I6,", En ",ES22.16, &
822 &", CFL ", F8.5, ", SL ",&
823 &es11.4,", M ",ES11.5,", S",f8.4,", T",f8.4,&
824 &", Me ",ES9.2,", Se ",ES9.2,", Te ",ES9.2)') &
825 trim(n_str), trim(day_str), cs%ntrunc, en_mass, max_cfl(1), &
826 -h_0ape(1), mass_tot, salin, temp, mass_anom/mass_tot, salin_anom, &
829 write(cs%fileenergy_ascii,
'(A,",",A,",", I6,", En ",ES22.16, &
830 &", CFL ", F8.5, ", SL ",&
831 &ES11.4,", Mass ",ES11.5,", Me ",ES9.2)') &
832 trim(n_str), trim(day_str), cs%ntrunc, en_mass, max_cfl(1), &
833 -h_0ape(1), mass_tot, mass_anom/mass_tot
836 if (cs%ntrunc > 0)
then
837 write(*,
'(A," Energy/Mass:",ES12.5," Truncations ",I0)') &
838 trim(date_str), en_mass, cs%ntrunc
841 if (cs%write_stocks)
then
842 write(*,
'(" Total Energy: ",Z16.16,ES24.16)') toten, toten
843 write(*,
'(" Total Mass: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
844 mass_tot, mass_chg, mass_anom, mass_anom/mass_tot
845 if (cs%use_temperature)
then
847 write(*,
'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') &
848 salt*0.001, salt_chg*0.001, salt_anom*0.001
850 write(*,
'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
851 salt*0.001, salt_chg*0.001, salt_anom*0.001, salt_anom/salt
854 write(*,
'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') &
855 heat, heat_chg, heat_anom
857 write(*,
'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
858 heat, heat_chg, heat_anom, heat_anom/heat
863 write(*,
'(" Total ",a,": ",ES24.16,X,a)') &
864 trim(tr_names(m)), tr_stocks(m), trim(tr_units(m))
866 if (tr_minmax_got(m))
then
867 write(*,
'(64X,"Global Min:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') &
868 tr_min(m),tr_min_x(m),tr_min_y(m),tr_min_z(m)
869 write(*,
'(64X,"Global Max:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') &
870 tr_max(m),tr_max_x(m),tr_max_y(m),tr_max_z(m)
877 var = real(cs%ntrunc)
878 call write_field(cs%fileenergy_nc, cs%fields(1), var, reday)
879 call write_field(cs%fileenergy_nc, cs%fields(2), toten, reday)
880 call write_field(cs%fileenergy_nc, cs%fields(3), pe, reday)
881 call write_field(cs%fileenergy_nc, cs%fields(4), ke, reday)
882 call write_field(cs%fileenergy_nc, cs%fields(5), h_0ape, reday)
883 call write_field(cs%fileenergy_nc, cs%fields(6), mass_lay, reday)
885 call write_field(cs%fileenergy_nc, cs%fields(7), mass_tot, reday)
886 call write_field(cs%fileenergy_nc, cs%fields(8), mass_chg, reday)
887 call write_field(cs%fileenergy_nc, cs%fields(9), mass_anom, reday)
888 call write_field(cs%fileenergy_nc, cs%fields(10), max_cfl(1), reday)
889 call write_field(cs%fileenergy_nc, cs%fields(11), max_cfl(1), reday)
890 if (cs%use_temperature)
then
891 call write_field(cs%fileenergy_nc, cs%fields(12), 0.001*salt, reday)
892 call write_field(cs%fileenergy_nc, cs%fields(13), 0.001*salt_chg, reday)
893 call write_field(cs%fileenergy_nc, cs%fields(14), 0.001*salt_anom, reday)
894 call write_field(cs%fileenergy_nc, cs%fields(15), heat, reday)
895 call write_field(cs%fileenergy_nc, cs%fields(16), heat_chg, reday)
896 call write_field(cs%fileenergy_nc, cs%fields(17), heat_anom, reday)
898 call write_field(cs%fileenergy_nc, cs%fields(17+m), tr_stocks(m), reday)
902 call write_field(cs%fileenergy_nc, cs%fields(11+m), tr_stocks(m), reday)
906 call flush_file(cs%fileenergy_nc)
909 if ((en_mass>cs%max_Energy) .or. &
910 ((en_mass>cs%max_Energy) .and. (en_mass<cs%max_Energy)))
then
911 write(mesg,
'("Energy per unit mass of ",ES11.4," exceeds ",ES11.4)') &
912 en_mass, cs%max_Energy
913 call mom_error(fatal, &
914 "write_energy : Excessive energy per unit mass or NaNs forced model termination.")
916 if (cs%ntrunc>cs%maxtrunc)
then
917 call mom_error(fatal,
"write_energy : Ocean velocity has been truncated too many times.")
920 cs%previous_calls = cs%previous_calls + 1
921 cs%mass_prev = mass_tot ; cs%fresh_water_input = 0.0
922 if (cs%use_temperature)
then
923 cs%salt_prev = salt ; cs%net_salt_input = 0.0
924 cs%heat_prev = heat ; cs%net_heat_input = 0.0
927 cs%mass_prev_EFP = mass_efp ; cs%fresh_water_in_EFP = real_to_efp(0.0)
928 if (cs%use_temperature)
then
929 cs%salt_prev_EFP = salt_efp ; cs%net_salt_in_EFP = real_to_efp(0.0)
930 cs%heat_prev_EFP = heat_efp ; cs%net_heat_in_EFP = real_to_efp(0.0)
932 end subroutine write_energy
936 subroutine accumulate_net_input(fluxes, sfc_state, dt, G, CS)
939 type(
surface),
intent(in) :: sfc_state
941 real,
intent(in) :: dt
946 real,
dimension(SZI_(G),SZJ_(G)) :: &
966 integer :: i, j, is, ie, js, je
968 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
971 fw_in(:,:) = 0.0 ; fw_input = 0.0
972 if (
associated(fluxes%evap))
then
973 if (
associated(fluxes%lprec) .and.
associated(fluxes%fprec))
then
974 do j=js,je ;
do i=is,ie
975 fw_in(i,j) = dt*g%areaT(i,j)*(fluxes%evap(i,j) + &
976 (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + fluxes%lrunoff(i,j)) + &
977 (fluxes%fprec(i,j) + fluxes%frunoff(i,j))))
980 call mom_error(warning, &
981 "accumulate_net_input called with associated evap field, but no precip field.")
985 if (
associated(fluxes%seaice_melt))
then ;
do j=js,je ;
do i=is,ie
986 fw_in(i,j) = fw_in(i,j) + dt * g%areaT(i,j) * fluxes%seaice_melt(i,j)
987 enddo ;
enddo ;
endif
989 salt_in(:,:) = 0.0 ; heat_in(:,:) = 0.0
990 if (cs%use_temperature)
then
992 if (
associated(fluxes%sw))
then ;
do j=js,je ;
do i=is,ie
993 heat_in(i,j) = heat_in(i,j) + dt*g%areaT(i,j) * (fluxes%sw(i,j) + &
994 (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(i,j))))
995 enddo ;
enddo ;
endif
997 if (
associated(fluxes%seaice_melt_heat))
then ;
do j=js,je ;
do i=is,ie
998 heat_in(i,j) = heat_in(i,j) + dt*g%areaT(i,j) * fluxes%seaice_melt_heat(i,j)
999 enddo ;
enddo ;
endif
1012 if (
associated(sfc_state%TempxPmE))
then
1013 do j=js,je ;
do i=is,ie
1014 heat_in(i,j) = heat_in(i,j) + (c_p * g%areaT(i,j)) * sfc_state%TempxPmE(i,j)
1016 elseif (
associated(fluxes%evap))
then
1017 do j=js,je ;
do i=is,ie
1018 heat_in(i,j) = heat_in(i,j) + (c_p * sfc_state%SST(i,j)) * fw_in(i,j)
1024 if (
associated(sfc_state%internal_heat))
then
1025 do j=js,je ;
do i=is,ie
1026 heat_in(i,j) = heat_in(i,j) + (c_p * g%areaT(i,j)) * &
1027 sfc_state%internal_heat(i,j)
1030 if (
associated(sfc_state%frazil))
then ;
do j=js,je ;
do i=is,ie
1031 heat_in(i,j) = heat_in(i,j) + g%areaT(i,j) * sfc_state%frazil(i,j)
1032 enddo ;
enddo ;
endif
1033 if (
associated(fluxes%heat_added))
then ;
do j=js,je ;
do i=is,ie
1034 heat_in(i,j) = heat_in(i,j) + dt*g%areaT(i,j)*fluxes%heat_added(i,j)
1035 enddo ;
enddo ;
endif
1040 if (
associated(fluxes%salt_flux))
then ;
do j=js,je ;
do i=is,ie
1042 salt_in(i,j) = dt*g%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j))
1043 enddo ;
enddo ;
endif
1046 if ((cs%use_temperature) .or.
associated(fluxes%lprec) .or. &
1047 associated(fluxes%evap))
then
1052 cs%fresh_water_input = cs%fresh_water_input + fw_input
1053 cs%net_salt_input = cs%net_salt_input + salt_input
1054 cs%net_heat_input = cs%net_heat_input + heat_input
1056 cs%fresh_water_in_EFP = cs%fresh_water_in_EFP + fw_in_efp
1057 cs%net_salt_in_EFP = cs%net_salt_in_EFP + salt_in_efp
1058 cs%net_heat_in_EFP = cs%net_heat_in_EFP + heat_in_efp
1061 end subroutine accumulate_net_input
1067 subroutine depth_list_setup(G, US, CS)
1075 if (cs%read_depth_list)
then
1077 call read_depth_list(g, us, cs, cs%depth_list_file)
1079 if (is_root_pe())
call mom_error(warning,
"depth_list_setup: "// &
1080 trim(cs%depth_list_file)//
" does not exist. Creating a new file.")
1081 call create_depth_list(g, cs)
1083 call write_depth_list(g, us, cs, cs%depth_list_file, cs%list_size+1)
1086 call create_depth_list(g, cs)
1090 cs%lH(k) = cs%list_size
1093 end subroutine depth_list_setup
1097 subroutine create_depth_list(G, CS)
1102 real,
dimension(G%Domain%niglobal*G%Domain%njglobal + 1) :: &
1103 Dlist, & !< The global list of bottom depths [Z ~> m].
1105 integer,
dimension(G%Domain%niglobal*G%Domain%njglobal+1) :: &
1112 logical :: add_to_list
1114 integer :: ir, indxt
1115 integer :: mls, list_size
1116 integer :: list_pos, i_global, j_global
1117 integer :: i, j, k, kl
1119 mls = g%Domain%niglobal*g%Domain%njglobal
1124 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1126 j_global = j + g%jdg_offset - (g%jsg-1)
1127 i_global = i + g%idg_offset - (g%isg-1)
1129 list_pos = (j_global-1)*g%Domain%niglobal + i_global
1130 dlist(list_pos) = g%bathyT(i,j)
1131 arealist(list_pos) = g%mask2dT(i,j) * g%areaT(i,j)
1135 call sum_across_pes(dlist, mls+1)
1136 call sum_across_pes(arealist, mls+1)
1138 do j=1,mls+1 ; indx2(j) = j ;
enddo
1139 k = mls / 2 + 1 ; ir = mls
1148 indx2(ir) = indx2(1)
1150 if (ir == 1)
then ; indx2(1) = indxt ;
exit ;
endif
1153 do ;
if (j > ir)
exit
1154 if (j < ir .AND. dlist(indx2(j)) < dlist(indx2(j+1))) j = j + 1
1155 if (dnow < dlist(indx2(j)))
then ; indx2(i) = indx2(j) ; i = j ; j = j + i
1156 else ; j = ir+1 ;
endif
1168 d_list_prev = dlist(indx2(mls))
1171 if (dlist(indx2(k)) < d_list_prev-cs%D_list_min_inc)
then
1172 list_size = list_size + 1
1173 d_list_prev = dlist(indx2(k))
1177 cs%list_size = list_size
1178 allocate(cs%DL(cs%list_size+1))
1180 vol = 0.0 ; area = 0.0
1181 dprev = dlist(indx2(mls))
1187 vol = vol + area * (dprev - dlist(i))
1188 area = area + arealist(i)
1190 add_to_list = .false.
1191 if ((kl == 0) .or. (k==1))
then
1192 add_to_list = .true.
1193 elseif (dlist(indx2(k-1)) < d_list_prev-cs%D_list_min_inc)
then
1194 add_to_list = .true.
1195 d_list_prev = dlist(indx2(k-1))
1198 if (add_to_list)
then
1200 cs%DL(kl)%depth = dlist(i)
1201 cs%DL(kl)%area = area
1202 cs%DL(kl)%vol_below = vol
1207 do while (kl < list_size)
1210 cs%DL(kl)%vol_below = cs%DL(kl-1)%vol_below * 1.000001
1211 cs%DL(kl)%area = cs%DL(kl-1)%area
1212 cs%DL(kl)%depth = cs%DL(kl-1)%depth
1215 cs%DL(cs%list_size+1)%vol_below = cs%DL(cs%list_size)%vol_below * 1000.0
1216 cs%DL(cs%list_size+1)%area = cs%DL(cs%list_size)%area
1217 cs%DL(cs%list_size+1)%depth = cs%DL(cs%list_size)%depth
1219 end subroutine create_depth_list
1222 subroutine write_depth_list(G, US, CS, filename, list_size)
1227 character(len=*),
intent(in) :: filename
1228 integer,
intent(in) :: list_size
1230 real,
allocatable :: tmp(:)
1231 integer :: ncid, dimid(1), Did, Aid, Vid, status, k
1232 character(len=16) :: depth_chksum, area_chksum
1235 call get_depth_list_checksums(g, depth_chksum, area_chksum)
1237 if (.not.is_root_pe())
return
1239 allocate(tmp(list_size)) ; tmp(:) = 0.0
1241 status = nf90_create(filename, 0, ncid)
1242 if (status /= nf90_noerr)
then
1243 call mom_error(warning, filename//trim(nf90_strerror(status)))
1247 status = nf90_def_dim(ncid,
"list", list_size, dimid(1))
1248 if (status /= nf90_noerr)
call mom_error(warning, &
1249 filename//trim(nf90_strerror(status)))
1251 status = nf90_def_var(ncid,
"depth", nf90_double, dimid, did)
1252 if (status /= nf90_noerr)
call mom_error(warning, &
1253 filename//
" depth "//trim(nf90_strerror(status)))
1254 status = nf90_put_att(ncid, did,
"long_name",
"Sorted depth")
1255 if (status /= nf90_noerr)
call mom_error(warning, &
1256 filename//
" depth "//trim(nf90_strerror(status)))
1257 status = nf90_put_att(ncid, did,
"units",
"m")
1258 if (status /= nf90_noerr)
call mom_error(warning, &
1259 filename//
" depth "//trim(nf90_strerror(status)))
1261 status = nf90_def_var(ncid,
"area", nf90_double, dimid, aid)
1262 if (status /= nf90_noerr)
call mom_error(warning, &
1263 filename//
" area "//trim(nf90_strerror(status)))
1264 status = nf90_put_att(ncid, aid,
"long_name",
"Open area at depth")
1265 if (status /= nf90_noerr)
call mom_error(warning, &
1266 filename//
" area "//trim(nf90_strerror(status)))
1267 status = nf90_put_att(ncid, aid,
"units",
"m2")
1268 if (status /= nf90_noerr)
call mom_error(warning, &
1269 filename//
" area "//trim(nf90_strerror(status)))
1271 status = nf90_def_var(ncid,
"vol_below", nf90_double, dimid, vid)
1272 if (status /= nf90_noerr)
call mom_error(warning, &
1273 filename//
" vol_below "//trim(nf90_strerror(status)))
1274 status = nf90_put_att(ncid, vid,
"long_name",
"Open volume below depth")
1275 if (status /= nf90_noerr)
call mom_error(warning, &
1276 filename//
" vol_below "//trim(nf90_strerror(status)))
1277 status = nf90_put_att(ncid, vid,
"units",
"m3")
1278 if (status /= nf90_noerr)
call mom_error(warning, &
1279 filename//
" vol_below "//trim(nf90_strerror(status)))
1282 status = nf90_put_att(ncid, nf90_global, depth_chksum_attr, depth_chksum)
1283 if (status /= nf90_noerr)
call mom_error(warning, &
1284 filename//
" "//depth_chksum_attr//
" "//trim(nf90_strerror(status)))
1286 status = nf90_put_att(ncid, nf90_global, area_chksum_attr, area_chksum)
1287 if (status /= nf90_noerr)
call mom_error(warning, &
1288 filename//
" "//area_chksum_attr//
" "//trim(nf90_strerror(status)))
1290 status = nf90_enddef(ncid)
1291 if (status /= nf90_noerr)
call mom_error(warning, &
1292 filename//trim(nf90_strerror(status)))
1294 do k=1,list_size ; tmp(k) = us%Z_to_m*cs%DL(k)%depth ;
enddo
1295 status = nf90_put_var(ncid, did, tmp)
1296 if (status /= nf90_noerr)
call mom_error(warning, &
1297 filename//
" depth "//trim(nf90_strerror(status)))
1299 do k=1,list_size ; tmp(k) = cs%DL(k)%area ;
enddo
1300 status = nf90_put_var(ncid, aid, tmp)
1301 if (status /= nf90_noerr)
call mom_error(warning, &
1302 filename//
" area "//trim(nf90_strerror(status)))
1304 do k=1,list_size ; tmp(k) = us%Z_to_m*cs%DL(k)%vol_below ;
enddo
1305 status = nf90_put_var(ncid, vid, tmp)
1306 if (status /= nf90_noerr)
call mom_error(warning, &
1307 filename//
" vol_below "//trim(nf90_strerror(status)))
1309 status = nf90_close(ncid)
1310 if (status /= nf90_noerr)
call mom_error(warning, &
1311 filename//trim(nf90_strerror(status)))
1313 end subroutine write_depth_list
1317 subroutine read_depth_list(G, US, CS, filename)
1322 character(len=*),
intent(in) :: filename
1324 character(len=32) :: mdl
1325 character(len=240) :: var_name, var_msg
1326 real,
allocatable :: tmp(:)
1327 integer :: ncid, status, varid, list_size, k
1328 integer :: ndim, len, var_dim_ids(NF90_MAX_VAR_DIMS)
1329 character(len=16) :: depth_file_chksum, depth_grid_chksum
1330 character(len=16) :: area_file_chksum, area_grid_chksum
1331 integer :: depth_attr_status, area_attr_status
1333 mdl =
"MOM_sum_output read_depth_list:"
1335 status = nf90_open(filename, nf90_nowrite, ncid)
1336 if (status /= nf90_noerr)
then
1337 call mom_error(fatal,mdl//
" Difficulties opening "//trim(filename)// &
1338 " - "//trim(nf90_strerror(status)))
1342 depth_attr_status = nf90_get_att(ncid, nf90_global, depth_chksum_attr, &
1344 area_attr_status = nf90_get_att(ncid, nf90_global, area_chksum_attr, &
1347 if (any([depth_attr_status, area_attr_status] == nf90_enotatt))
then
1348 var_msg = trim(cs%depth_list_file) //
" checksums are missing;"
1349 if (cs%require_depth_list_chksum)
then
1350 call mom_error(fatal, trim(var_msg) //
" aborting.")
1351 elseif (cs%update_depth_list_chksum)
then
1352 call mom_error(warning, trim(var_msg) //
" updating file.")
1353 call create_depth_list(g, cs)
1354 call write_depth_list(g, us, cs, cs%depth_list_file, cs%list_size+1)
1357 call mom_error(warning, &
1358 trim(var_msg) //
" some diagnostics may not be reproducible.")
1362 if (depth_attr_status /= nf90_noerr)
then
1363 var_msg = mdl //
"Failed to read " // trim(filename) //
":" &
1364 // depth_chksum_attr
1365 call mom_error(fatal, &
1366 trim(var_msg) //
" - " // nf90_strerror(depth_attr_status))
1369 if (area_attr_status /= nf90_noerr)
then
1370 var_msg = mdl //
"Failed to read " // trim(filename) //
":" &
1372 call mom_error(fatal, &
1373 trim(var_msg) //
" - " // nf90_strerror(area_attr_status))
1376 call get_depth_list_checksums(g, depth_grid_chksum, area_grid_chksum)
1378 if (depth_grid_chksum /= depth_file_chksum &
1379 .or. area_grid_chksum /= area_file_chksum)
then
1380 var_msg = trim(cs%depth_list_file) //
" checksums do not match;"
1381 if (cs%require_depth_list_chksum)
then
1382 call mom_error(fatal, trim(var_msg) //
" aborting.")
1383 elseif (cs%update_depth_list_chksum)
then
1384 call mom_error(warning, trim(var_msg) //
" updating file.")
1385 call create_depth_list(g, cs)
1386 call write_depth_list(g, us, cs, cs%depth_list_file, cs%list_size+1)
1389 call mom_error(warning, &
1390 trim(var_msg) //
" some diagnostics may not be reproducible.")
1396 var_msg = trim(var_name)//
" in "//trim(filename)//
" - "
1397 status = nf90_inq_varid(ncid, var_name, varid)
1398 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1399 " Difficulties finding variable "//trim(var_msg)//&
1400 trim(nf90_strerror(status)))
1402 status = nf90_inquire_variable(ncid, varid, ndims=ndim, dimids=var_dim_ids)
1403 if (status /= nf90_noerr)
then
1404 call mom_error(fatal,mdl//
" cannot inquire about "//trim(var_msg)//&
1405 trim(nf90_strerror(status)))
1406 elseif (ndim > 1)
then
1407 call mom_error(fatal,mdl//
" "//trim(var_msg)//&
1408 " has too many or too few dimensions.")
1412 status = nf90_inquire_dimension(ncid, var_dim_ids(1), len=list_size)
1413 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1414 " cannot inquire about dimension(1) of "//trim(var_msg)//&
1415 trim(nf90_strerror(status)))
1417 cs%list_size = list_size-1
1418 allocate(cs%DL(list_size))
1419 allocate(tmp(list_size))
1421 status = nf90_get_var(ncid, varid, tmp)
1422 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1423 " Difficulties reading variable "//trim(var_msg)//&
1424 trim(nf90_strerror(status)))
1426 do k=1,list_size ; cs%DL(k)%depth = us%m_to_Z*tmp(k) ;
enddo
1429 var_msg = trim(var_name)//
" in "//trim(filename)//
" - "
1430 status = nf90_inq_varid(ncid, var_name, varid)
1431 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1432 " Difficulties finding variable "//trim(var_msg)//&
1433 trim(nf90_strerror(status)))
1434 status = nf90_get_var(ncid, varid, tmp)
1435 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1436 " Difficulties reading variable "//trim(var_msg)//&
1437 trim(nf90_strerror(status)))
1439 do k=1,list_size ; cs%DL(k)%area = tmp(k) ;
enddo
1441 var_name =
"vol_below"
1442 var_msg = trim(var_name)//
" in "//trim(filename)
1443 status = nf90_inq_varid(ncid, var_name, varid)
1444 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1445 " Difficulties finding variable "//trim(var_msg)//&
1446 trim(nf90_strerror(status)))
1447 status = nf90_get_var(ncid, varid, tmp)
1448 if (status /= nf90_noerr)
call mom_error(fatal,mdl// &
1449 " Difficulties reading variable "//trim(var_msg)//&
1450 trim(nf90_strerror(status)))
1452 do k=1,list_size ; cs%DL(k)%vol_below = us%m_to_Z*tmp(k) ;
enddo
1454 status = nf90_close(ncid)
1455 if (status /= nf90_noerr)
call mom_error(warning, mdl// &
1456 " Difficulties closing "//trim(filename)//
" - "//trim(nf90_strerror(status)))
1460 end subroutine read_depth_list
1473 subroutine get_depth_list_checksums(G, depth_chksum, area_chksum)
1475 character(len=16),
intent(out) :: depth_chksum
1476 character(len=16),
intent(out) :: area_chksum
1479 real,
allocatable :: field(:,:)
1481 allocate(field(g%isc:g%iec, g%jsc:g%jec))
1484 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1485 field(i,j) = g%bathyT(i,j)
1487 write(depth_chksum,
'(Z16)') mpp_chksum(field(:,:))
1490 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1491 field(i,j) = g%mask2dT(i,j) * g%areaT(i,j)
1493 write(area_chksum,
'(Z16)') mpp_chksum(field(:,:))
1496 end subroutine get_depth_list_checksums