MOM6
Functions/Subroutines
neverworld_initialization Module Reference

Detailed Description

Initialization for the "Neverworld" configuration.

Functions/Subroutines

subroutine, public neverworld_initialize_topography (D, G, param_file, max_depth)
 This subroutine sets up the Neverworld test case topography. More...
 
real function cosbell (x, L)
 Returns the value of a cosine-bell function evaluated at x/L. More...
 
real function spike (x, L)
 Returns the value of a sin-spike function evaluated at x/L. More...
 
real function cone (x, x0, L, clip)
 Returns the value of a triangular function centered at x=x0 with value 1 and linearly decreasing to 0 at x=x0+/-L, and 0 otherwise. If clip is present the top of the cone is cut off at "clip", which effectively defaults to 1. More...
 
real function scurve (x, x0, L)
 Returns an s-curve s(x) s.t. s(x0)<=0, s(x0+L)>=1 and cubic in between. More...
 
real function cstprof (x, x0, L, lf, bf, sf, sh)
 Returns a "coastal" profile. More...
 
real function dist_line_fixed_x (x, y, x0, y0, y1)
 Distance between points x,y and a line segment (x0,y0) and (x0,y1). More...
 
real function dist_line_fixed_y (x, y, x0, x1, y0)
 Distance between points x,y and a line segment (x0,y0) and (x1,y0). More...
 
real function ns_coast (lon, lat, lon0, lat0, lat1, dlon, sh)
 A "coast profile" applied in an N-S line from lon0,lat0 to lon0,lat1. More...
 
real function ew_coast (lon, lat, lon0, lon1, lat0, dlat, sh)
 A "coast profile" applied in an E-W line from lon0,lat0 to lon1,lat0. More...
 
real function ns_ridge (lon, lat, lon0, lat0, lat1, dlon, rh)
 A NS ridge. More...
 
real function circ_ridge (lon, lat, lon0, lat0, ring_radius, ring_thickness, ridge_height)
 A circular ridge. More...
 
subroutine, public neverworld_initialize_thickness (h, G, GV, US, param_file, eqn_of_state, P_ref)
 This subroutine initializes layer thicknesses for the Neverworld test case, by finding the depths of interfaces in a specified latitude-dependent temperature profile with an exponentially decaying thermocline on top of a linear stratification. More...
 

Function/Subroutine Documentation

◆ circ_ridge()

real function neverworld_initialization::circ_ridge ( real, intent(in)  lon,
real, intent(in)  lat,
real, intent(in)  lon0,
real, intent(in)  lat0,
real, intent(in)  ring_radius,
real, intent(in)  ring_thickness,
real, intent(in)  ridge_height 
)
private

A circular ridge.

Parameters
[in]lonLongitude [degrees_E]
[in]latLatitude [degrees_N]
[in]lon0Longitude of center of ring [degrees_E]
[in]lat0Latitude of center of ring [degrees_N]
[in]ring_radiusRadius of ring [degrees]
[in]ring_thicknessRadial thickness of ring [degrees]
[in]ridge_heightRidge height as fraction of full depth [nondim]

Definition at line 223 of file Neverworld_initialization.F90.

223  real, intent(in) :: lon !< Longitude [degrees_E]
224  real, intent(in) :: lat !< Latitude [degrees_N]
225  real, intent(in) :: lon0 !< Longitude of center of ring [degrees_E]
226  real, intent(in) :: lat0 !< Latitude of center of ring [degrees_N]
227  real, intent(in) :: ring_radius !< Radius of ring [degrees]
228  real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees]
229  real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim]
230  real :: r
231 
232  r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point
233  r = abs( r - ring_radius) ! Pseudo-distance from a circle
234  r = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height
235  circ_ridge = 1. - r ! nondim depths (1-frac_ridge_height) .. 1

◆ cone()

real function neverworld_initialization::cone ( real, intent(in)  x,
real, intent(in)  x0,
real, intent(in)  L,
real, intent(in), optional  clip 
)
private

