Skip to content

Commit aab5f48

Browse files
Raphael DussinRaphael Dussin
Raphael Dussin
authored and
Raphael Dussin
committed
add restart rescaling by 2**power
1 parent ceb92ee commit aab5f48

File tree

1 file changed

+41
-39
lines changed

1 file changed

+41
-39
lines changed

‎src/parameterizations/lateral/MOM_internal_tides.F90‎

Lines changed: 41 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ module MOM_internal_tides
2020
use MOM_grid, only : ocean_grid_type
2121
use MOM_int_tide_input, only: int_tide_input_CS, get_input_TKE, get_barotropic_tidal_vel
2222
use MOM_io, only : slasher, MOM_read_data, file_exists, axis_info
23-
use MOM_io, only : set_axis_info, get_axis_info
23+
use MOM_io, only : set_axis_info, get_axis_info, stdout
2424
use MOM_restart, only : register_restart_field, MOM_restart_CS, restart_init, save_restart
2525
use MOM_restart, only : lock_check, restart_registry_lock
2626
use MOM_spatial_means, only : global_area_integral
@@ -148,6 +148,7 @@ module MOM_internal_tides
148148
real :: q_itides !< fraction of local dissipation [nondim]
149149
real :: En_sum !< global sum of energy for use in debugging, in MKS units [H Z2 T-2 L2 ~> m5 s-2 or J]
150150
real :: En_underflow !< A minuscule amount of energy [H Z2 T-2 ~> m3 s-2 or J m-2]
151+
integer :: En_restart_power !< A power factor of 2 by which to multiply the energy in restart [nondim]
151152
type(time_type), pointer :: Time => NULL() !< A pointer to the model's clock.
152153
character(len=200) :: inputdir !< directory to look for coastline angle file
153154
real :: decay_rate !< A constant rate at which internal tide energy is
@@ -308,6 +309,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
308309
real :: I_D_here ! The inverse of the local water column thickness [H-1 ~> m-1 or m2 kg-1]
309310
real :: I_mass ! The inverse of the local water mass [R-1 Z-1 ~> m2 kg-1]
310311
real :: I_dt ! The inverse of the timestep [T-1 ~> s-1]
312+
real :: En_restart_factor ! A multiplicative factor of the form 2**En_restart_power [nondim]
313+
real :: I_En_restart_factor ! The inverse of the restart mult factor [nondim]
311314
real :: freq2 ! The frequency squared [T-2 ~> s-2]
312315
real :: PE_term ! total potential energy of profile [R Z ~> kg m-2]
313316
real :: KE_term ! total kinetic energy of profile [R Z ~> kg m-2]
@@ -347,6 +350,8 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
347350
en_subRO = 1e-30*J_m2_to_HZ2_T2
348351

349352
I_dt = 1.0 / dt
353+
En_restart_factor = 2**CS%En_restart_power
354+
I_En_restart_factor = 1.0 / En_restart_factor
350355

351356
! initialize local arrays
352357
TKE_itidal_input(:,:,:) = 0.
@@ -360,30 +365,30 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
360365

361366
! Rebuild energy density array from multiple restarts
362367
do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
363-
CS%En(i,j,a,fr,1) = CS%En_restart_mode1(i,j,a,fr)
368+
CS%En(i,j,a,fr,1) = CS%En_restart_mode1(i,j,a,fr) * I_En_restart_factor
364369
enddo ; enddo ; enddo ; enddo
365370

366371
if (CS%nMode >= 2) then
367372
do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
368-
CS%En(i,j,a,fr,2) = CS%En_restart_mode2(i,j,a,fr)
373+
CS%En(i,j,a,fr,2) = CS%En_restart_mode2(i,j,a,fr) * I_En_restart_factor
369374
enddo ; enddo ; enddo ; enddo
370375
endif
371376

372377
if (CS%nMode >= 3) then
373378
do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
374-
CS%En(i,j,a,fr,3) = CS%En_restart_mode3(i,j,a,fr)
379+
CS%En(i,j,a,fr,3) = CS%En_restart_mode3(i,j,a,fr) * I_En_restart_factor
375380
enddo ; enddo ; enddo ; enddo
376381
endif
377382

