SCM_CVMix_tests.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!> Initial conditions and forcing for the single column model (SCM) CVMix test set.
7
8use mom_domains, only : pass_var, pass_vector, to_all
9use mom_error_handler, only : mom_error, fatal
10use mom_file_parser, only : get_param, log_version, param_file_type
12use mom_grid, only : ocean_grid_type
14use mom_safe_alloc, only : safe_alloc_ptr
16use mom_time_manager, only : time_type, operator(+), operator(/), time_type_to_real
19
20implicit none ; private
21
22#include <MOM_memory.h>
23
29
30! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34
35!> Container for surface forcing parameters
36type scm_cvmix_tests_cs ; private
37 logical :: usewindstress !< True to use wind stress
38 logical :: useheatflux !< True to use heat flux
39 logical :: useevaporation !< True to use evaporation
40 logical :: usediurnalsw !< True to use diurnal sw radiation
41 real :: tau_x !< (Constant) Wind stress, X [R L Z T-2 ~> Pa]
42 real :: tau_y !< (Constant) Wind stress, Y [R L Z T-2 ~> Pa]
43 real :: surf_hf !< (Constant) Heat flux [C Z T-1 ~> m degC s-1]
44 real :: surf_evap !< (Constant) Evaporation rate [Z T-1 ~> m s-1]
45 real :: max_sw !< maximum of diurnal sw radiation [C Z T-1 ~> degC m s-1]
46 real :: rho0 !< reference density [R ~> kg m-3]
47 real :: rho_restore !< The density that is used to convert piston velocities
48 !! into salt or heat fluxes [R ~> kg m-3]
49end type
50
51! This include declares and sets the variable "version".
52#include "version_variable.h"
53
54character(len=40) :: mdl = "SCM_CVMix_tests" !< This module's name.
55
56contains
57
58!> Initializes temperature and salinity for the SCM CVMix test example
59subroutine scm_cvmix_tests_ts_init(T, S, h, G, GV, US, param_file, just_read)
60 type(ocean_grid_type), intent(in) :: g !< Grid structure
61 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
62 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Potential temperature [C ~> degC]
63 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Salinity [S ~> ppt]
64 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m]
65 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
66 type(param_file_type), intent(in) :: param_file !< Input parameter structure
67 logical, intent(in) :: just_read !< If present and true, this call
68 !! will only read parameters without changing T & S.
69 ! Local variables
70 real :: upperlayertempmld !< Upper layer Temp MLD thickness [Z ~> m].
71 real :: upperlayersaltmld !< Upper layer Salt MLD thickness [Z ~> m].
72 real :: upperlayertemp !< Upper layer temperature (SST if thickness 0) [C ~> degC]
73 real :: upperlayersalt !< Upper layer salinity (SSS if thickness 0) [S ~> ppt]
74 real :: lowerlayertemp !< Temp at top of lower layer [C ~> degC]
75 real :: lowerlayersalt !< Salt at top of lower layer [S ~> ppt]
76 real :: lowerlayerdtdz !< Temp gradient in lower layer [C Z-1 ~> degC m-1].
77 real :: lowerlayerdsdz !< Salt gradient in lower layer [S Z-1 ~> ppt m-1].
78 real :: lowerlayermintemp !< Minimum temperature in lower layer [C ~> degC]
79 real :: zc, dz, top, bottom ! Depths and thicknesses [Z ~> m].
80 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
81
82 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
83 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
84
85 if (.not.just_read) call log_version(param_file, mdl, version)
86 call get_param(param_file, mdl, "SCM_TEMP_MLD", upperlayertempmld, &
87 'Initial temp mixed layer depth', &
88 units='m', default=0.0, scale=us%m_to_Z, do_not_log=just_read)
89 call get_param(param_file, mdl, "SCM_SALT_MLD", upperlayersaltmld, &
90 'Initial salt mixed layer depth', &
91 units='m', default=0.0, scale=us%m_to_Z, do_not_log=just_read)
92 call get_param(param_file, mdl, "SCM_L1_SALT", upperlayersalt, &
93 'Layer 2 surface salinity', units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=just_read)
94 call get_param(param_file, mdl, "SCM_L1_TEMP", upperlayertemp, &
95 'Layer 1 surface temperature', units="degC", default=20.0, scale=us%degC_to_C, do_not_log=just_read)
96 call get_param(param_file, mdl, "SCM_L2_SALT", lowerlayersalt, &
97 'Layer 2 surface salinity', units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=just_read)
98 call get_param(param_file, mdl, "SCM_L2_TEMP", lowerlayertemp, &
99 'Layer 2 surface temperature', units="degC", default=20.0, scale=us%degC_to_C, do_not_log=just_read)
100 call get_param(param_file, mdl, "SCM_L2_DTDZ", lowerlayerdtdz, &
101 'Initial temperature stratification in layer 2', &
102 units='C/m', default=0.0, scale=us%degC_to_C*us%Z_to_m, do_not_log=just_read)
103 call get_param(param_file, mdl, "SCM_L2_DSDZ", lowerlayerdsdz, &
104 'Initial salinity stratification in layer 2', &
105 units='PPT/m', default=0.0, scale=us%ppt_to_S*us%Z_to_m, do_not_log=just_read)
106 call get_param(param_file, mdl, "SCM_L2_MINTEMP",lowerlayermintemp, &
107 'Layer 2 minimum temperature', units="degC", default=4.0, scale=us%degC_to_C, do_not_log=just_read)
108
109 if (just_read) return ! All run-time parameters have been read, so return.
110
111 do j=js,je ; do i=is,ie
112 top = 0. ! Reference to surface
113 bottom = 0.
114 do k=1,nz
115 bottom = bottom - h(i,j,k) ! Interface below layer [Z ~> m]
116 zc = 0.5*( top + bottom ) ! Z of middle of layer [Z ~> m]
117 dz = min(0., zc + upperlayertempmld)
118 t(i,j,k) = max(lowerlayermintemp,lowerlayertemp + lowerlayerdtdz * dz)
119 dz = min(0., zc + upperlayersaltmld)
120 s(i,j,k) = lowerlayersalt + lowerlayerdsdz * dz
121 top = bottom
122 enddo ! k
123 enddo ; enddo
124
125end subroutine scm_cvmix_tests_ts_init
126
127!> Initializes surface forcing for the CVMix test case suite
128subroutine scm_cvmix_tests_surface_forcing_init(Time, G, param_file, CS)
129 type(time_type), intent(in) :: time !< Model time
130 type(ocean_grid_type), intent(in) :: g !< Grid structure
131 type(param_file_type), intent(in) :: param_file !< Input parameter structure
132 type(scm_cvmix_tests_cs), pointer :: cs !< Parameter container
133
134
135 ! This include declares and sets the variable "version".
136# include "version_variable.h"
137 type(unit_scale_type), pointer :: us => null() !< A dimensional unit scaling type
138
139 us => g%US
140
141 if (associated(cs)) then
142 call mom_error(fatal, "SCM_CVMix_tests_surface_forcing_init called with an associated "// &
143 "control structure.")
144 return
145 endif
146 allocate(cs)
147
148 ! Read all relevant parameters and write them to the model log.
149 call log_version(param_file, mdl, version, "")
150 call get_param(param_file, mdl, "SCM_USE_WIND_STRESS", cs%UseWindStress, &
151 "Wind Stress switch used in the SCM CVMix surface forcing.", &
152 default=.false.)
153 call get_param(param_file, mdl, "SCM_USE_HEAT_FLUX", cs%UseHeatFlux, &
154 "Heat flux switch used in the SCM CVMix test surface forcing.", &
155 default=.false.)
156 call get_param(param_file, mdl, "SCM_USE_EVAPORATION", cs%UseEvaporation, &
157 "Evaporation switch used in the SCM CVMix test surface forcing.", &
158 default=.false.)
159 call get_param(param_file, mdl, "SCM_USE_DIURNAL_SW", cs%UseDiurnalSW, &
160 "Diurnal sw radation switch used in the SCM CVMix test surface forcing.", &
161 default=.false.)
162 if (cs%UseWindStress) then
163 call get_param(param_file, mdl, "SCM_TAU_X", cs%tau_x, &
164 "Constant X-dir wind stress used in the SCM CVMix test surface forcing.", &
165 units='N/m2', scale=us%kg_m2s_to_RZ_T*us%m_s_to_L_T, fail_if_missing=.true.)
166 call get_param(param_file, mdl, "SCM_TAU_Y", cs%tau_y, &
167 "Constant y-dir wind stress used in the SCM CVMix test surface forcing.", &
168 units='N/m2', scale=us%kg_m2s_to_RZ_T*us%m_s_to_L_T, fail_if_missing=.true.)
169 endif
170 if (cs%UseHeatFlux) then
171 call get_param(param_file, mdl, "SCM_HEAT_FLUX", cs%surf_HF, &
172 "Constant surface heat flux used in the SCM CVMix test surface forcing.", &
173 units='m K/s', scale=us%m_to_Z*us%degC_to_C*us%T_to_s, fail_if_missing=.true.)
174 endif
175 if (cs%UseEvaporation) then
176 call get_param(param_file, mdl, "SCM_EVAPORATION", cs%surf_evap, &
177 "Constant surface evaporation used in the SCM CVMix test surface forcing.", &
178 units='m/s', scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
179 endif
180 if (cs%UseDiurnalSW) then
181 call get_param(param_file, mdl, "SCM_DIURNAL_SW_MAX", cs%Max_sw, &
182 "Maximum diurnal sw radiation used in the SCM CVMix test surface forcing.", &
183 units='m K/s', scale=us%m_to_Z*us%degC_to_C*us%T_to_s, fail_if_missing=.true.)
184 endif
185 call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
186 "The mean ocean density used with BOUSSINESQ true to "//&
187 "calculate accelerations and the mass for conservation "//&
188 "properties, or with BOUSSINESQ false to convert some "//&
189 "parameters from vertical units of m to kg m-2.", &
190 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
191 call get_param(param_file, mdl, "RESTORE_FLUX_RHO", cs%rho_restore, &
192 "The density that is used to convert piston velocities into salt or heat fluxes.", &
193 units="kg m-3", default=cs%Rho0*us%R_to_kg_m3, scale=us%kg_m3_to_R)
194
196
197subroutine scm_cvmix_tests_wind_forcing(sfc_state, forces, day, G, US, CS)
198 type(surface), intent(in) :: sfc_state !< Surface state structure
199 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
200 type(time_type), intent(in) :: day !< Time in days
201 type(ocean_grid_type), intent(inout) :: g !< Grid structure
202 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
203 type(scm_cvmix_tests_cs), pointer :: cs !< Container for SCM parameters
204 ! Local variables
205 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
206 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
207 real :: mag_tau ! The magnitude of the wind stress [R Z2 T-2 ~> Pa]
208 ! Bounds for loops and memory allocation
209 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
210 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
211 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
212 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
213
214 do j=js,je ; do i=isq,ieq
215 forces%taux(i,j) = cs%tau_x
216 enddo ; enddo
217 do j=jsq,jeq ; do i=is,ie
218 forces%tauy(i,j) = cs%tau_y
219 enddo ; enddo
220 call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)
221
222 mag_tau = us%L_to_Z * sqrt((cs%tau_x*cs%tau_x) + (cs%tau_y*cs%tau_y))
223 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
224 forces%ustar(i,j) = sqrt( mag_tau / cs%Rho0 )
225 enddo ; enddo ; endif
226
227 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
228 forces%tau_mag(i,j) = mag_tau
229 enddo ; enddo ; endif
230
231end subroutine scm_cvmix_tests_wind_forcing
232
233
234subroutine scm_cvmix_tests_buoyancy_forcing(sfc_state, fluxes, day, G, US, CS)
235 type(surface), intent(in) :: sfc_state !< Surface state structure
236 type(forcing), intent(inout) :: fluxes !< Surface fluxes structure
237 type(time_type), intent(in) :: day !< Current model time
238 type(ocean_grid_type), intent(inout) :: g !< Grid structure
239 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
240 type(scm_cvmix_tests_cs), pointer :: cs !< Container for SCM parameters
241
242 ! Local variables
243 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
244 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
245 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
246
247 pi = 4.0*atan(1.0)
248
249 ! Bounds for loops and memory allocation
250 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
251 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
252 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
253 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
254
255 if (cs%UseHeatFlux) then
256 ! Note CVMix test inputs give Heat flux in [Z C T-1 ~> m K s-1]
257 ! therefore must convert to [Q R Z T-1 ~> W m-2] by multiplying
258 ! by Rho0*Cp
259 do j=jsq,jeq ; do i=is,ie
260 fluxes%sens(i,j) = cs%surf_HF * cs%rho_restore * fluxes%C_p
261 enddo ; enddo
262 endif
263
264 if (cs%UseEvaporation) then
265 do j=jsq,jeq ; do i=is,ie
266 ! Note CVMix test inputs give evaporation in [Z T-1 ~> m s-1]
267 ! This therefore must be converted to mass flux in [R Z T-1 ~> kg m-2 s-1]
268 ! by multiplying by density and some unit conversion factors.
269 fluxes%evap(i,j) = cs%surf_evap * cs%rho_restore
270 enddo ; enddo
271 endif
272
273 if (cs%UseDiurnalSW) then
274 do j=jsq,jeq ; do i=is,ie
275 ! Note CVMix test inputs give max sw rad in [Z C T-1 ~> m degC s-1]
276 ! therefore must convert to [Q R Z T-1 ~> W m-2] by multiplying by Rho0*Cp
277 ! Note diurnal cycle peaks at Noon.
278 fluxes%sw(i,j) = cs%Max_sw * max(0.0, cos(2*pi*(time_type_to_real(day)/86400.0 - 0.5))) * &
279 cs%rho_restore * fluxes%C_p
280 enddo ; enddo
281 endif
282
284
285end module scm_cvmix_tests