MOM6
Neverland_surface_forcing.F90
Go to the documentation of this file.
1 !> Wind and buoyancy forcing for the Neverland configurations
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
7 use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr
8 use mom_domains, only : pass_var, pass_vector, agrid
9 use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
13 use mom_grid, only : ocean_grid_type
14 use mom_io, only : file_exists, read_data, slasher
15 use mom_time_manager, only : time_type, operator(+), operator(/)
17 use mom_variables, only : surface
18 
19 implicit none ; private
20 
24 
25 !> This control structure should be used to store any run-time variables
26 !! associated with the Neverland forcing.
27 !!
28 !! It can be readily modified for a specific case, and because it is private there
29 !! will be no changes needed in other code (although they will have to be recompiled).
30 type, public :: neverland_surface_forcing_cs ; private
31 
32  logical :: use_temperature !< If true, use temperature and salinity.
33  logical :: restorebuoy !< If true, use restoring surface buoyancy forcing.
34  real :: rho0 !< The density used in the Boussinesq approximation [R ~> kg m-3].
35  real :: g_earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
36  real :: flux_const !< The restoring rate at the surface [Z T-1 ~> m s-1].
37  real, dimension(:,:), pointer :: &
38  buoy_restore(:,:) => null() !< The pattern to restore buoyancy to.
39  character(len=200) :: inputdir !< The directory where NetCDF input files are.
40  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
41  !! timing of diagnostic output.
42  logical :: first_call = .true. !< True until Neverland_buoyancy_forcing has been called
44 
45 contains
46 
47 !> Sets the surface wind stresses, forces%taux and forces%tauy for the
48 !! Neverland forcing configuration.
49 subroutine neverland_wind_forcing(sfc_state, forces, day, G, US, CS)
50  type(surface), intent(inout) :: sfc_state !< A structure containing fields that
51  !! describe the surface state of the ocean.
52  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
53  type(time_type), intent(in) :: day !< Time used for determining the fluxes.
54  type(ocean_grid_type), intent(inout) :: g !< Grid structure.
55  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
56  type(neverland_surface_forcing_cs), pointer :: cs !< Control structure for this module.
57 
58  ! Local variables
59  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
60  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
61  real :: x, y
62  real :: pi
63  real :: tau_max ! The magnitude of the wind stress [R Z L T-2 ~> Pa]
64  real :: off
65 
66  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
67  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
68  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
69  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
70 
71  ! Allocate the forcing arrays, if necessary.
72  call allocate_mech_forcing(g, forces, stress=.true.)
73 
74  ! Set the surface wind stresses, in units of Pa. A positive taux
75  ! accelerates the ocean to the (pseudo-)east.
76 
77  ! The i-loop extends to is-1 so that taux can be used later in the
78  ! calculation of ustar - otherwise the lower bound would be Isq.
79  pi = 4.0*atan(1.0)
80  forces%taux(:,:) = 0.0
81  tau_max = 0.2 * us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z
82  off = 0.02
83  do j=js,je ; do i=is-1,ieq
84 ! x = (G%geoLonT(i,j)-G%west_lon)/G%len_lon
85  y = (g%geoLatT(i,j)-g%south_lat)/g%len_lat
86 ! forces%taux(I,j) = G%mask2dCu(I,j) * 0.0
87 
88  if (y <= 0.29) then
89  forces%taux(i,j) = forces%taux(i,j) + tau_max * ( (1/0.29)*y - ( 1/(2*pi) )*sin( (2*pi*y) / 0.29 ) )
90  endif
91  if ((y > 0.29) .and. (y <= (0.8-off))) then
92  forces%taux(i,j) = forces%taux(i,j) + tau_max *(0.35+0.65*cos(pi*(y-0.29)/(0.51-off)) )
93  endif
94  if ((y > (0.8-off)) .and. (y <= (1-off))) then
95  forces%taux(i,j) = forces%taux(i,j) + tau_max *( 1.5*( (y-1+off) - (0.1/pi)*sin(10.0*pi*(y-0.8+off)) ) )
96  endif
97  enddo ; enddo
98 
99  do j=js-1,jeq ; do i=is,ie
100  forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0 ! Change this to the desired expression.
101  enddo ; enddo
102 
103  ! Set the surface friction velocity, in units of m s-1. ustar
104  ! is always positive.
105 ! if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
106 ! ! This expression can be changed if desired, but need not be.
107 ! forces%ustar(i,j) = G%mask2dT(i,j) * sqrt((CS%gust_const + &
108 ! sqrt(0.5*(forces%taux(I-1,j)**2 + forces%taux(I,j)**2) + &
109 ! 0.5*(forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2))) * &
110 ! (US%L_to_Z / CS%Rho0) )
111 ! enddo ; enddo ; endif
112 
113 end subroutine neverland_wind_forcing
114 
115 !> Returns the value of a cosine-bell function evaluated at x/L
116 real function cosbell(x,L)
117 
118  real , intent(in) :: x !< non-dimensional position
119  real , intent(in) :: l !< non-dimensional width
120  real :: pi !< 3.1415926... calculated as 4*atan(1)
121 
122  pi = 4.0*atan(1.0)
123  cosbell = 0.5 * (1 + cos(pi*min(abs(x/l),1.0)))
124 end function cosbell
125 
126 !> Returns the value of a sin-spike function evaluated at x/L
127 real function spike(x,L)
128 
129  real , intent(in) :: x !< non-dimensional position
130  real , intent(in) :: l !< non-dimensional width
131  real :: pi !< 3.1415926... calculated as 4*atan(1)
132 
133  pi = 4.0*atan(1.0)
134  spike = (1 - sin(pi*min(abs(x/l),0.5)))
135 end function spike
136 
137 
138 !> Surface fluxes of buoyancy for the Neverland configurations.
139 subroutine neverland_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
140  type(surface), intent(inout) :: sfc_state !< A structure containing fields that
141  !! describe the surface state of the ocean.
142  type(forcing), intent(inout) :: fluxes !< Forcing fields.
143  type(time_type), intent(in) :: day !< Time used for determining the fluxes.
144  real, intent(in) :: dt !< Forcing time step (s).
145  type(ocean_grid_type), intent(inout) :: g !< Grid structure.
146  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
147  type(neverland_surface_forcing_cs), pointer :: cs !< Control structure for this module.
148  ! Local variables
149  real :: buoy_rest_const ! A constant relating density anomalies to the
150  ! restoring buoyancy flux [L2 T-3 R-1 ~> m5 s-3 kg-1].
151  real :: density_restore ! De
152  integer :: i, j, is, ie, js, je
153  integer :: isd, ied, jsd, jed
154 
155  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
156  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
157 
158  ! When modifying the code, comment out this error message. It is here
159  ! so that the original (unmodified) version is not accidentally used.
160 
161  ! Allocate and zero out the forcing arrays, as necessary. This portion is
162  ! usually not changed.
163  if (cs%use_temperature) then
164  call mom_error(fatal, "Neverland_buoyancy_forcing: " // &
165  "Temperature and salinity mode not coded!" )
166  else
167  ! This is the buoyancy only mode.
168  call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
169  endif
170 
171 
172  ! MODIFY THE CODE IN THE FOLLOWING LOOPS TO SET THE BUOYANCY FORCING TERMS.
173  if (cs%restorebuoy .and. cs%first_call) then
174  call safe_alloc_ptr(cs%buoy_restore, isd, ied, jsd, jed)
175  cs%first_call = .false.
176  ! Set CS%buoy_restore(i,j) here
177  endif
178 
179  if ( cs%use_temperature ) then
180  call mom_error(fatal, "Neverland_buoyancy_surface_forcing: " // &
181  "Temperature/salinity restoring not coded!" )
182  else ! This is the buoyancy only mode.
183  do j=js,je ; do i=is,ie
184  ! fluxes%buoy is the buoyancy flux into the ocean [L2 T-3 ~> m2 s-3]. A positive
185  ! buoyancy flux is of the same sign as heating the ocean.
186  fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
187  enddo ; enddo
188  endif
189 
190  if (cs%restorebuoy) then
191  if (cs%use_temperature) then
192  call mom_error(fatal, "Neverland_buoyancy_surface_forcing: " // &
193  "Temperature/salinity restoring not coded!" )
194  else
195  ! When modifying the code, comment out this error message. It is here
196  ! so that the original (unmodified) version is not accidentally used.
197 
198  ! The -1 is because density has the opposite sign to buoyancy.
199  buoy_rest_const = -1.0 * (cs%G_Earth * cs%Flux_const) / cs%Rho0
200  do j=js,je ; do i=is,ie
201  ! Set density_restore to an expression for the surface potential
202  ! density [kg m-3] that is being restored toward.
203  density_restore = 1030.0
204 
205  fluxes%buoy(i,j) = g%mask2dT(i,j) * buoy_rest_const * &
206  us%kg_m3_to_R*(density_restore - sfc_state%sfc_density(i,j))
207  enddo ; enddo
208  endif
209  endif ! end RESTOREBUOY
210 
211 end subroutine neverland_buoyancy_forcing
212 
213 !> Initializes the Neverland control structure.
214 subroutine neverland_surface_forcing_init(Time, G, US, param_file, diag, CS)
215  type(time_type), intent(in) :: time !< The current model time.
216  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
217  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
218  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to parse for
219  !! model parameter values.
220  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate diagnostic output.
221  type(neverland_surface_forcing_cs), pointer :: cs !< A pointer that is set to point to the control structure
222  !! for this module
223  ! This include declares and sets the variable "version".
224 #include "version_variable.h"
225  ! Local variables
226  character(len=40) :: mdl = "Neverland_surface_forcing" ! This module's name.
227 
228  if (associated(cs)) then
229  call mom_error(warning, "Neverland_surface_forcing_init called with an associated "// &
230  "control structure.")
231  return
232  endif
233  allocate(cs)
234  cs%diag => diag
235 
236  ! Read all relevant parameters and write them to the model log.
237  call log_version(param_file, mdl, version, "")
238  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
239  "If true, Temperature and salinity are used as state "//&
240  "variables.", default=.true.)
241 
242  call get_param(param_file, mdl, "G_EARTH", cs%G_Earth, &
243  "The gravitational acceleration of the Earth.", &
244  units="m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
245  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
246  "The mean ocean density used with BOUSSINESQ true to "//&
247  "calculate accelerations and the mass for conservation "//&
248  "properties, or with BOUSSINSEQ false to convert some "//&
249  "parameters from vertical units of m to kg m-2.", &
250  units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
251 ! call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
252 ! "The background gustiness in the winds.", units="Pa", &
253 ! default=0.02)
254 
255  call get_param(param_file, mdl, "RESTOREBUOY", cs%restorebuoy, &
256  "If true, the buoyancy fluxes drive the model back "//&
257  "toward some specified surface state with a rate "//&
258  "given by FLUXCONST.", default= .false.)
259 
260  if (cs%restorebuoy) then
261  call get_param(param_file, mdl, "FLUXCONST", cs%flux_const, &
262  "The constant that relates the restoring surface fluxes "//&
263  "to the relative surface anomalies (akin to a piston "//&
264  "velocity). Note the non-MKS units.", &
265  units="m day-1", scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
266  ! Convert CS%flux_const from m day-1 to m s-1.
267  cs%flux_const = cs%flux_const / 86400.0
268  endif
269 
270 end subroutine neverland_surface_forcing_init
271 
272 end module neverland_surface_forcing
neverland_surface_forcing::spike
real function spike(x, L)
Returns the value of a sin-spike function evaluated at x/L.
Definition: Neverland_surface_forcing.F90:128
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
mom_forcing_type::allocate_forcing_type
subroutine, public allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, iceberg, salt)
Conditionally allocate fields within the forcing type.
Definition: MOM_forcing_type.F90:2811
mom_forcing_type::mech_forcing
Structure that contains pointers to the mechanical forcing at the surface used to drive the liquid oc...
Definition: MOM_forcing_type.F90:187
mom_variables::surface
Pointers to various fields which may be used describe the surface state of MOM, and which will be ret...
Definition: MOM_variables.F90:38
neverland_surface_forcing::neverland_buoyancy_forcing
subroutine, public neverland_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
Surface fluxes of buoyancy for the Neverland configurations.
Definition: Neverland_surface_forcing.F90:140
neverland_surface_forcing::neverland_surface_forcing_cs
This control structure should be used to store any run-time variables associated with the Neverland f...
Definition: Neverland_surface_forcing.F90:30
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Call this subroutine to determine whether the averaging is currently enabled. .true....
Definition: MOM_diag_mediator.F90:1850
neverland_surface_forcing::neverland_surface_forcing_init
subroutine, public neverland_surface_forcing_init(Time, G, US, param_file, diag, CS)
Initializes the Neverland control structure.
Definition: Neverland_surface_forcing.F90:215
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_forcing_type::allocate_mech_forcing
subroutine, public allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg)
Conditionally allocate fields within the mechanical forcing type.
Definition: MOM_forcing_type.F90:2879
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
neverland_surface_forcing
Wind and buoyancy forcing for the Neverland configurations.
Definition: Neverland_surface_forcing.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
neverland_surface_forcing::cosbell
real function cosbell(x, L)
Returns the value of a cosine-bell function evaluated at x/L.
Definition: Neverland_surface_forcing.F90:117
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_domains::pass_vector
Do a halo update on a pair of arrays representing the two components of a vector.
Definition: MOM_domains.F90:54
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_forcing_type::forcing
Structure that contains pointers to the boundary forcing used to drive the liquid ocean simulated by ...
Definition: MOM_forcing_type.F90:50
mom_io::file_exists
Indicate whether a file exists, perhaps with domain decomposition.
Definition: MOM_io.F90:68
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
neverland_surface_forcing::neverland_wind_forcing
subroutine, public neverland_wind_forcing(sfc_state, forces, day, G, US, CS)
Sets the surface wind stresses, forcestaux and forcestauy for the Neverland forcing configuration.
Definition: Neverland_surface_forcing.F90:50
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239