378383
if (CS%nMode >= 4) then
379384
do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
380-
CS%En(i,j,a,fr,4) = CS%En_restart_mode4(i,j,a,fr)
385+
CS%En(i,j,a,fr,4) = CS%En_restart_mode4(i,j,a,fr) * I_En_restart_factor
381386
enddo ; enddo ; enddo ; enddo
382387
endif
383388

384389
if (CS%nMode >= 5) then
385390
do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
386-
CS%En(i,j,a,fr,5) = CS%En_restart_mode5(i,j,a,fr)
391+
CS%En(i,j,a,fr,5) = CS%En_restart_mode5(i,j,a,fr) * I_En_restart_factor
387392
enddo ; enddo ; enddo ; enddo
388393
endif
389394

@@ -501,7 +506,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
501506
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
502507
do m=1,CS%nMode ; do fr=1,CS%Nfreq
503508
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after forcing')
504-
if (is_root_pe()) print *, 'prop_int_tide: after forcing', CS%En_sum
509+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after forcing', CS%En_sum
505510
enddo ; enddo
506511
endif
507512

@@ -528,7 +533,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
528533
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
529534
do m=1,CS%nMode ; do fr=1,CS%Nfreq
530535
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after 1/2 refraction')
531-
if (is_root_pe()) print *, 'prop_int_tide: after 1/2 refraction', CS%En_sum
536+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after 1/2 refraction', CS%En_sum
532537
enddo ; enddo
533538
! Check for En<0 - for debugging, delete later
534539
do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle
@@ -558,7 +563,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
558563
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo R", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
559564
do m=1,CS%nMode ; do fr=1,CS%Nfreq
560565
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after correct halo rotation')
561-
if (is_root_pe()) print *, 'prop_int_tide: after correct halo rotation', CS%En_sum
566+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after correct halo rotation', CS%En_sum
562567
enddo ; enddo
563568
endif
564569

@@ -589,7 +594,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
589594
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af prop", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
590595
do m=1,CS%nMode ; do fr=1,CS%Nfreq
591596
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after propagate')
592-
if (is_root_pe()) print *, 'prop_int_tide: after propagate', CS%En_sum
597+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after propagate', CS%En_sum
593598
enddo ; enddo
594599
! Check for En<0 - for debugging, delete later
595600
do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle
@@ -631,7 +636,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
631636
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr2", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
632637
do m=1,CS%nMode ; do fr=1,CS%Nfreq
633638
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after 2/2 refraction')
634-
if (is_root_pe()) print *, 'prop_int_tide: after 2/2 refraction', CS%En_sum
639+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after 2/2 refraction', CS%En_sum
635640
enddo ; enddo
636641
! Check for En<0 - for debugging, delete later
637642
do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle
@@ -687,9 +692,9 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
687692
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after leak", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
688693
do m=1,CS%nMode ; do fr=1,CS%Nfreq
689694
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after background drag')
690-
if (is_root_pe()) print *, 'prop_int_tide: after background drag', CS%En_sum
695+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after background drag', CS%En_sum
691696
call sum_En(G, GV, US, CS, CS%TKE_leak_loss(:,:,:,fr,m) * dt, 'prop_int_tide: loss after background drag')
692-
if (is_root_pe()) print *, 'prop_int_tide: loss after background drag', CS%En_sum
697+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: loss after background drag', CS%En_sum
693698
enddo ; enddo
694699
! Check for En<0 - for debugging, delete later
695700
do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle
@@ -858,7 +863,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
858863
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after wave", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
859864
do m=1,CS%nMode ; do fr=1,CS%Nfreq
860865
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: before Froude drag')
861-
if (is_root_pe()) print *, 'prop_int_tide: before Froude drag', CS%En_sum
866+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: before Froude drag', CS%En_sum
862867
enddo ; enddo
863868
! save loss term for online budget, may want to add a debug flag later
864869
do m=1,CS%nMode ; do fr=1,CS%nFreq
@@ -932,9 +937,9 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
932937
call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after froude", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2)
933938
do m=1,CS%nMode ; do fr=1,CS%Nfreq
934939
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after Froude drag')
935-
if (is_root_pe()) print *, 'prop_int_tide: after Froude drag', CS%En_sum
940+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after Froude drag', CS%En_sum
936941
call sum_En(G, GV, US, CS, CS%TKE_Froude_loss(:,:,:,fr,m) * dt, 'prop_int_tide: loss after Froude drag')
937-
if (is_root_pe()) print *, 'prop_int_tide: loss after Froude drag', CS%En_sum
942+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: loss after Froude drag', CS%En_sum
938943
enddo ; enddo
939944
! save loss term for online budget, may want to add a debug flag later
940945
do m=1,CS%nMode ; do fr=1,CS%nFreq
@@ -1015,7 +1020,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
10151020
CS%TKE_quad_loss_glo_dt(fr,m) - CS%TKE_itidal_loss_glo_dt(fr,m) - &
10161021
CS%TKE_Froude_loss_glo_dt(fr,m) - CS%TKE_residual_loss_glo_dt(fr,m) - &
10171022
CS%En_end_glo(fr,m)
1018-
if (is_root_pe()) print *, "error in Energy budget", CS%error_mode(fr,m)
1023+
if (is_root_pe()) write(stdout,'(A,F18.10)'), "error in Energy budget", CS%error_mode(fr,m)
10191024
enddo ; enddo
10201025
endif
10211026