Returns the value of a triangular function centered at x=x0 with value 1 and linearly decreasing to 0 at x=x0+/-L, and 0 otherwise. If clip is present the top of the cone is cut off at "clip", which effectively defaults to 1.

Parameters
[in]xnon-dimensional coordinate [nondim]
[in]x0position of peak [nondim]
[in]lhalf-width of base of cone [nondim]
[in]clipclipping height of cone [nondim]

Definition at line 113 of file Neverworld_initialization.F90.

113  real, intent(in) :: x !< non-dimensional coordinate [nondim]
114  real, intent(in) :: x0 !< position of peak [nondim]
115  real, intent(in) :: l !< half-width of base of cone [nondim]
116  real, optional, intent(in) :: clip !< clipping height of cone [nondim]
117 
118  cone = max( 0., 1. - abs(x - x0) / l )
119  if (present(clip)) cone = min(clip, cone)

◆ cosbell()

real function neverworld_initialization::cosbell ( real, intent(in)  x,
real, intent(in)  L 
)
private

Returns the value of a cosine-bell function evaluated at x/L.

Parameters
[in]xnon-dimensional position
[in]lnon-dimensional width

Definition at line 89 of file Neverworld_initialization.F90.

89  real , intent(in) :: x !< non-dimensional position
90  real , intent(in) :: l !< non-dimensional width
91  real :: pi !< 3.1415926... calculated as 4*atan(1)
92 
93  pi = 4.0*atan(1.0)
94  cosbell = 0.5 * (1 + cos(pi*min(abs(x/l),1.0)))

◆ cstprof()

real function neverworld_initialization::cstprof ( real, intent(in)  x,
real, intent(in)  x0,
real, intent(in)  L,
real, intent(in)  lf,
real, intent(in)  bf,
real, intent(in)  sf,
real, intent(in)  sh 
)
private

Returns a "coastal" profile.

Parameters
[in]xnon-dimensional coordinate [nondim]
[in]x0position of peak [nondim]
[in]lwidth of profile [nondim]
[in]lffraction of width that is "land" [nondim]
[in]bffraction of width that is "beach" [nondim]
[in]sffraction of width that is "continental slope" [nondim]
[in]shdepth of shelf as fraction of full depth [nondim]

Definition at line 135 of file Neverworld_initialization.F90.

135  real, intent(in) :: x !< non-dimensional coordinate [nondim]
136  real, intent(in) :: x0 !< position of peak [nondim]
137  real, intent(in) :: l !< width of profile [nondim]
138  real, intent(in) :: lf !< fraction of width that is "land" [nondim]
139  real, intent(in) :: bf !< fraction of width that is "beach" [nondim]
140  real, intent(in) :: sf !< fraction of width that is "continental slope" [nondim]
141  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
142  real :: s
143 
144  s = max( 0., min( 1.,( x - x0 ) / l ) )
145  cstprof = sh * scurve(s-lf,0.,bf) + (1.-sh) * scurve(s - (1.-sf),0.,sf)

◆ dist_line_fixed_x()

real function neverworld_initialization::dist_line_fixed_x ( real, intent(in)  x,
real, intent(in)  y,
real, intent(in)  x0,
real, intent(in)  y0,
real, intent(in)  y1 
)
private

Distance between points x,y and a line segment (x0,y0) and (x0,y1).

Parameters
[in]xnon-dimensional x-coordinate [nondim]
[in]ynon-dimensional y-coordinate [nondim]
[in]x0x-position of line segment [nondim]
[in]y0y-position of line segment end[nondim]
[in]y1y-position of line segment end[nondim]

Definition at line 150 of file Neverworld_initialization.F90.

