MOM6
Functions/Subroutines | Variables
basin_builder Module Reference

Detailed Description

An idealized topography building system.

Functions/Subroutines

subroutine, public basin_builder_topography (D, G, param_file, max_depth)
 Constructs idealized topography from simple functions. 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, lonC, lat0, lat1, dlon, sh)
 A "coast profile" applied in an N-S line from lonC,lat0 to lonC,lat1. More...
 
real function ew_coast (lon, lat, latC, lon0, lon1, dlat, sh)
 A "coast profile" applied in an E-W line from lon0,latC to lon1,latC. More...
 
real function ns_conic_ridge (lon, lat, lonC, lat0, lat1, dlon, rh)
 A NS ridge with a cone profile. More...
 
real function ns_scurve_ridge (lon, lat, lonC, lat0, lat1, dlon, rh)
 A NS ridge with an scurve profile. More...
 
real function circ_conic_ridge (lon, lat, lon0, lat0, ring_radius, ring_thickness, ridge_height)
 A circular ridge with cutoff conic profile. More...
 
real function circ_scurve_ridge (lon, lat, lon0, lat0, ring_radius, ring_thickness, ridge_height)
 A circular ridge with cutoff scurve profile. More...
 

Variables

character(len=40) mdl = "basin_builder"
 This module's name.
 

Function/Subroutine Documentation

◆ basin_builder_topography()

subroutine, public basin_builder::basin_builder_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 
)

Constructs idealized topography from simple functions.

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 28 of file basin_builder.F90.

