tidal_bay_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Configures the model for the "tidal_bay" experiment.
6!! tidal_bay = Tidally resonant bay from Zygmunt Kowalik's class on tides.
8
9use mom_coms, only : reproducing_sum
11use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
12use mom_file_parser, only : get_param, log_version, param_file_type
13use mom_grid, only : ocean_grid_type
19use mom_time_manager, only : time_type, time_to_real
20
21implicit none ; private
22
23#include <MOM_memory.h>
24
27
28!> Control structure for tidal bay open boundaries.
29type, public :: tidal_bay_obc_cs ; private
30 real :: tide_flow = 3.0e6 !< Maximum tidal flux with the tidal bay configuration [L2 Z T-1 ~> m3 s-1]
31 real :: tide_period !< The period associated with the tidal bay configuration [T ~> s]
32 real :: tide_ssh_amp !< The magnitude of the sea surface height anomalies at the inflow
33 !! with the tidal bay configuration [Z ~> m]
34end type tidal_bay_obc_cs
35
36contains
37
38!> Add tidal bay to OBC registry.
39function register_tidal_bay_obc(param_file, CS, US, OBC_Reg)
40 type(param_file_type), intent(in) :: param_file !< parameter file.
41 type(tidal_bay_obc_cs), intent(inout) :: cs !< tidal bay control structure.
42 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
43 type(obc_registry_type), pointer :: obc_reg !< OBC registry.
44 logical :: register_tidal_bay_obc
45 character(len=32) :: casename = "tidal bay" !< This case's name.
46 character(len=40) :: mdl = "tidal_bay_initialization" ! This module's name.
47
48 call get_param(param_file, mdl, "TIDAL_BAY_FLOW", cs%tide_flow, &
49 "Maximum total tidal volume flux.", &
50 units="m3 s-1", default=3.0e6, scale=us%m_s_to_L_T*us%m_to_L*us%m_to_Z)
51 call get_param(param_file, mdl, "TIDAL_BAY_PERIOD", cs%tide_period, &
52 "Period of the inflow in the tidal bay configuration.", &
53 units="s", default=12.0*3600.0, scale=us%s_to_T)
54 call get_param(param_file, mdl, "TIDAL_BAY_SSH_ANOM", cs%tide_ssh_amp, &
55 "Magnitude of the sea surface height anomalies at the inflow with the "//&
56 "tidal bay configuration.", &
57 units="m", default=0.1, scale=us%m_to_Z)
58
59 ! Register the open boundaries.
60 call register_obc(casename, param_file, obc_reg)
62
63end function register_tidal_bay_obc
64
65!> This subroutine sets the properties of flow at open boundary conditions.
66subroutine tidal_bay_set_obc_data(OBC, CS, G, GV, US, h, Time)
67 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
68 !! whether, where, and what open boundary
69 !! conditions are used.
70 type(tidal_bay_obc_cs), intent(in) :: cs !< tidal bay control structure.
71 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
72 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
73 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
74 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
75 type(time_type), intent(in) :: time !< model time.
76
77 ! The following variables are used to set up the transport in the tidal_bay example.
78 real :: time_sec ! Elapsed model time [T ~> s]
79 real :: cff_eta ! The sea surface height anomalies associated with the inflow [Z ~> m]
80 real :: my_flux ! The volume flux through the face [L2 Z T-1 ~> m3 s-1]
81 real :: total_area ! The total face area of the OBCs [L Z ~> m2]
82 real :: normal_vel ! The normal velocity through the inflow face [L T-1 ~> m s-1]
83 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
84 real, allocatable :: my_area(:,:) ! The total OBC inflow area [L Z ~> m2]
85 integer :: turns ! Number of index quarter turns
86 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, n
87 integer :: isdb, iedb, jsdb, jedb
88 type(obc_segment_type), pointer :: segment => null()
89
90 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
91 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
92 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
93
94 pi = 4.0*atan(1.0)
95
96 turns = modulo(g%HI%turns, 4)
97
98 if (.not.associated(obc)) return
99
100 time_sec = time_to_real(time, scale=us%s_to_T)
101 cff_eta = cs%tide_ssh_amp * sin(2.0*pi*time_sec / cs%tide_period)
102
103 segment => obc%segment(1)
104
105 if (turns == 0) then
106 allocate(my_area(1:1,js:je), source=0.0)
107 do j=segment%HI%jsc,segment%HI%jec ; do i=segment%HI%IscB,segment%HI%IecB
108 if (obc%segnum_u(i,j) > 0) then ! (segment%direction == OBC_DIRECTION_E)
109 do k=1,nz
110 my_area(1,j) = my_area(1,j) + h(i,j,k)*(gv%H_to_m*us%m_to_Z)*g%dyCu(i,j)
111 enddo
112 endif
113 enddo ; enddo
114 elseif (turns == 1) then
115 allocate(my_area(is:ie,1:1), source=0.0)
116 do j=segment%HI%JscB,segment%HI%JecB ; do i=segment%HI%isc,segment%HI%iec
117 if (obc%segnum_v(i,j) > 0) then ! (segment%direction == OBC_DIRECTION_N)
118 do k=1,nz
119 my_area(i,1) = my_area(i,1) + h(i,j,k)*(gv%H_to_m*us%m_to_Z)*g%dxCv(i,j)
120 enddo
121 endif
122 enddo ; enddo
123 elseif (turns == 2) then
124 allocate(my_area(1:1,js:je), source=0.0)
125 do j=segment%HI%jsc,segment%HI%jec ; do i=segment%HI%IscB,segment%HI%IecB
126 if (obc%segnum_u(i,j) < 0) then ! (segment%direction == OBC_DIRECTION_W)
127 do k=1,nz
128 my_area(1,j) = my_area(1,j) + h(i+1,j,k)*(gv%H_to_m*us%m_to_Z)*g%dyCu(i,j)
129 enddo
130 endif
131 enddo ; enddo
132 elseif (turns == 3) then
133 allocate(my_area(is:ie,1:1), source=0.0)
134 do j=segment%HI%JscB,segment%HI%JecB ; do i=segment%HI%isc,segment%HI%iec
135 if (obc%segnum_v(i,j) < 0) then ! (segment%direction == OBC_DIRECTION_S)
136 do k=1,nz
137 my_area(i,1) = my_area(i,1) + h(i,j+1,k)*(gv%H_to_m*us%m_to_Z)*g%dxCv(i,j)
138 enddo
139 endif
140 enddo ; enddo
141 endif
142
143 total_area = reproducing_sum(my_area, unscale=us%Z_to_m*us%L_to_m)
144 my_flux = - cs%tide_flow * sin(2.0*pi*time_sec / cs%tide_period)
145 normal_vel = my_flux / total_area
146 if ((turns==2) .or. (turns==3)) normal_vel = -1.0 * normal_vel
147
148 do n = 1, obc%number_of_segments
149 segment => obc%segment(n)
150 if (.not. segment%on_pe) cycle
151
152 segment%normal_vel_bt(:,:) = normal_vel
153 segment%SSH(:,:) = cff_eta
154
155 enddo ! end segment loop
156
157end subroutine tidal_bay_set_obc_data
158