dumbbell_surface_forcing.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!> Surface forcing for the dumbbell test case
7
10use mom_domains, only : pass_var, pass_vector, agrid
11use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
12use mom_file_parser, only : get_param, param_file_type, log_version
13use mom_forcing_type, only : forcing, allocate_forcing_type
14use mom_grid, only : ocean_grid_type
15use mom_safe_alloc, only : safe_alloc_ptr
16use mom_time_manager, only : time_type, operator(+), operator(/), get_time
20use mom_variables, only : surface
21
22implicit none ; private
23
25
26!> Control structure for the dumbbell test case forcing
27type, public :: dumbbell_surface_forcing_cs ; private
28 logical :: use_temperature !< If true, temperature and salinity are used as state variables.
29 logical :: restorebuoy !< If true, use restoring surface buoyancy forcing.
30 real :: g_earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
31 real :: flux_const !< The restoring rate at the surface [R Z T-1 ~> kg m-2 s-1].
32! real :: gust_const !< A constant unresolved background gustiness
33! !! that contributes to ustar [R L Z T-2 ~> Pa].
34 real :: slp_amplitude !< The amplitude of pressure loading [R L2 T-2 ~> Pa] applied
35 !! to the reservoirs
36 real :: slp_period !< Period of sinusoidal pressure wave [days]
37 real, dimension(:,:), allocatable :: &
38 forcing_mask !< A mask regulating where forcing occurs [nondim]
39 real, dimension(:,:), allocatable :: &
40 s_restore !< The surface salinity field toward which to restore [S ~> ppt].
41 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the
42 !! timing of diagnostic output.
44
45contains
46
47!> Surface buoyancy (heat and fresh water) fluxes for the dumbbell test case
48subroutine dumbbell_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
49 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
50 !! describe the surface state of the ocean.
51 type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
52 !! possible forcing fields. Unused fields
53 !! have NULL ptrs.
54 type(time_type), intent(in) :: day !< Time of the fluxes.
55 real, intent(in) :: dt !< The amount of time over which
56 !! the fluxes apply [T ~> s]
57 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
58 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
59 type(dumbbell_surface_forcing_cs), pointer :: cs !< A control structure returned by a previous
60 !! call to dumbbell_surface_forcing_init
61 ! Local variables
62 integer :: i, j, is, ie, js, je
63 integer :: isd, ied, jsd, jed
64
65 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
66 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
67
68
69 ! Allocate and zero out the forcing arrays, as necessary.
70 if (cs%use_temperature) then
71 call safe_alloc_ptr(fluxes%evap, isd, ied, jsd, jed)
72 call safe_alloc_ptr(fluxes%lprec, isd, ied, jsd, jed)
73 call safe_alloc_ptr(fluxes%fprec, isd, ied, jsd, jed)
74 call safe_alloc_ptr(fluxes%lrunoff, isd, ied, jsd, jed)
75 call safe_alloc_ptr(fluxes%frunoff, isd, ied, jsd, jed)
76 call safe_alloc_ptr(fluxes%vprec, isd, ied, jsd, jed)
77
78 call safe_alloc_ptr(fluxes%sw, isd, ied, jsd, jed)
79 call safe_alloc_ptr(fluxes%lw, isd, ied, jsd, jed)
80 call safe_alloc_ptr(fluxes%latent, isd, ied, jsd, jed)
81 call safe_alloc_ptr(fluxes%sens, isd, ied, jsd, jed)
82 else ! This is the buoyancy only mode.
83 call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
84 endif
85
86
87 ! MODIFY THE CODE IN THE FOLLOWING LOOPS TO SET THE BUOYANCY FORCING TERMS.
88
89 if ( cs%use_temperature ) then
90 ! Set whichever fluxes are to be used here. Any fluxes that
91 ! are always zero do not need to be changed here.
92 do j=js,je ; do i=is,ie
93 ! Fluxes of fresh water through the surface are in units of [R Z T-1 ~> kg m-2 s-1]
94 ! and are positive downward - i.e. evaporation should be negative.
95 fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)
96 fluxes%lprec(i,j) = 0.0 * g%mask2dT(i,j)
97
98 ! vprec will be set later, if it is needed for salinity restoring.
99 fluxes%vprec(i,j) = 0.0
100
101 ! Heat fluxes are in units of [Q R Z T-1 ~> W m-2] and are positive into the ocean.
102 fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)
103 fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)
104 fluxes%sens(i,j) = 0.0 * g%mask2dT(i,j)
105 fluxes%sw(i,j) = 0.0 * g%mask2dT(i,j)
106 enddo ; enddo
107 else ! This is the buoyancy only mode.
108 do j=js,je ; do i=is,ie
109 ! fluxes%buoy is the buoyancy flux into the ocean [L2 T-3 ~> m2 s-3]. A positive
110 ! buoyancy flux is of the same sign as heating the ocean.
111 fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
112 enddo ; enddo
113 endif
114
115 if (cs%use_temperature .and. cs%restorebuoy) then
116 do j=js,je ; do i=is,ie
117 if (cs%forcing_mask(i,j)>0.) then
118 fluxes%vprec(i,j) = - (g%mask2dT(i,j) * cs%Flux_const) * &
119 ((cs%S_restore(i,j) - sfc_state%SSS(i,j)) / (0.5 * (cs%S_restore(i,j) + sfc_state%SSS(i,j))))
120
121 endif
122 enddo ; enddo
123 endif
124
125end subroutine dumbbell_buoyancy_forcing
126
127!> Dynamic forcing for the dumbbell test case
128subroutine dumbbell_dynamic_forcing(sfc_state, fluxes, day, dt, G, US, CS)
129 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
130 !! describe the surface state of the ocean.
131 type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
132 !! possible forcing fields. Unused fields
133 !! have NULL ptrs.
134 type(time_type), intent(in) :: day !< Time of the fluxes.
135 real, intent(in) :: dt !< The amount of time over which
136 !! the fluxes apply [T ~> s]
137 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
138 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
139 type(dumbbell_surface_forcing_cs), pointer :: cs !< A control structure returned by a previous
140 !! call to dumbbell_surface_forcing_init
141 ! Local variables
142 integer :: i, j, is, ie, js, je
143 integer :: isd, ied, jsd, jed
144 integer :: idays, isecs
145 real :: deg_rad ! A conversion factor from degrees to radians [nondim]
146 real :: rdays ! The elapsed time [days]
147
148
149 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
150 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
151
152 deg_rad = atan(1.0)*4.0/180.
153
154 call get_time(day,isecs,idays)
155 rdays = real(idays) + real(isecs)/8.64e4
156 ! This could be: rdays = time_type_to_real(day)/8.64e4
157
158 ! Allocate and zero out the forcing arrays, as necessary.
159 call safe_alloc_ptr(fluxes%p_surf, isd, ied, jsd, jed)
160 call safe_alloc_ptr(fluxes%p_surf_full, isd, ied, jsd, jed)
161
162 do j=js,je ; do i=is,ie
163 fluxes%p_surf(i,j) = cs%forcing_mask(i,j)* cs%slp_amplitude * &
164 g%mask2dT(i,j) * sin(deg_rad*(rdays/cs%slp_period))
165 fluxes%p_surf_full(i,j) = cs%forcing_mask(i,j) * cs%slp_amplitude * &
166 g%mask2dT(i,j) * sin(deg_rad*(rdays/cs%slp_period))
167 enddo ; enddo
168
169end subroutine dumbbell_dynamic_forcing
170
171!> Reads and sets up the forcing for the dumbbell test case
172subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS)
173 type(time_type), intent(in) :: time !< The current model time.
174 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
175 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
176 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
177 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to
178 !! regulate diagnostic output.
180 pointer :: cs !< A pointer to the control structure for this module
181 ! Local variables
182 real :: s_surf ! Initial surface salinity [S ~> ppt]
183 real :: s_range ! Range of the initial vertical distribution of salinity [S ~> ppt]
184 real :: x ! Latitude normalized by the domain size [nondim]
185 real :: rho0 ! The density used in the Boussinesq approximation [R ~> kg m-3]
186 real :: rho_restore ! The density that is used to convert piston velocities into salt
187 ! or heat fluxes with salinity or temperature restoring [R ~> kg m-3]
188 integer :: i, j
189 logical :: dbrotate ! If true, rotate the domain.
190# include "version_variable.h"
191 character(len=40) :: mdl = "dumbbell_surface_forcing" ! This module's name.
192
193 if (associated(cs)) then
194 call mom_error(warning, "dumbbell_surface_forcing_init called with an associated "// &
195 "control structure.")
196 return
197 endif
198 allocate(cs)
199 cs%diag => diag
200
201 ! Read all relevant parameters and write them to the model log.
202 call log_version(param_file, mdl, version, "")
203 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
204 "If true, Temperature and salinity are used as state variables.", default=.true.)
205
206 call get_param(param_file, mdl, "G_EARTH", cs%G_Earth, &
207 "The gravitational acceleration of the Earth.", &
208 units="m s-2", default=9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
209 call get_param(param_file, mdl, "RHO_0", rho0, &
210 "The mean ocean density used with BOUSSINESQ true to "//&
211 "calculate accelerations and the mass for conservation "//&
212 "properties, or with BOUSSINESQ false to convert some "//&
213 "parameters from vertical units of m to kg m-2.", &
214 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
215 call get_param(param_file, mdl, "DUMBBELL_SLP_AMP", cs%slp_amplitude, &
216 "Amplitude of SLP forcing in reservoirs.", &
217 units="Pa", default=10000.0, scale=us%Pa_to_RL2_T2)
218 call get_param(param_file, mdl, "DUMBBELL_SLP_PERIOD", cs%slp_period, &
219 "Periodicity of SLP forcing in reservoirs.", &
220 units="days", default=1.0)
221 call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
222 'Logical for rotation of dumbbell domain.',&
223 default=.false., do_not_log=.true.)
224 call get_param(param_file, mdl,"INITIAL_SSS", s_surf, &
225 "Initial surface salinity", &
226 units="ppt", default=34.0, scale=us%ppt_to_S, do_not_log=.true.)
227 call get_param(param_file, mdl,"INITIAL_S_RANGE", s_range, &
228 "Initial salinity range (bottom - surface)", &
229 units="ppt", default=2., scale=us%ppt_to_S, do_not_log=.true.)
230
231 call get_param(param_file, mdl, "RESTOREBUOY", cs%restorebuoy, &
232 "If true, the buoyancy fluxes drive the model back "//&
233 "toward some specified surface state with a rate "//&
234 "given by FLUXCONST.", default=.false.)
235 if (cs%restorebuoy) then
236 call get_param(param_file, mdl, "FLUXCONST", cs%Flux_const, &
237 "The constant that relates the restoring surface fluxes to the relative "//&
238 "surface anomalies (akin to a piston velocity). Note the non-MKS units.", &
239 default=0.0, units="m day-1", scale=us%m_to_Z*us%T_to_s)
240 call get_param(param_file, mdl, "RESTORE_FLUX_RHO", rho_restore, &
241 "The density that is used to convert piston velocities into salt or heat "//&
242 "fluxes with RESTORE_SALINITY or RESTORE_TEMPERATURE.", &
243 units="kg m-3", default=rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R, &
244 do_not_log=(cs%Flux_const==0.0))
245 ! Convert FLUXCONST from m day-1 to m s-1 and Flux_const to [R Z T-1 ~> kg m-2 s-1]
246 cs%Flux_const = rho_restore * (cs%Flux_const / 86400.0)
247
248
249 allocate(cs%forcing_mask(g%isd:g%ied, g%jsd:g%jed), source=0.0)
250 allocate(cs%S_restore(g%isd:g%ied, g%jsd:g%jed))
251
252 do j=g%jsc,g%jec
253 do i=g%isc,g%iec
254 ! Compute normalized zonal coordinates (x,y=0 at center of domain)
255! x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5
256! y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5
257 if (dbrotate) then
258 ! This is really y in the rotated case
259 x = ( g%geoLatT(i,j) - g%south_lat ) / g%len_lat - 0.5
260 else
261 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon - 0.5
262 endif
263 cs%forcing_mask(i,j)=0
264 cs%S_restore(i,j) = s_surf
265 if ((x>0.25)) then
266 cs%forcing_mask(i,j) = 1
267 cs%S_restore(i,j) = s_surf + s_range
268 elseif ((x<-0.25)) then
269 cs%forcing_mask(i,j) = 1
270 cs%S_restore(i,j) = s_surf - s_range
271 endif
272 enddo
273 enddo
274 endif
275
277