150  real, intent(in) :: x !< non-dimensional x-coordinate [nondim]
151  real, intent(in) :: y !< non-dimensional y-coordinate [nondim]
152  real, intent(in) :: x0 !< x-position of line segment [nondim]
153  real, intent(in) :: y0 !< y-position of line segment end[nondim]
154  real, intent(in) :: y1 !< y-position of line segment end[nondim]
155  real :: dx, yr, dy
156 
157  dx = x - x0
158  yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1
159  dy = y - yr ! =0 within y0<y<y1, =y0-y for y<y0, =y-y1 for y>y1
160  dist_line_fixed_x = sqrt( dx*dx + dy*dy )

◆ dist_line_fixed_y()

real function neverworld_initialization::dist_line_fixed_y ( real, intent(in)  x,
real, intent(in)  y,
real, intent(in)  x0,
real, intent(in)  x1,
real, intent(in)  y0 
)
private

Distance between points x,y and a line segment (x0,y0) and (x1,y0).

Parameters
[in]xnon-dimensional x-coordinate [nondim]
[in]ynon-dimensional y-coordinate [nondim]
[in]x0x-position of line segment end[nondim]
[in]x1x-position of line segment end[nondim]
[in]y0y-position of line segment [nondim]

Definition at line 165 of file Neverworld_initialization.F90.

165  real, intent(in) :: x !< non-dimensional x-coordinate [nondim]
166  real, intent(in) :: y !< non-dimensional y-coordinate [nondim]
167  real, intent(in) :: x0 !< x-position of line segment end[nondim]
168  real, intent(in) :: x1 !< x-position of line segment end[nondim]
169  real, intent(in) :: y0 !< y-position of line segment [nondim]
170  real :: dx, yr, dy
171 
172  dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1)

◆ ew_coast()

real function neverworld_initialization::ew_coast ( real, intent(in)  lon,
real, intent(in)  lat,
real, intent(in)  lon0,
real, intent(in)  lon1,
real, intent(in)  lat0,
real, intent(in)  dlat,
real, intent(in)  sh 
)
private

A "coast profile" applied in an E-W line from lon0,lat0 to lon1,lat0.

Parameters
[in]lonLongitude [degrees_E]
[in]latLatitude [degrees_N]
[in]lon0Longitude of coast end [degrees_E]
[in]lon1Longitude of coast end [degrees_E]
[in]lat0Latitude of coast [degrees_N]
[in]dlat"Radius" of coast profile [degrees]
[in]shdepth of shelf as fraction of full depth [nondim]

Definition at line 192 of file Neverworld_initialization.F90.

192  real, intent(in) :: lon !< Longitude [degrees_E]
193  real, intent(in) :: lat !< Latitude [degrees_N]
194  real, intent(in) :: lon0 !< Longitude of coast end [degrees_E]
195  real, intent(in) :: lon1 !< Longitude of coast end [degrees_E]
196  real, intent(in) :: lat0 !< Latitude of coast [degrees_N]
197  real, intent(in) :: dlat !< "Radius" of coast profile [degrees]
198  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
199  real :: r
200 
201  r = dist_line_fixed_y( lon, lat, lon0, lon1, lat0 )
202  ew_coast = cstprof(r, 0., dlat, 0.125, 0.125, 0.5, sh)

◆ neverworld_initialize_thickness()

subroutine, public neverworld_initialization::neverworld_initialize_thickness ( real, dimension(szi_(g),szj_(g),szk_(gv)), intent(out)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(eos_type), pointer  eqn_of_state,
real, intent(in)  P_ref 
)

This subroutine initializes layer thicknesses for the Neverworld test case, by finding the depths of interfaces in a specified latitude-dependent temperature profile with an exponentially decaying thermocline on top of a linear stratification.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[out]hThe thickness that is being initialized [H ~> m or kg m-2].
[in]param_fileA structure indicating the open file to parse for model parameter values.
eqn_of_stateEquation of state structure
[in]p_refThe coordinate-density reference pressure [R L2 T-2 ~> Pa].

Definition at line 243 of file Neverworld_initialization.F90.

