RGC_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 models sponges for the Rotating Gravity Current (RGC) experiment.
7
8!***********************************************************************
9!* By Elizabeth Yankovsky, May 2018 *
10!***********************************************************************
11
12use mom_ale_sponge, only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge
13use mom_ale_sponge, only : set_up_ale_sponge_vel_field
14use mom_domains, only : pass_var
16use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe, warning
17use mom_file_parser, only : get_param, log_version, param_file_type
18use mom_get_input, only : directories
19use mom_grid, only : ocean_grid_type
20use mom_io, only : file_exists, mom_read_data, slasher
26use mom_eos, only : calculate_density, eos_domain
27implicit none ; private
28
29#include <MOM_memory.h>
30
32
33! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
34! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
35! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
36! vary with the Boussinesq approximation, the Boussinesq variant is given first.
37
38contains
39
40!> Sets up the inverse restoration time, and the values towards which the interface heights,
41!! velocities and tracers should be restored within the sponges for the RGC test case.
42subroutine rgc_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, use_ALE, CSp, ACSp)
43 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
44 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
45 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
46 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers
47 !! to any available thermodynamic
48 !! fields, potential temperature and
49 !! salinity or mixed layer density.
50 !! Absent fields have NULL pointers.
51 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
52 target, intent(in) :: u !< Array with the u velocity [L T-1 ~> m s-1]
53 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
54 target, intent(in) :: v !< Array with the v velocity [L T-1 ~> m s-1]
55 real, dimension(SZI_(G),SZJ_(G)), &
56 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
57 type(param_file_type), intent(in) :: pf !< A structure to parse for model parameter values.
58 logical, intent(in) :: use_ale !< If true, indicates model is in ALE mode
59 type(sponge_cs), pointer :: csp !< Layer-mode sponge structure
60 type(ale_sponge_cs), pointer :: acsp !< ALE-mode sponge structure
61
62 ! Local variables
63 real :: t(szi_(g),szj_(g),szk_(gv)) ! A temporary array for temperature [C ~> degC]
64 real :: s(szi_(g),szj_(g),szk_(gv)) ! A temporary array for salinity [S ~> ppt]
65 real :: u1(szib_(g),szj_(g),szk_(gv)) ! A temporary array for u [L T-1 ~> m s-1]
66 real :: v1(szi_(g),szjb_(g),szk_(gv)) ! A temporary array for v [L T-1 ~> m s-1]
67 real :: rho(szi_(g),szj_(g)) ! A temporary array for mixed layer density [R ~> kg m-3].
68 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Sponge layer thicknesses in height units [Z ~> m]
69 real :: idamp(szi_(g),szj_(g)) ! The sponge damping rate at h points [T-1 ~> s-1]
70 real :: tnudg ! Nudging time scale [T ~> s]
71 real :: pres(szi_(g)) ! An array of the reference pressure [R L2 T-2 ~> Pa]
72 real :: eta(szi_(g),szj_(g),szk_(gv)+1) ! A temporary array for eta, positive upward [Z ~> m]
73 logical :: sponge_uv ! Nudge velocities (u and v) towards zero
74 real :: min_depth ! The minimum depth of the ocean [Z ~> m]
75 real :: dummy1 ! The position relative to the sponge width [nondim]
76 real :: min_thickness ! A minimum layer thickness [H ~> m or kg m-2] (unused)
77 real :: lensponge ! The width of the sponge in axis units, [km] or [m]
78 character(len=40) :: filename, state_file
79 character(len=40) :: temp_var, salt_var, eta_var, inputdir, h_var
80
81 character(len=40) :: mdl = "RGC_initialize_sponges" ! This subroutine's name.
82 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
83 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, nz, iscb, iecb, jscb, jecb
84
85 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
86 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
87 iscb = g%iscB ; iecb = g%iecB ; jscb = g%jscB ; jecb = g%jecB
88
89 ! The variable min_thickness is unused, and can probably be eliminated.
90 call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
91 units='m', default=1.e-3, scale=gv%m_to_H)
92
93 call get_param(pf, mdl, "RGC_TNUDG", tnudg, 'Nudging time scale for sponge layers', &
94 units='days', default=0.0, scale=86400.0*us%s_to_T)
95
96 call get_param(pf, mdl, "LENSPONGE", lensponge, &
97 "The length of the sponge layer.", &
98 units=g%x_ax_unit_short, default=10.0)
99
100 call get_param(pf, mdl, "SPONGE_UV", sponge_uv, &
101 "Nudge velocities (u and v) towards zero in the sponge layer.", &
102 default=.false., do_not_log=.true.)
103
104 t(:,:,:) = 0.0 ; s(:,:,:) = 0.0 ; idamp(:,:) = 0.0
105
106 call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
107 "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)
108
109 if (associated(csp)) call mom_error(fatal, &
110 "RGC_initialize_sponges called with an associated control structure.")
111 if (associated(acsp)) call mom_error(fatal, &
112 "RGC_initialize_sponges called with an associated ALE-sponge control structure.")
113
114 ! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0
115 ! wherever there is no sponge, and the subroutines that are called
116 ! will automatically set up the sponges only where Idamp is positive
117 ! and mask2dT is 1.
118
119 do j=js,je ; do i=is,ie
120 if ((depth_tot(i,j) <= min_depth) .or. (g%geoLonT(i,j) <= lensponge)) then
121 idamp(i,j) = 0.0
122 elseif (g%geoLonT(i,j) >= (g%len_lon - lensponge) .AND. g%geoLonT(i,j) <= g%len_lon) then
123 dummy1 = (g%geoLonT(i,j)-(g%len_lon - lensponge))/(lensponge)
124 idamp(i,j) = (1.0/tnudg) * max(0.0,dummy1)
125 else
126 idamp(i,j) = 0.0
127 endif
128 enddo ; enddo
129
130
131 ! 1) Read eta, salt and temp from IC file
132 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
133 inputdir = slasher(inputdir)
134 call get_param(pf, mdl, "RGC_SPONGE_FILE", state_file, &
135 "The name of the file with temps., salts. and interfaces to \n"// &
136 " damp toward.", fail_if_missing=.true.)
137 call get_param(pf, mdl, "SPONGE_PTEMP_VAR", temp_var, &
138 "The name of the potential temperature variable in \n"//&
139 "SPONGE_STATE_FILE.", default="Temp")
140 call get_param(pf, mdl, "SPONGE_SALT_VAR", salt_var, &
141 "The name of the salinity variable in \n"//&
142 "SPONGE_STATE_FILE.", default="Salt")
143 call get_param(pf, mdl, "SPONGE_ETA_VAR", eta_var, &
144 "The name of the interface height variable in \n"//&
145 "SPONGE_STATE_FILE.", default="eta")
146 call get_param(pf, mdl, "SPONGE_H_VAR", h_var, &
147 "The name of the layer thickness variable in \n"//&
148 "SPONGE_STATE_FILE.", default="h")
149
150 !read temp and eta
151 filename = trim(inputdir)//trim(state_file)
152 if (.not.file_exists(filename, g%Domain)) &
153 call mom_error(fatal, " RGC_initialize_sponges: Unable to open "//trim(filename))
154 call mom_read_data(filename, temp_var, t(:,:,:), g%Domain, scale=us%degC_to_C)
155 call mom_read_data(filename, salt_var, s(:,:,:), g%Domain, scale=us%ppt_to_S)
156 if (use_ale) then
157
158 call mom_read_data(filename, h_var, dz(:,:,:), g%Domain, scale=us%m_to_Z)
159 call pass_var(dz, g%domain)
160
161 call initialize_ale_sponge(idamp, g, gv, pf, acsp, dz, nz, data_h_is_z=.true.)
162
163 ! The remaining calls to set_up_sponge_field can be in any order.
164 if ( associated(tv%T) ) call set_up_ale_sponge_field(t, g, gv, tv%T, acsp, 'temp', &
165 sp_long_name='temperature', sp_unit='degC s-1')
166 if ( associated(tv%S) ) call set_up_ale_sponge_field(s, g, gv, tv%S, acsp, 'salt', &
167 sp_long_name='salinity', sp_unit='g kg-1 s-1')
168
169 if (sponge_uv) then
170 u1(:,:,:) = 0.0 ; v1(:,:,:) = 0.0
171 call set_up_ale_sponge_vel_field(u1, v1, g, gv, u, v, acsp)
172 endif
173
174
175 else ! layer mode
176
177 !read eta
178 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
179
180 ! Set the sponge damping rates so that the model will know where to
181 ! apply the sponges, along with the interface heights.
182 call initialize_sponge(idamp, eta, g, pf, csp, gv)
183
184 if ( gv%nkml>0 ) then
185 ! This call to set_up_sponge_ML_density registers the target values of the
186 ! mixed layer density, which is used in determining which layers can be
187 ! inflated without causing static instabilities.
188 do i=is,ie ; pres(i) = tv%P_Ref ; enddo
189 eosdom(:) = eos_domain(g%HI)
190 do j=js,je
191 call calculate_density(t(:,j,1), s(:,j,1), pres, rho(:,j), tv%eqn_of_state, eosdom)
192 enddo
193
194 call set_up_sponge_ml_density(rho, g, csp)
195 endif
196
197 ! Apply sponge in tracer fields
198 call set_up_sponge_field(t, tv%T, g, gv, nz, csp)
199 call set_up_sponge_field(s, tv%S, g, gv, nz, csp)
200
201 endif
202
203end subroutine rgc_initialize_sponges
204
205end module rgc_initialization