@@ -1048,37 +1053,32 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
10481053
call post_data(CS%id_En_mode(fr,m), tot_En, CS%diag)
10491054
endif ; enddo ; enddo
10501055

1051-
! fix underflows
1052-
do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
1053-
if (abs(CS%En(i,j,a,fr,m)) < CS%En_underflow) CS%En(i,j,a,fr,m) = 0.0
1054-
enddo ; enddo ; enddo ; enddo ; enddo
1055-
10561056
! split energy array into multiple restarts
10571057
do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
1058-
CS%En_restart_mode1(i,j,a,fr) = CS%En(i,j,a,fr,1)
1058+
CS%En_restart_mode1(i,j,a,fr) = CS%En(i,j,a,fr,1) * En_restart_factor
10591059
enddo ; enddo ; enddo ; enddo
10601060

10611061
if (CS%nMode >= 2) then
10621062
do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
1063-
CS%En_restart_mode2(i,j,a,fr) = CS%En(i,j,a,fr,2)
1063+
CS%En_restart_mode2(i,j,a,fr) = CS%En(i,j,a,fr,2) * En_restart_factor
10641064
enddo ; enddo ; enddo ; enddo
10651065
endif
10661066

10671067
if (CS%nMode >= 3) then
10681068
do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
1069-
CS%En_restart_mode3(i,j,a,fr) = CS%En(i,j,a,fr,3)
1069+
CS%En_restart_mode3(i,j,a,fr) = CS%En(i,j,a,fr,3) * En_restart_factor
10701070
enddo ; enddo ; enddo ; enddo
10711071
endif
10721072

10731073
if (CS%nMode >= 4) then
10741074
do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
1075-
CS%En_restart_mode4(i,j,a,fr) = CS%En(i,j,a,fr,4)
1075+
CS%En_restart_mode4(i,j,a,fr) = CS%En(i,j,a,fr,4) * En_restart_factor
10761076
enddo ; enddo ; enddo ; enddo
10771077
endif
10781078

10791079
if (CS%nMode >= 5) then
10801080
do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
1081-
CS%En_restart_mode5(i,j,a,fr) = CS%En(i,j,a,fr,5)
1081+
CS%En_restart_mode5(i,j,a,fr) = CS%En(i,j,a,fr,5) * En_restart_factor
10821082
enddo ; enddo ; enddo ; enddo
10831083
endif
10841084

@@ -1608,23 +1608,23 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
16081608
enddo
16091609