243  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
244  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
245  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
246  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being
247  !! initialized [H ~> m or kg m-2].
248  type(param_file_type), intent(in) :: param_file !< A structure indicating the open
249  !! file to parse for model
250  !! parameter values.
251  type(eos_type), pointer :: eqn_of_state !< Equation of state structure
252  real, intent(in) :: p_ref !< The coordinate-density
253  !! reference pressure [R L2 T-2 ~> Pa].
254  ! Local variables
255  real :: e0(szk_(g)+1) ! The resting interface heights, in depth units [Z ~> m],
256  ! usually negative because it is positive upward.
257  real, dimension(SZK_(G)) :: h_profile ! Vector of initial thickness profile [Z ~> m].
258  real :: e_interface ! Current interface position [Z ~> m].
259  real :: x,y,r1,r2 ! x,y and radial coordinates for computation of initial pert.
260  real :: pert_amp ! Amplitude of perturbations measured in Angstrom_H
261  real :: h_noise ! Amplitude of noise to scale h by
262  real :: noise ! Noise
263  type(randomnumberstream) :: rns ! Random numbers for stochastic tidal parameterization
264  character(len=40) :: mdl = "Neverworld_initialize_thickness" ! This subroutine's name.
265  integer :: i, j, k, k1, is, ie, js, je, nz, itt
266 
267  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
268 
269  call mom_mesg(" Neverworld_initialization.F90, Neverworld_initialize_thickness: setting thickness", 5)
270  call get_param(param_file, mdl, "INIT_THICKNESS_PROFILE", h_profile, &
271  "Profile of initial layer thicknesses.", units="m", scale=us%m_to_Z, &
272  fail_if_missing=.true.)
273  call get_param(param_file, mdl, "NL_THICKNESS_PERT_AMP", pert_amp, &
274  "Amplitude of finite scale perturbations as fraction of depth.", &
275  units="nondim", default=0.)
276  call get_param(param_file, mdl, "NL_THICKNESS_NOISE_AMP", h_noise, &
277  "Amplitude of noise to scale layer by.", units="nondim", default=0.)
278 
279  ! e0 is the notional position of interfaces
280  e0(1) = 0. ! The surface
281  do k=1,nz
282  e0(k+1) = e0(k) - h_profile(k)
283  enddo
284 
285  do j=js,je ; do i=is,ie
286  e_interface = -g%bathyT(i,j)
287  do k=nz,2,-1
288  h(i,j,k) = gv%Z_to_H * (e0(k) - e_interface) ! Nominal thickness
289  x=(g%geoLonT(i,j)-g%west_lon)/g%len_lon
290  y=(g%geoLatT(i,j)-g%south_lat)/g%len_lat
291  r1=sqrt((x-0.7)**2+(y-0.2)**2)
292  r2=sqrt((x-0.3)**2+(y-0.25)**2)
293  h(i,j,k) = h(i,j,k) + pert_amp * (e0(k) - e0(nz+1)) * gv%Z_to_H * &
294  (spike(r1,0.15)-spike(r2,0.15)) ! Prescribed perturbation
295  if (h_noise /= 0.) then
296  rns = initializerandomnumberstream( int( 4096*(x + (y+1.)) ) )
297  call getrandomnumbers(rns, noise) ! x will be in (0,1)
298  noise = h_noise * 2. * ( noise - 0.5 ) ! range -h_noise to h_noise
299  h(i,j,k) = ( 1. + noise ) * h(i,j,k)
300  endif
301  h(i,j,k) = max( gv%Angstrom_H, h(i,j,k) ) ! Limit to non-negative
302  e_interface = e_interface + gv%H_to_Z * h(i,j,k) ! Actual position of upper interface
303  enddo
304  h(i,j,1) = gv%Z_to_H * (e0(1) - e_interface) ! Nominal thickness
305  h(i,j,1) = max( gv%Angstrom_H, h(i,j,1) ) ! Limit to non-negative
306  enddo ; enddo
307 

◆ neverworld_initialize_topography()

