BFB_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 boundary-forced-basin (BFB) configuration
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(/)
20use mom_variables, only : surface
21
22implicit none ; private
23
25
26!> Control structure for BFB_surface_forcing
27type, public :: bfb_surface_forcing_cs ; private
28
29 logical :: use_temperature !< If true, temperature and salinity are used as state variables.
30 logical :: restorebuoy !< If true, use restoring surface buoyancy forcing.
31 real :: rho0 !< The density used in the Boussinesq approximation [R ~> kg m-3].
32 real :: g_earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
33 real :: flux_const !< The restoring rate at the surface [Z T-1 ~> m s-1].
34 real :: rho_restore !< The density that is used to convert piston velocities into salt
35 !! or heat fluxes with salinity or temperature restoring [R ~> kg m-3]
36 real :: sst_s !< SST at the southern edge of the linear forcing ramp [C ~> degC]
37 real :: sst_n !< SST at the northern edge of the linear forcing ramp [C ~> degC]
38 real :: s_ref !< Reference salinity used throughout the domain [S ~> ppt]
39 real :: lfrslat !< Southern latitude where the linear forcing ramp begins [degrees_N] or [km]
40 real :: lfrnlat !< Northern latitude where the linear forcing ramp ends [degrees_N] or [km]
41 real :: rho_t0_s0 !< The density at T=0, S=0 [R ~> kg m-3]
42 real :: drho_dt !< The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
43 real :: drho_ds !< The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
44 !! Note that temperature and salinity are being used as dummy variables here.
45 !! All temperatures are converted into density.
46
47 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
48 !! regulate the timing of diagnostic output.
50
51contains
52
53!> Bouyancy forcing for the boundary-forced-basin (BFB) configuration
54subroutine bfb_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
55 type(surface), intent(inout) :: sfc_state !< A structure containing fields that
56 !! describe the surface state of the ocean.
57 type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
58 !! possible forcing fields. Unused fields
59 !! have NULL ptrs.
60 type(time_type), intent(in) :: day !< Time of the fluxes.
61 real, intent(in) :: dt !< The amount of time over which
62 !! the fluxes apply [T ~> s]
63 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
64 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
65 type(bfb_surface_forcing_cs), pointer :: cs !< A pointer to the control structure
66 !! returned by a previous call to
67 !! BFB_surface_forcing_init.
68 ! Local variables
69 real :: temp_restore ! The temperature that is being restored toward [C ~> degC].
70 real :: salin_restore ! The salinity that is being restored toward [S ~> ppt].
71 real :: density_restore ! The potential density that is being restored
72 ! toward [R ~> kg m-3].
73 real :: rhoxcp ! Reference density times heat capacity times unit scaling
74 ! factors [Q R C-1 ~> J m-3 degC-1]
75 real :: buoy_rest_const ! A constant relating density anomalies to the
76 ! restoring buoyancy flux [L2 T-3 R-1 ~> m5 s-3 kg-1].
77 integer :: i, j, is, ie, js, je
78 integer :: isd, ied, jsd, jed
79
80 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
81 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
82
83 ! Allocate and zero out the forcing arrays, as necessary. This portion is
84 ! usually not changed.
85 if (cs%use_temperature) then
86 call safe_alloc_ptr(fluxes%evap, isd, ied, jsd, jed)
87 call safe_alloc_ptr(fluxes%lprec, isd, ied, jsd, jed)
88 call safe_alloc_ptr(fluxes%fprec, isd, ied, jsd, jed)
89 call safe_alloc_ptr(fluxes%lrunoff, isd, ied, jsd, jed)
90 call safe_alloc_ptr(fluxes%frunoff, isd, ied, jsd, jed)
91 call safe_alloc_ptr(fluxes%vprec, isd, ied, jsd, jed)
92
93 call safe_alloc_ptr(fluxes%sw, isd, ied, jsd, jed)
94 call safe_alloc_ptr(fluxes%lw, isd, ied, jsd, jed)
95 call safe_alloc_ptr(fluxes%latent, isd, ied, jsd, jed)
96 call safe_alloc_ptr(fluxes%sens, isd, ied, jsd, jed)
97 else ! This is the buoyancy only mode.
98 call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
99 endif
100
101 if ( cs%use_temperature ) then
102 ! Set whichever fluxes are to be used here. Any fluxes that
103 ! are always zero do not need to be changed here.
104 do j=js,je ; do i=is,ie
105 ! Fluxes of fresh water through the surface are in units of [R Z T-1 ~> kg m-2 s-1]
106 ! and are positive downward - i.e. evaporation should be negative.
107 fluxes%evap(i,j) = -0.0 * g%mask2dT(i,j)
108 fluxes%lprec(i,j) = 0.0 * g%mask2dT(i,j)
109
110 ! vprec will be set later, if it is needed for salinity restoring.
111 fluxes%vprec(i,j) = 0.0
112
113 ! Heat fluxes are in units of [Q R Z T-1 ~> W m-2] and are positive into the ocean.
114 fluxes%lw(i,j) = 0.0 * g%mask2dT(i,j)
115 fluxes%latent(i,j) = 0.0 * g%mask2dT(i,j)
116 fluxes%sens(i,j) = 0.0 * g%mask2dT(i,j)
117 fluxes%sw(i,j) = 0.0 * g%mask2dT(i,j)
118 enddo ; enddo
119 else ! This is the buoyancy only mode.
120 do j=js,je ; do i=is,ie
121 ! fluxes%buoy is the buoyancy flux into the ocean [L2 T-3 ~> m2 s-3]. A positive
122 ! buoyancy flux is of the same sign as heating the ocean.
123 fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
124 enddo ; enddo
125 endif
126
127 if (cs%restorebuoy) then
128 if (cs%use_temperature) then
129 call safe_alloc_ptr(fluxes%heat_added, isd, ied, jsd, jed)
130 ! When modifying the code, comment out this error message. It is here
131 ! so that the original (unmodified) version is not accidentally used.
132 call mom_error(fatal, "User_buoyancy_surface_forcing: " // &
133 "Temperature and salinity restoring used without modification." )
134
135 rhoxcp = cs%rho_restore * fluxes%C_p
136 do j=js,je ; do i=is,ie
137 ! Set Temp_restore and Salin_restore to the temperature (in [C ~> degC]) and
138 ! salinity (in [S ~> ppt]) that are being restored toward.
139 temp_restore = 0.0
140 salin_restore = 0.0
141
142 fluxes%heat_added(i,j) = (g%mask2dT(i,j) * (rhoxcp * cs%Flux_const)) * &
143 (temp_restore - sfc_state%SST(i,j))
144 fluxes%vprec(i,j) = - (g%mask2dT(i,j) * (cs%rho_restore*cs%Flux_const)) * &
145 ((salin_restore - sfc_state%SSS(i,j)) / (0.5 * (salin_restore + sfc_state%SSS(i,j))))
146 enddo ; enddo
147 else
148 ! When modifying the code, comment out this error message. It is here
149 ! so that the original (unmodified) version is not accidentally used.
150 ! call MOM_error(FATAL, "User_buoyancy_surface_forcing: " // &
151 ! "Buoyancy restoring used without modification." )
152
153 ! The -1 is because density has the opposite sign to buoyancy.
154 buoy_rest_const = -1.0 * (cs%G_Earth * cs%Flux_const) / cs%rho_restore
155 temp_restore = 0.0
156 do j=js,je ; do i=is,ie
157 ! Set density_restore to an expression for the surface potential
158 ! density [R ~> kg m-3] that is being restored toward.
159 if (g%geoLatT(i,j) < cs%lfrslat) then
160 temp_restore = cs%SST_s
161 elseif (g%geoLatT(i,j) > cs%lfrnlat) then
162 temp_restore = cs%SST_n
163 else
164 temp_restore = (cs%SST_s - cs%SST_n)/(cs%lfrslat - cs%lfrnlat) * &
165 (g%geoLatT(i,j) - cs%lfrslat) + cs%SST_s
166 endif
167
168 density_restore = (cs%Rho_T0_S0 + cs%dRho_dS*cs%S_ref) + cs%dRho_dT*temp_restore
169 fluxes%buoy(i,j) = g%mask2dT(i,j) * buoy_rest_const * &
170 (density_restore - sfc_state%sfc_density(i,j))
171 enddo ; enddo
172 endif
173 endif ! end RESTOREBUOY
174
175end subroutine bfb_buoyancy_forcing
176
177!> Initialization for forcing the boundary-forced-basin (BFB) configuration
178subroutine bfb_surface_forcing_init(Time, G, US, param_file, diag, CS)
179 type(time_type), intent(in) :: time !< The current model time.
180 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
181 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
182 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
183 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to
184 !! regulate diagnostic output.
185 type(bfb_surface_forcing_cs), pointer :: cs !< A pointer to the control structure for this module
186
187 ! This include declares and sets the variable "version".
188# include "version_variable.h"
189 character(len=40) :: mdl = "BFB_surface_forcing" ! This module's name.
190
191 if (associated(cs)) then
192 call mom_error(warning, "BFB_surface_forcing_init called with an associated "// &
193 "control structure.")
194 return
195 endif
196 allocate(cs)
197 cs%diag => diag
198
199 ! Read all relevant parameters and write them to the model log.
200 call log_version(param_file, mdl, version, "")
201 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
202 "If true, Temperature and salinity are used as state variables.", default=.true.)
203
204 call get_param(param_file, mdl, "G_EARTH", cs%G_Earth, &
205 "The gravitational acceleration of the Earth.", &
206 units="m s-2", default=9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
207 call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
208 "The mean ocean density used with BOUSSINESQ true to "//&
209 "calculate accelerations and the mass for conservation "//&
210 "properties, or with BOUSSINESQ false to convert some "//&
211 "parameters from vertical units of m to kg m-2.", &
212 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
213 call get_param(param_file, mdl, "LFR_SLAT", cs%lfrslat, &
214 "Southern latitude where the linear forcing ramp begins.", &
215 units=g%y_ax_unit_short, default=20.0)
216 call get_param(param_file, mdl, "LFR_NLAT", cs%lfrnlat, &
217 "Northern latitude where the linear forcing ramp ends.", &
218 units=g%y_ax_unit_short, default=40.0)
219 call get_param(param_file, mdl, "SST_S", cs%SST_s, &
220 "SST at the southern edge of the linear forcing ramp.", &
221 units="degC", default=20.0, scale=us%degC_to_C)
222 call get_param(param_file, mdl, "SST_N", cs%SST_n, &
223 "SST at the northern edge of the linear forcing ramp.", &
224 units="degC", default=10.0, scale=us%degC_to_C)
225 call get_param(param_file, mdl, "DRHO_DT", cs%dRho_dT, &
226 "The partial derivative of density with temperature.", &
227 units="kg m-3 K-1", default=-0.2, scale=us%kg_m3_to_R*us%C_to_degC)
228 call get_param(param_file, mdl, "DRHO_DS", cs%dRho_dS, &
229 "The partial derivative of density with salinity.", &
230 units="kg m-3 PSU-1", default=0.8, scale=us%kg_m3_to_R*us%S_to_ppt)
231 call get_param(param_file, mdl, "RHO_T0_S0", cs%Rho_T0_S0, &
232 "The density at T=0, S=0.", units="kg m-3", default=1000.0, scale=us%kg_m3_to_R)
233 call get_param(param_file, mdl, "S_REF", cs%S_ref, &
234 "The reference salinity used here throughout the domain.", &
235 units="PSU", default=35.0, scale=us%ppt_to_S)
236
237 call get_param(param_file, mdl, "RESTOREBUOY", cs%restorebuoy, &
238 "If true, the buoyancy fluxes drive the model back "//&
239 "toward some specified surface state with a rate "//&
240 "given by FLUXCONST.", default=.false.)
241 if (cs%restorebuoy) then
242 call get_param(param_file, mdl, "FLUXCONST", cs%Flux_const, &
243 "The constant that relates the restoring surface fluxes to the relative "//&
244 "surface anomalies (akin to a piston velocity). Note the non-MKS units.", &
245 default=0.0, units="m day-1", scale=us%m_to_Z*us%T_to_s)
246 ! Convert CS%Flux_const from m day-1 to m s-1.
247 cs%Flux_const = cs%Flux_const / 86400.0
248 call get_param(param_file, mdl, "RESTORE_FLUX_RHO", cs%rho_restore, &
249 "The density that is used to convert piston velocities into salt or heat "//&
250 "fluxes with RESTORE_SALINITY or RESTORE_TEMPERATURE.", &
251 units="kg m-3", default=cs%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R, &
252 do_not_log=(cs%Flux_const==0.0))
253 endif
254
255end subroutine bfb_surface_forcing_init
256
257end module bfb_surface_forcing