28  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
29  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
30  intent(out) :: D !< Ocean bottom depth in the units of depth_max
31  type(param_file_type), intent(in) :: param_file !< Parameter file structure
32  real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units
33  ! Local variables
34  character(len=17) :: pname1, pname2 ! For construction of parameter names
35  character(len=20) :: funcs ! Basin build function
36  real, dimension(20) :: pars ! Parameters for each function
37  real :: lon ! Longitude [degrees_E}
38  real :: lat ! Latitude [degrees_N]
39  integer :: i, j, n, n_funcs
40 
41  call mom_mesg(" basin_builder.F90, basin_builder_topography: setting topography", 5)
42  call log_version(param_file, mdl, version, "")
43 
44  do j=g%jsc,g%jec ; do i=g%isc,g%iec
45  d(i,j) = 1.0
46  enddo ; enddo
47 
48  call get_param(param_file, mdl, "BBUILDER_N", n_funcs, &
49  "Number of pieces of topography to use.", fail_if_missing=.true.)
50 
51  do n=1,n_funcs
52  write( pname1, "('BBUILDER_',i3.3,'_FUNC')" ) n
53  write( pname2, "('BBUILDER_',i3.3,'_PARS')" ) n
54  call get_param(param_file, mdl, pname1, funcs, &
55  "The basin builder function to apply with parameters "//&
56  trim(pname2)//". Choices are: NS_COAST, EW_COAST, "//&
57  "CIRC_CONIC_RIDGE, NS_CONIC_RIDGE, CIRC_SCURVE_RIDGE, "//&
58  "NS_SCURVE_RIDGE.", &
59  fail_if_missing=.true.)
60  pars(:) = 0.
61  if (trim(lowercase(funcs)) == 'ns_coast') then
62  call get_param(param_file, mdl, pname2, pars(1:5), &
63  "NS_COAST parameters: longitude, starting latitude, "//&
64  "ending latitude, footprint radius, shelf depth.", &
65  units="degrees_E,degrees_N,degrees_N,degrees,m", &
66  fail_if_missing=.true.)
67  pars(5) = pars(5) / max_depth
68  do j=g%jsc,g%jec ; do i=g%isc,g%iec
69  lon = g%geoLonT(i,j)
70  lat = g%geoLatT(i,j)
71  d(i,j) = min( d(i,j), ns_coast(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
72  enddo ; enddo
73  elseif (trim(lowercase(funcs)) == 'ns_conic_ridge') then
74  call get_param(param_file, mdl, pname2, pars(1:5), &
75  "NS_CONIC_RIDGE parameters: longitude, starting latitude, "//&
76  "ending latitude, footprint radius, ridge height.", &
77  units="degrees_E,degrees_N,degrees_N,degrees,m", &
78  fail_if_missing=.true.)
79  pars(5) = pars(5) / max_depth
80  do j=g%jsc,g%jec ; do i=g%isc,g%iec
81  lon = g%geoLonT(i,j)
82  lat = g%geoLatT(i,j)
83  d(i,j) = min( d(i,j), ns_conic_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
84  enddo ; enddo
85  elseif (trim(lowercase(funcs)) == 'ns_scurve_ridge') then
86  call get_param(param_file, mdl, pname2, pars(1:5), &
87  "NS_SCURVE_RIDGE parameters: longitude, starting latitude, "//&
88  "ending latitude, footprint radius, ridge height.", &
89  units="degrees_E,degrees_N,degrees_N,degrees,m", &
90  fail_if_missing=.true.)
91  pars(5) = pars(5) / max_depth
92  do j=g%jsc,g%jec ; do i=g%isc,g%iec
93  lon = g%geoLonT(i,j)
94  lat = g%geoLatT(i,j)
95  d(i,j) = min( d(i,j), ns_scurve_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
96  enddo ; enddo
97  elseif (trim(lowercase(funcs)) == 'ew_coast') then
98  call get_param(param_file, mdl, pname2, pars(1:5), &
99  "EW_COAST parameters: latitude, starting longitude, "//&
100  "ending longitude, footprint radius, shelf depth.", &
101  units="degrees_N,degrees_E,degrees_E,degrees,m", &
102  fail_if_missing=.true.)
103  pars(5) = pars(5) / max_depth
104  do j=g%jsc,g%jec ; do i=g%isc,g%iec
105  lon = g%geoLonT(i,j)
106  lat = g%geoLatT(i,j)
107  d(i,j) = min( d(i,j), ew_coast(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
108  enddo ; enddo
109  elseif (trim(lowercase(funcs)) == 'circ_conic_ridge') then
110  call get_param(param_file, mdl, pname2, pars(1:5), &
111  "CIRC_CONIC_RIDGE parameters: center longitude, center latitude, "//&
112  "ring radius, footprint radius, ridge height.", &
113  units="degrees_E,degrees_N,degrees,degrees,m", &
114  fail_if_missing=.true.)
115  pars(5) = pars(5) / max_depth
116  do j=g%jsc,g%jec ; do i=g%isc,g%iec
117  lon = g%geoLonT(i,j)
118  lat = g%geoLatT(i,j)
119  d(i,j) = min( d(i,j), circ_conic_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
120  enddo ; enddo
121  elseif (trim(lowercase(funcs)) == 'circ_scurve_ridge') then
122  call get_param(param_file, mdl, pname2, pars(1:5), &
123  "CIRC_SCURVe_RIDGE parameters: center longitude, center latitude, "//&
124  "ring radius, footprint radius, ridge height.", &
125  units="degrees_E,degrees_N,degrees,degrees,m", &
126  fail_if_missing=.true.)
127  pars(5) = pars(5) / max_depth
128  do j=g%jsc,g%jec ; do i=g%isc,g%iec
129  lon = g%geoLonT(i,j)
130  lat = g%geoLatT(i,j)
131  d(i,j) = min( d(i,j), circ_scurve_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
132  enddo ; enddo
133  else
134  call mom_error(fatal, "basin_builder.F90, basin_builer_topography:\n"//&
135  "Unrecognized function "//trim(funcs))
136  endif
137 
138  enddo ! n
139 
140  do j=g%jsc,g%jec ; do i=g%isc,g%iec
141  ! Dimensionalize by scaling 1 to max_depth
142  d(i,j) = d(i,j) * max_depth
143  enddo ; enddo
144 

◆ circ_conic_ridge()

real function basin_builder::circ_conic_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 with cutoff conic profile.

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 276 of file basin_builder.F90.

276  real, intent(in) :: lon !< Longitude [degrees_E]
277  real, intent(in) :: lat !< Latitude [degrees_N]
278  real, intent(in) :: lon0 !< Longitude of center of ring [degrees_E]
279  real, intent(in) :: lat0 !< Latitude of center of ring [degrees_N]
280  real, intent(in) :: ring_radius !< Radius of ring [degrees]
281  real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees]
282  real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim]
283  real :: r
284 
285  r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point
286  r = abs( r - ring_radius) ! Pseudo-distance from a circle
287  r = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height
288  circ_conic_ridge = 1. - r ! nondim depths (1-frac_ridge_height) .. 1

◆ circ_scurve_ridge()

real function basin_builder::circ_scurve_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 with cutoff scurve profile.

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 293 of file basin_builder.F90.

293  real, intent(in) :: lon !< Longitude [degrees_E]
294  real, intent(in) :: lat !< Latitude [degrees_N]
295  real, intent(in) :: lon0 !< Longitude of center of ring [degrees_E]
296  real, intent(in) :: lat0 !< Latitude of center of ring [degrees_N]
297  real, intent(in) :: ring_radius !< Radius of ring [degrees]
298  real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees]
299  real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim]
300  real :: r
301 
302  r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point
303  r = abs( r - ring_radius) ! Pseudo-distance from a circle
304  r = 1. - scurve(r, 0., ring_thickness) ! 0 .. 1
305  r = r * ridge_height ! 0 .. frac_ridge_height
306  circ_scurve_ridge = 1. - r ! nondim depths (1-frac_ridge_height) .. 1

◆ cone()

real function basin_builder::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 152 of file basin_builder.F90.

152  real, intent(in) :: x !< non-dimensional coordinate [nondim]
153  real, intent(in) :: x0 !< position of peak [nondim]
154  real, intent(in) :: L !< half-width of base of cone [nondim]
155  real, optional, intent(in) :: clip !< clipping height of cone [nondim]
156 
157  cone = max( 0., 1. - abs(x - x0) / l )
158  if (present(clip)) cone = min(clip, cone)

◆ cstprof()

real function basin_builder::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 174 of file basin_builder.F90.

174  real, intent(in) :: x !< non-dimensional coordinate [nondim]
175  real, intent(in) :: x0 !< position of peak [nondim]
176  real, intent(in) :: L !< width of profile [nondim]
177  real, intent(in) :: lf !< fraction of width that is "land" [nondim]
178  real, intent(in) :: bf !< fraction of width that is "beach" [nondim]
179  real, intent(in) :: sf !< fraction of width that is "continental slope" [nondim]
180  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
181  real :: s
182 
183  s = max( 0., min( 1.,( x - x0 ) / l ) )
184  cstprof = sh * scurve(s-lf,0.,bf) + (1.-sh) * scurve(s - (1.-sf),0.,sf)

◆ dist_line_fixed_x()

real function basin_builder::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 189 of file basin_builder.F90.

189  real, intent(in) :: x !< non-dimensional x-coordinate [nondim]
190  real, intent(in) :: y !< non-dimensional y-coordinate [nondim]
191  real, intent(in) :: x0 !< x-position of line segment [nondim]
192  real, intent(in) :: y0 !< y-position of line segment end[nondim]
193  real, intent(in) :: y1 !< y-position of line segment end[nondim]
194  real :: dx, yr, dy
195 
196  dx = x - x0
197  yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1
198  dy = y - yr ! =0 within y0<y<y1, =y0-y for y<y0, =y-y1 for y>y1
199  dist_line_fixed_x = sqrt( dx*dx + dy*dy )

◆ dist_line_fixed_y()

real function basin_builder::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 204 of file basin_builder.F90.

204  real, intent(in) :: x !< non-dimensional x-coordinate [nondim]
205  real, intent(in) :: y !< non-dimensional y-coordinate [nondim]
206  real, intent(in) :: x0 !< x-position of line segment end[nondim]
207  real, intent(in) :: x1 !< x-position of line segment end[nondim]
208  real, intent(in) :: y0 !< y-position of line segment [nondim]
209  real :: dx, yr, dy
210 
211  dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1)

◆ ew_coast()

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

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

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

Definition at line 231 of file basin_builder.F90.

231  real, intent(in) :: lon !< Longitude [degrees_E]
232  real, intent(in) :: lat !< Latitude [degrees_N]
233  real, intent(in) :: latC !< Latitude of coast [degrees_N]
234  real, intent(in) :: lon0 !< Longitude of coast end [degrees_E]
235  real, intent(in) :: lon1 !< Longitude of coast end [degrees_E]
236  real, intent(in) :: dlat !< "Radius" of coast profile [degrees]
237  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
238  real :: r
239 
240  r = dist_line_fixed_y( lon, lat, lon0, lon1, latc )
241  ew_coast = cstprof(r, 0., dlat, 0.125, 0.125, 0.5, sh)

◆ ns_coast()

real function basin_builder::ns_coast ( real, intent(in)  lon,
real, intent(in)  lat,
real, intent(in)  lonC,
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 lonC,lat0 to lonC,lat1.

Parameters
[in]lonLongitude [degrees_E]
[in]latLatitude [degrees_N]
[in]loncLongitude 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 216 of file basin_builder.F90.

216  real, intent(in) :: lon !< Longitude [degrees_E]
217  real, intent(in) :: lat !< Latitude [degrees_N]
218  real, intent(in) :: lonC !< Longitude of coast [degrees_E]
219  real, intent(in) :: lat0 !< Latitude of coast end [degrees_N]
220  real, intent(in) :: lat1 !< Latitude of coast end [degrees_N]
221  real, intent(in) :: dlon !< "Radius" of coast profile [degrees]
222  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
223  real :: r
224 
225  r = dist_line_fixed_x( lon, lat, lonc, lat0, lat1 )
226  ns_coast = cstprof(r, 0., dlon, 0.125, 0.125, 0.5, sh)

◆ ns_conic_ridge()

real function basin_builder::ns_conic_ridge ( real, intent(in)  lon,
real, intent(in)  lat,
real, intent(in)  lonC,
real, intent(in)  lat0,
real, intent(in)  lat1,
real, intent(in)  dlon,
real, intent(in)  rh 
)
private

A NS ridge with a cone profile.

Parameters
[in]lonLongitude [degrees_E]
[in]latLatitude [degrees_N]
[in]loncLongitude 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 246 of file basin_builder.F90.

246  real, intent(in) :: lon !< Longitude [degrees_E]
247  real, intent(in) :: lat !< Latitude [degrees_N]
248  real, intent(in) :: lonC !< Longitude of ridge center [degrees_E]
249  real, intent(in) :: lat0 !< Latitude of ridge end [degrees_N]
250  real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N]
251  real, intent(in) :: dlon !< "Radius" of ridge profile [degrees]
252  real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim]
253  real :: r
254 
255  r = dist_line_fixed_x( lon, lat, lonc, lat0, lat1 )
256  ns_conic_ridge = 1. - rh * cone(r, 0., dlon)

◆ ns_scurve_ridge()

real function basin_builder::ns_scurve_ridge ( real, intent(in)  lon,
real, intent(in)  lat,
real, intent(in)  lonC,
real, intent(in)  lat0,
real, intent(in)  lat1,
real, intent(in)  dlon,
real, intent(in)  rh 
)
private

A NS ridge with an scurve profile.

Parameters
[in]lonLongitude [degrees_E]
[in]latLatitude [degrees_N]
[in]loncLongitude 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 261 of file basin_builder.F90.

261  real, intent(in) :: lon !< Longitude [degrees_E]
262  real, intent(in) :: lat !< Latitude [degrees_N]
263  real, intent(in) :: lonC !< Longitude of ridge center [degrees_E]
264  real, intent(in) :: lat0 !< Latitude of ridge end [degrees_N]
265  real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N]
266  real, intent(in) :: dlon !< "Radius" of ridge profile [degrees]
267  real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim]
268  real :: r
269 
270  r = dist_line_fixed_x( lon, lat, lonc, lat0, lat1 )
271  ns_scurve_ridge = 1. - rh * (1. - scurve(r, 0., dlon) )

◆ scurve()

real function basin_builder::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 163 of file basin_builder.F90.

163  real, intent(in) :: x !< non-dimensional coordinate [nondim]
164  real, intent(in) :: x0 !< position of peak [nondim]
165  real, intent(in) :: L !< half-width of base of cone [nondim]
166  real :: s
167 
168  s = max( 0., min( 1.,( x - x0 ) / l ) )
169  scurve = ( 3. - 2.*s ) * ( s * s )