subroutine, public neverworld_initialization::neverworld_initialize_topography ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)  D,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
real, intent(in)  max_depth 
)

This subroutine sets up the Neverworld test case topography.

Parameters
[in]gThe dynamic horizontal grid type
[out]dOcean bottom depth in the units of depth_max
[in]param_fileParameter file structure
[in]max_depthMaximum ocean depth in arbitrary units

Definition at line 36 of file Neverworld_initialization.F90.

36  type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
37  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
38  intent(out) :: d !< Ocean bottom depth in the units of depth_max
39  type(param_file_type), intent(in) :: param_file !< Parameter file structure
40  real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units
41 
42  ! Local variables
43  real :: pi ! 3.1415926... calculated as 4*atan(1)
44  real :: d0 ! A constant to make the maximum !
45  ! basin depth MAXIMUM_DEPTH. !
46  real :: x, y
47  ! This include declares and sets the variable "version".
48 # include "version_variable.h"
49  character(len=40) :: mdl = "Neverworld_initialize_topography" ! This subroutine's name.
50  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
51  real :: nl_roughness_amp, nl_top_amp
52  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
53  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
54 
55  call mom_mesg(" Neverworld_initialization.F90, Neverworld_initialize_topography: setting topography", 5)
56 
57  call log_version(param_file, mdl, version, "")
58  call get_param(param_file, mdl, "NL_ROUGHNESS_AMP", nl_roughness_amp, &
59  "Amplitude of wavy signal in bathymetry.", default=0.05)
60  call get_param(param_file, mdl, "NL_CONTINENT_AMP", nl_top_amp, &
61  "Scale factor for topography - 0.0 for no continents.", default=1.0)
62 
63  pi = 4.0*atan(1.0)
64 
65 ! Calculate the depth of the bottom.
66  do j=js,je ; do i=is,ie
67  x = (g%geoLonT(i,j)-g%west_lon) / g%len_lon
68  y =( g%geoLatT(i,j)-g%south_lat) / g%len_lat
69 ! This sets topography that has a reentrant channel to the south.
70  d(i,j) = 1.0 - 1.1 * spike(y-1,0.12) - 1.1 * spike(y,0.12) - & !< The great northern wall and Antarctica
71  nl_top_amp*( &
72  (1.2 * spike(x,0.2) + 1.2 * spike(x-1.0,0.2)) * spike(min(0.0,y-0.3),0.2) & !< South America
73  + 1.2 * spike(x-0.5,0.2) * spike(min(0.0,y-0.55),0.2) & !< Africa
74  + 1.2 * (spike(x,0.12) + spike(x-1,0.12)) * spike(max(0.0,y-0.06),0.12) & !< Antarctic Peninsula
75  + 0.1 * (cosbell(x,0.1) + cosbell(x-1,0.1)) & !< Drake Passage ridge
76  + 0.5 * cosbell(x-0.16,0.05) * (cosbell(y-0.18,0.13)**0.4) & !< Scotia Arc East
77  + 0.4 * (cosbell(x-0.09,0.08)**0.4) * cosbell(y-0.26,0.05) & !< Scotia Arc North
78  + 0.4 * (cosbell(x-0.08,0.08)**0.4) * cosbell(y-0.1,0.05)) & !< Scotia Arc South
79  - nl_roughness_amp * cos(14*pi*x) * sin(14*pi*y) & !< roughness
80  - nl_roughness_amp * cos(20*pi*x) * cos(20*pi*y) !< roughness
81  if (d(i,j) < 0.0) d(i,j) = 0.0
82  d(i,j) = d(i,j) * max_depth
83  enddo ; enddo
84 

◆ ns_coast()

real function neverworld_initialization::ns_coast ( real, intent(in)  lon,
real, intent(in)  lat,
real, intent(in)  lon0,
real, intent(in)  lat0,
real, intent(in)  lat1,
real, intent(in)  dlon,
real, intent(in)  sh 
)
private