16101610
if (abs(verif_N -1.0) > threshold_verif) then
1611-
print *, i, j, verif_N
1611+
write(stdout,'(I5,I5,F18.10)'), i, j, verif_N
16121612
call MOM_error(FATAL, "mismatch integral for N profile")
16131613
endif
16141614
if (abs(verif_N2 -1.0) > threshold_verif) then
1615-
print *, i, j, verif_N2
1615+
write(stdout,'(I5,I5,F18.10)'), i, j, verif_N2
16161616
call MOM_error(FATAL, "mismatch integral for N2 profile")
16171617
endif
16181618
if (abs(verif_bbl -1.0) > threshold_verif) then
1619-
print *, i, j, verif_bbl
1619+
write(stdout,'(I5,I5,F18.10)'), i, j, verif_bbl
16201620
call MOM_error(FATAL, "mismatch integral for bbl profile")
16211621
endif
16221622
if (abs(verif_stl1 -1.0) > threshold_verif) then
1623-
print *, i, j, verif_stl1
1623+
write(stdout,'(I5,I5,F18.10)'), i, j, verif_stl1
16241624
call MOM_error(FATAL, "mismatch integral for stl1 profile")
16251625
endif
16261626
if (abs(verif_stl2 -1.0) > threshold_verif) then
1627-
print *, i, j, verif_stl2
1627+
write(stdout,'(I5,I5,F18.10)'), i, j, verif_stl2
16281628
call MOM_error(FATAL, "mismatch integral for stl2 profile")
16291629
endif
16301630

@@ -2104,7 +2104,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
21042104
if (CS%debug) then
21052105
do m=1,CS%nMode ; do fr=1,CS%Nfreq
21062106
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: top of routine')
2107-
if (is_root_pe()) print *, 'propagate: top of routine', CS%En_sum
2107+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: top of routine', CS%En_sum
21082108
enddo ; enddo
21092109
endif
21102110

@@ -2176,7 +2176,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
21762176
if (CS%debug) then
21772177
do m=1,CS%nMode ; do fr=1,CS%Nfreq
21782178
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: after propagate_x')
2179-
if (is_root_pe()) print *, 'propagate: after propagate_x', CS%En_sum
2179+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: after propagate_x', CS%En_sum
21802180
enddo ; enddo
21812181
endif
21822182

@@ -2187,7 +2187,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
21872187
if (CS%debug) then
21882188
do m=1,CS%nMode ; do fr=1,CS%Nfreq
21892189
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: after halo update')
2190-
if (is_root_pe()) print *, 'propagate: after halo update', CS%En_sum
2190+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: after halo update', CS%En_sum
21912191
enddo ; enddo
21922192
endif
21932193
! Apply propagation in y-direction (reflection included)
@@ -2206,7 +2206,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
22062206
if (CS%debug) then
22072207
do m=1,CS%nMode ; do fr=1,CS%Nfreq
22082208
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: after propagate_y')
2209-
if (is_root_pe()) print *, 'propagate: after propagate_y', CS%En_sum
2209+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: after propagate_y', CS%En_sum
22102210
enddo ; enddo
22112211
endif
22122212

@@ -2215,7 +2215,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
22152215
if (CS%debug) then
22162216
do m=1,CS%nMode ; do fr=1,CS%Nfreq
22172217
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: bottom of routine')
2218-
if (is_root_pe()) print *, 'propagate: bottom of routine', CS%En_sum
2218+
if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: bottom of routine', CS%En_sum
22192219
enddo ; enddo
22202220
endif
22212221

@@ -3573,8 +3573,10 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
35733573
do_not_log=.not.CS%apply_Froude_drag)
35743574
call get_param(param_file, mdl, "EN_UNDERFLOW", CS%En_underflow, &
35753575
"A small energy density below which Energy is set to zero.", &
3576-
units="J m-2", default=1.0e-100, scale=J_m2_to_HZ2_T2, &
3577-
do_not_log=.not.CS%apply_Froude_drag)
3576+
units="J m-2", default=1.0e-100, scale=J_m2_to_HZ2_T2)
3577+
call get_param(param_file, mdl, "EN_RESTART_POWER", CS%En_restart_power, &
3578+
"A power factor to save larger values x 2**(power) in restart files.", &
3579+
units="nondim", default=0)
35783580
call get_param(param_file, mdl, "CDRAG", CS%cdrag, &
35793581
"CDRAG is the drag coefficient relating the magnitude of "//&
35803582
"the velocity field to the bottom stress.", &

0 commit comments

Comments
 (0)