A "coast profile" applied in an N-S line from lon0,lat0 to lon0,lat1.

Parameters
[in]lonLongitude [degrees_E]
[in]latLatitude [degrees_N]
[in]lon0Longitude of coast [degrees_E]
[in]lat0Latitude of coast end [degrees_N]
[in]lat1Latitude of coast end [degrees_N]
[in]dlon"Radius" of coast profile [degrees]
[in]shdepth of shelf as fraction of full depth [nondim]

Definition at line 177 of file Neverworld_initialization.F90.

177  real, intent(in) :: lon !< Longitude [degrees_E]
178  real, intent(in) :: lat !< Latitude [degrees_N]
179  real, intent(in) :: lon0 !< Longitude of coast [degrees_E]
180  real, intent(in) :: lat0 !< Latitude of coast end [degrees_N]
181  real, intent(in) :: lat1 !< Latitude of coast end [degrees_N]
182  real, intent(in) :: dlon !< "Radius" of coast profile [degrees]
183  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
184  real :: r
185 
186  r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 )
187  ns_coast = cstprof(r, 0., dlon, 0.125, 0.125, 0.5, sh)

◆ ns_ridge()

real function neverworld_initialization::ns_ridge ( real, intent(in)  lon,
real, intent(in)  lat,
real, intent(in)  lon0,
real, intent(in)  lat0,
real, intent(in)  lat1,
real, intent(in)  dlon,
real, intent(in)  rh 
)
private

A NS ridge.

Parameters
[in]lonLongitude [degrees_E]
[in]latLatitude [degrees_N]
[in]lon0Longitude of ridge center [degrees_E]
[in]lat0Latitude of ridge end [degrees_N]
[in]lat1Latitude of ridge end [degrees_N]
[in]dlon"Radius" of ridge profile [degrees]
[in]rhdepth of ridge as fraction of full depth [nondim]

Definition at line 207 of file Neverworld_initialization.F90.

207  real, intent(in) :: lon !< Longitude [degrees_E]
208  real, intent(in) :: lat !< Latitude [degrees_N]
209  real, intent(in) :: lon0 !< Longitude of ridge center [degrees_E]
210  real, intent(in) :: lat0 !< Latitude of ridge end [degrees_N]
211  real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N]
212  real, intent(in) :: dlon !< "Radius" of ridge profile [degrees]
213  real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim]
214  real :: r
215 
216  r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 )
217  ns_ridge = 1. - rh * cone(r, 0., dlon)

◆ scurve()

real function neverworld_initialization::scurve ( real, intent(in)  x,
real, intent(in)  x0,
real, intent(in)  L 
)
private

Returns an s-curve s(x) s.t. s(x0)<=0, s(x0+L)>=1 and cubic in between.

Parameters
[in]xnon-dimensional coordinate [nondim]
[in]x0position of peak [nondim]
[in]lhalf-width of base of cone [nondim]

Definition at line 124 of file Neverworld_initialization.F90.

124  real, intent(in) :: x !< non-dimensional coordinate [nondim]
125  real, intent(in) :: x0 !< position of peak [nondim]
126  real, intent(in) :: l !< half-width of base of cone [nondim]
127  real :: s
128 
129  s = max( 0., min( 1.,( x - x0 ) / l ) )
130  scurve = ( 3. - 2.*s ) * ( s * s )

◆ spike()

real function neverworld_initialization::spike ( real, intent(in)  x,
real, intent(in)  L 
)
private

Returns the value of a sin-spike function evaluated at x/L.

Parameters
[in]xnon-dimensional position
[in]lnon-dimensional width

Definition at line 99 of file Neverworld_initialization.F90.

99 
100  real , intent(in) :: x !< non-dimensional position
101  real , intent(in) :: l !< non-dimensional width
102  real :: pi !< 3.1415926... calculated as 4*atan(1)
103 
104  pi = 4.0*atan(1.0)
105  spike = (1 - sin(pi*min(abs(x/l),0.5)))