1
! Copyright (C) 2007 Imperial College London and others.
3
! Please see the AUTHORS file in the main source directory for a full list
4
! of copyright holders.
7
! Applied Modelling and Computation Group
8
! Department of Earth Science and Engineering
9
! Imperial College London
11
! amcgsoftware@imperial.ac.uk
13
! This library is free software; you can redistribute it and/or
14
! modify it under the terms of the GNU Lesser General Public
15
! License as published by the Free Software Foundation,
16
! version 2.1 of the License.
18
! This library is distributed in the hope that it will be useful,
19
! but WITHOUT ANY WARRANTY; without even the implied warranty of
20
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21
! Lesser General Public License for more details.
23
! You should have received a copy of the GNU Lesser General Public
24
! License along with this library; if not, write to the Free Software
25
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
32
module vertical_extrapolation_module
33
!!< Module containing routines for vertical extrapolation on
34
!!< fully unstructured 3D meshes. Also contains routines for updating
35
!!< distance to top and bottom fields.
41
use transform_elements
42
use boundary_conditions
44
use global_parameters, only: real_4, real_8
46
use dynamic_bin_sort_module
48
use integer_set_module
52
interface VerticalExtrapolation
53
module procedure VerticalExtrapolationScalar, &
54
VerticalExtrapolationMultiple, VerticalExtrapolationVector
55
end interface VerticalExtrapolation
57
interface vertical_integration
58
module procedure vertical_integration_scalar, &
59
vertical_integration_multiple, vertical_integration_vector
60
end interface vertical_integration
62
!! element i is considered to be above j iff the inner product of the face
63
!! normal point from i to j and the gravity normal is bigger than this
64
real, parameter:: VERTICAL_INTEGRATION_EPS=1.0e-8
66
! the two hemispheres are both projected to the unit circle, but the projected southern
67
! hemisphere will have its origin shifted. The projected hemispheres are slightly bigger
68
! than the unit circle, as a bit of overlap with the opposite hemisphere is included around
69
! the equator. So this shift needs to be big enough to not make both projected hemispheres
71
real, parameter:: HEMISPHERE_SHIFT=10.0
75
public CalculateTopBottomDistance
76
public VerticalExtrapolation, vertical_integration
77
public vertical_element_ordering
78
public VerticalProlongationOperator
79
public vertical_extrapolation_module_check_options
83
subroutine UpdateDistanceField(state, name, vertical_coordinate)
84
! This sub calculates the vertical distance to the free surface
85
! and bottom of the ocean to all nodes. The results are stored
86
! in the 'DistanceToBottom/FreeSurface' fields from state.
87
type(state_type), intent(inout):: state
88
character(len=*), intent(in):: name
89
type(scalar_field), intent(in):: vertical_coordinate
92
type(vector_field), pointer:: positions, vertical_normal
93
type(scalar_field), pointer:: distance
95
integer, pointer, dimension(:):: surface_element_list
97
! the distance field to compute:
98
distance => extract_scalar_field(state, name)
99
! its first boundary condition is on the related top or bottom mesh
100
call get_boundary_condition(distance, 1, &
101
surface_element_list=surface_element_list)
102
positions => extract_vector_field(state, "Coordinate")
103
vertical_normal => extract_vector_field(state, "GravityDirection")
105
! in each node of the mesh, set "distance" to the vertical coordinate
106
! of this node projected to the above/below surface mesh
107
call VerticalExtrapolation(vertical_coordinate, distance, positions, &
108
vertical_normal, surface_element_list, surface_name=name)
110
! the distance is then calculated by subtracting its own vertical coordinate
111
call addto(distance, vertical_coordinate, scale=-1.0)
113
if (name=="DistanceToBottom") then
114
! make distance to bottom positive
115
call scale(distance, -1.0)
118
end subroutine UpdateDistanceField
120
subroutine CalculateTopBottomDistance(state)
121
!! This sub calculates the vertical distance to the free surface
122
!! and bottom of the ocean to all nodes. The results are stored
123
!! in the 'DistanceToBottom/Top' fields from state.
124
type(state_type), intent(inout):: state
126
type(mesh_type), pointer:: xmesh
127
type(scalar_field):: vertical_coordinate
129
xmesh => extract_mesh(state, "CoordinateMesh")
130
call allocate(vertical_coordinate, xmesh, "VerticalCoordinate")
131
call calculate_vertical_coordinate(state, vertical_coordinate)
132
call UpdateDistanceField(state, "DistanceToTop", vertical_coordinate)
133
call UpdateDistanceField(state, "DistanceToBottom", vertical_coordinate)
134
call deallocate(vertical_coordinate)
136
end subroutine CalculateTopBottomDistance
138
subroutine calculate_vertical_coordinate(state, vertical_coordinate)
139
!! Computes a vertical coordinate, i.e. a scalar field such that
140
!! for each 2 nodes above each other, the difference of the field
141
!! in these nodes gives the distance between them.
142
type(state_type), intent(inout):: state
143
type(scalar_field), intent(inout):: vertical_coordinate
145
type(vector_field), pointer:: positions, gravity_normal
146
type(scalar_field):: positions_magnitude
148
positions => extract_vector_field(state, "Coordinate")
149
if(have_option('/geometry/spherical_earth')) then
150
! use the radius as vertical coordinate
151
! that is, the l2-norm of the coordinate field
152
positions_magnitude=magnitude(positions)
153
call set(vertical_coordinate, positions_magnitude)
154
call deallocate(positions_magnitude)
156
gravity_normal => extract_vector_field(state, "GravityDirection")
157
assert(gravity_normal%field_type==FIELD_TYPE_CONSTANT)
158
call inner_product(vertical_coordinate, gravity_normal, positions)
159
! gravity points down, we want a vertical coordinate that increases upward
160
call scale(vertical_coordinate, -1.0)
163
end subroutine calculate_vertical_coordinate
165
subroutine VerticalExtrapolationScalar(from_field, to_field, &
166
positions, vertical_normal, surface_element_list, surface_name)
167
!!< This sub extrapolates the values on a horizontal 2D surface
168
!!< in the vertical direction to 3D fields
169
!! The from_fields should be 3D fields of which only the values on the
170
!! 2D horizontal surface are used.
171
type(scalar_field), intent(in):: from_field
172
!! Resulting extrapolated field. May be the same field or a field on
173
!! a different mesh (different degree).
174
type(scalar_field), intent(inout):: to_field
175
!! positions, and upward normal vector on the whole domain
176
type(vector_field), target, intent(inout):: positions
177
type(vector_field), target, intent(in):: vertical_normal
178
!! the surface elements (faces numbers) that make up the surface
179
integer, dimension(:), intent(in):: surface_element_list
180
!! If provided the projected surface mesh onto horizontal coordinates
181
!! and its associated rtree/pickers are cached under this name and
182
!! attached to the 'positions'. In this case when called again with
183
!! the same 'positions' and the same surface_name,
184
!! the same surface_element_list should again be provided.
185
character(len=*), optional, intent(in):: surface_name
187
type(scalar_field), dimension(1):: to_fields
189
to_fields=(/ to_field /)
190
call VerticalExtrapolationMultiple( (/ from_field /) , to_fields, &
191
positions, vertical_normal, surface_element_list, &
192
surface_name=surface_name)
194
end subroutine VerticalExtrapolationScalar
196
subroutine VerticalExtrapolationVector(from_field, to_field, &
197
positions, vertical_normal, surface_element_list, surface_name)
198
!!< This sub extrapolates the values on a horizontal 2D surface
199
!!< in the vertical direction to 3D fields
200
!! The from_fields should be 3D fields of which only the values on the
201
!! 2D horizontal surface are used.
202
type(vector_field), intent(in):: from_field
203
!! Resulting extrapolated field. May be the same field or a field on
204
!! a different mesh (different degree).
205
type(vector_field), intent(inout):: to_field
206
!! positions, and upward normal vector on the whole domain
207
type(vector_field), target, intent(inout):: positions
208
type(vector_field), target, intent(in):: vertical_normal
209
!! the surface elements (faces numbers) that make up the surface
210
integer, dimension(:), intent(in):: surface_element_list
211
!! If provided the projected surface mesh onto horizontal coordinates
212
!! and its associated rtree/pickers are cached under this name and
213
!! attached to the 'positions'. In this case when called again with
214
!! the same 'positions' and the same surface_name,
215
!! the same surface_element_list should again be provided.
216
character(len=*), optional, intent(in):: surface_name
218
type(scalar_field), dimension(from_field%dim):: from_field_components, to_field_components
221
assert(from_field%dim==to_field%dim)
223
do i=1, from_field%dim
224
from_field_components(i)=extract_scalar_field(from_field, i)
225
to_field_components(i)=extract_scalar_field(to_field, i)
228
call VerticalExtrapolationMultiple( from_field_components, to_field_components, &
229
positions, vertical_normal, surface_element_list, &
230
surface_name=surface_name)
232
end subroutine VerticalExtrapolationVector
234
subroutine VerticalExtrapolationMultiple(from_fields, to_fields, &
235
positions, vertical_normal, surface_element_list, surface_name)
236
!!< This sub extrapolates the values on a horizontal 2D surface
237
!!< in the vertical direction to 3D fields
238
!! The from_fields should be 3D fields of which only the values on the
239
!! 2D horizontal surface are used.
240
!! This version takes multiple from_fields at the same time and extrapolates
241
!! to to_fields, such that the surface search only has to be done once. This
242
!! will only work if all the from_fields are on the same mesh, and on all the
243
!! to_field are on the same (possibly a different) mesh.
244
type(scalar_field), dimension(:), intent(in):: from_fields
245
!! Resulting extrapolated field. May be the same field or a field on
246
!! a different mesh (different degree).
247
type(scalar_field), dimension(:), intent(inout):: to_fields
248
!! positions, and upward normal vector on the whole domain
249
type(vector_field), target, intent(inout):: positions
250
type(vector_field), target, intent(in):: vertical_normal
252
!! the surface elements (faces numbers) that make up the surface
253
integer, dimension(:), intent(in):: surface_element_list
254
!! If provided the projected surface mesh onto horizontal coordinates
255
!! and its associated rtree/pickers are cached under this name and
256
!! attached to the 'positions'. In this case when called again with
257
!! the same 'positions' and the same surface_name,
258
!! the same surface_element_list should again be provided.
259
character(len=*), optional, intent(in):: surface_name
261
character(len=FIELD_NAME_LEN):: lsurface_name
262
real, dimension(:,:), allocatable:: loc_coords
263
integer, dimension(:), allocatable:: seles
264
integer i, j, to_nodes, face
266
assert(size(from_fields)==size(to_fields))
267
assert(element_count(from_fields(1))==element_count(to_fields(1)))
268
do i=2, size(to_fields)
269
assert(to_fields(1)%mesh==to_fields(i)%mesh)
270
assert(from_fields(1)%mesh==from_fields(i)%mesh)
273
to_nodes=nowned_nodes(to_fields(1))
274
! local coordinates is one more than horizontal coordinate dim
275
allocate( seles(to_nodes), loc_coords(1:positions%dim, 1:to_nodes) )
277
if (present(surface_name)) then
278
lsurface_name=surface_name
280
lsurface_name="TempSurfaceName"
283
! project the positions of to_fields(1)%mesh into the horizontal plane
284
! and returns 'seles' (indices in surface_element_list) and loc_coords
285
! to tell where these projected nodes are found in the surface mesh
286
call horizontal_picker(to_fields(1)%mesh, positions, vertical_normal, &
287
surface_element_list, lsurface_name, &
290
! interpolate using the returned faces and loc_coords
291
do i=1, size(to_fields)
293
face=surface_element_list(seles(j))
294
call set(to_fields(i), j, &
295
dot_product( eval_shape( face_shape(from_fields(i), face), loc_coords(:,j) ), &
296
face_val( from_fields(i), face ) ))
300
if (IsParallel()) then
301
do i=1, size(to_fields)
302
call halo_update(to_fields(i))
306
if (.not. present(surface_name)) then
307
call remove_boundary_condition(positions, "TempSurfaceName")
310
end subroutine VerticalExtrapolationMultiple
312
subroutine horizontal_picker(mesh, positions, vertical_normal, &
313
surface_element_list, surface_name, &
315
!! Searches the nodes of 'mesh' in the surface mesh above.
316
!! Returns the surface elements 'seles' that each node lies under
317
!! and the loc_coords in this element of this node projected
318
!! upward (radially on the sphere) onto the surface mesh
321
type(mesh_type), intent(in):: mesh
322
!! a valid positions field for the whole domain, not necessarily on 'mesh'
323
!! for instance in a periodic domain, mesh is periodic and positions should not be
324
type(vector_field), target, intent(inout):: positions
325
!! upward normal vector on the whole domain
326
type(vector_field), target, intent(in):: vertical_normal
327
!! the surface elements (faces numbers) that make up the surface
328
integer, dimension(:), intent(in):: surface_element_list
329
!! The projected surface mesh onto horizontal coordinates
330
!! and its associated rtree/pickers are cached under this name and
331
!! attached to the 'positions'. When called again with
332
!! the same 'positions' and the same surface_name,
333
!! the same surface_element_list should again be provided.
334
character(len=*), intent(in):: surface_name
336
!! returned surface elements (face numbers in 'mesh')
337
!! and loc coords each node has been found in
338
!! size(seles)==size(loc_coords,2)==nowned_nodes(mesh)
339
integer, dimension(:), intent(out):: seles
340
real, dimension(:,:), intent(out):: loc_coords
342
type(vector_field):: mesh_positions
343
type(vector_field), pointer:: horizontal_positions
344
integer, dimension(:), pointer:: horizontal_mesh_list
345
real, dimension(:,:), allocatable:: horizontal_coordinate
346
real, dimension(vertical_normal%dim):: normal_vector
347
real, dimension(positions%dim):: xyz
348
integer:: i, stat, nodes
350
assert(.not. mesh_periodic(positions))
352
! search only the first owned nodes
353
nodes=nowned_nodes(mesh)
355
assert( size(seles)==nodes )
356
assert(size(loc_coords,1)==positions%dim)
357
assert( size(loc_coords,2)==nodes )
359
if (mesh==positions%mesh) then
360
mesh_positions=positions
361
! make mesh_positions indep. ref. of the field, so we can deallocate it
362
! safely without destroying positions
363
call incref(mesh_positions)
365
call allocate(mesh_positions, positions%dim, mesh, &
366
name='ToPositions_VerticalExtrapolation')
367
call remap_field(positions, mesh_positions, stat)
368
if (stat/=0 .and. stat/=REMAP_ERR_HIGHER_LOWER_CONTINUOUS .and. &
369
stat/=REMAP_ERR_UNPERIODIC_PERIODIC) then
370
! Mapping from higher order to lower order is allowed for coordinates
371
! (well depends on how the higher order is derived from the lower order)
373
! Mapping to periodic coordinates is ok in this case as we only need
374
! locations for the nodes individually (i.e. we don't care about elements
375
! in 'mesh') - the created horizontal_positions will be non-periodic
376
! So that we should be able to find nodes on the periodic boundary on
377
! either side. Using this to interpolate is consistent as long as
378
! the interpolated from field is indeed periodic.
379
FLAbort("Unknown error in remmaping positions in horizontal_picker.")
384
! create an array of the horizontally projected coordinates of the 'mesh' nodes
385
allocate( horizontal_coordinate(1:positions%dim-1, 1:nodes) )
386
if (have_option('/geometry/spherical_earth')) then
387
assert(mesh_positions%dim==3)
389
xyz=node_val(mesh_positions, i)
391
horizontal_coordinate(:,i) = &
392
map2horizontal_sphere( node_val(mesh_positions, i), +1.0)
394
horizontal_coordinate(:,i) = &
395
map2horizontal_sphere( node_val(mesh_positions, i), -1.0)
396
horizontal_coordinate(1,i)=horizontal_coordinate(1,i)+HEMISPHERE_SHIFT
400
assert( vertical_normal%field_type==FIELD_TYPE_CONSTANT )
401
normal_vector=node_val(vertical_normal,1)
404
horizontal_coordinate(:,i) = &
405
map2horizontal( node_val(mesh_positions, i), normal_vector )
409
call get_horizontal_positions(positions, surface_element_list, &
410
vertical_normal, surface_name, &
411
horizontal_positions, horizontal_mesh_list)
413
call picker_inquire( horizontal_positions, horizontal_coordinate, &
414
seles, loc_coords, global=.false. )
416
! in the spherical case some of the surface elements may be duplicated
417
! within horizontal positions, the returned seles should refer to entries
418
! in surface_element_list however - also check for nodes not found
421
seles(i)=horizontal_mesh_list(seles(i))
423
ewrite(-1,*) "For node with coordinate", node_val(mesh_positions, i)
424
ewrite(-1,*) "no top surface node was found."
425
FLAbort("Something wrong with the geometry.")
429
call deallocate(mesh_positions)
430
deallocate(horizontal_mesh_list)
432
end subroutine horizontal_picker
434
subroutine get_horizontal_positions(positions, surface_element_list, vertical_normal, surface_name, &
435
horizontal_positions, horizontal_mesh_list)
436
! returns a horizontal positions field over the surface mesh indicated by
437
! 'surface_element_list'. This field will be created and cached on 'positions'
438
! as a surface field attached to a dummy boundary condition under the name surface_name
439
type(vector_field), intent(inout):: positions
440
type(vector_field), intent(in):: vertical_normal
441
integer, dimension(:), intent(in):: surface_element_list
442
character(len=*), intent(in):: surface_name
444
! Returns the horizontal positions on a horizontal mesh, this mesh
445
! may have some of the faces in surface_element_list duplicated -
446
! this is only in the spherical case where the horizontal mesh
447
! consists of two disjoint regions representing the two hemispheres
448
! and surface elements near the equator may be represented in both.
449
! Therefore we also return a map between elements in the horizontal positions mesh
450
! and entries in surface_element_list. If ele is an element in the horizontal
451
! position mesh, then surface_element_list(horizontal_mesh(ele))
452
! is the face number in the full mesh
453
! horizontal_positions is a borrowed reference, don't allocate
454
! horizontal_mesh_list does need to be deallocated
455
type(vector_field), pointer:: horizontal_positions
456
integer, dimension(:), pointer:: horizontal_mesh_list
458
integer, dimension(:), pointer:: horizontal_sele_list
459
integer:: i, j, k, sele
461
if (.not. has_boundary_condition_name(positions, surface_name)) then
462
if (have_option('/geometry/spherical_earth')) then
463
call create_horizontal_positions_sphere(positions, &
464
surface_element_list, surface_name)
466
call create_horizontal_positions_flat(positions, &
467
surface_element_list, vertical_normal, surface_name)
471
horizontal_positions => extract_surface_field(positions, surface_name, &
472
trim(surface_name)//"HorizontalCoordinate")
474
! call vtk_write_fields('horizontal_mesh', 0, horizontal_positions, &
475
! horizontal_positions%mesh)
478
allocate( horizontal_mesh_list(1:element_count(horizontal_positions)) )
479
if (have_option('/geometry/spherical_earth')) then
480
call get_boundary_condition(positions, surface_name, &
481
surface_element_list=horizontal_sele_list)
482
! by construction the map can be obtained by running
483
! through surface_element_list twice (for each hemisphere)
484
i=1 ! position in horizontal_sele_list
485
sele=horizontal_sele_list(i) ! face we're looking for
488
do k=1, size(surface_element_list)
489
if (surface_element_list(k)==sele) then
490
! found matching position in surface_element_list
491
horizontal_mesh_list(i)=k
494
if (i>size(horizontal_sele_list)) exit outer_loop
495
sele=horizontal_sele_list(i)
500
if (i<=size(horizontal_sele_list)) then
501
! not all were found in 2 loops through surface_element_list
503
FLAbort("Internal error in horizontal mesh administration")
507
! no duplication: horizontal_mesh_list is simply the identity map
508
do i=1, size(horizontal_mesh_list)
509
horizontal_mesh_list(i)=i
513
end subroutine get_horizontal_positions
515
subroutine create_horizontal_positions_flat(positions, surface_element_list, vertical_normal, surface_name)
516
! adds a "boundary condition" to 'positions' with an associated vector surface field containing a dim-1
517
! horizontal coordinate field that can be used to map from the surface mesh specified
518
! by 'surface_element_list'. This "boundary condition" will be stored under the name 'surface_name'
519
! The horizontal coordinates are created by projecting out the component
520
! in the direction of 'vertical_normal' and then throwing out the x, y or z
521
! coordinate that is most aligned with 'vertical_normal'
522
type(vector_field), intent(inout):: positions
523
type(vector_field), intent(in):: vertical_normal
524
integer, dimension(:), intent(in):: surface_element_list
525
character(len=*), intent(in):: surface_name
527
type(mesh_type), pointer:: surface_mesh
528
type(vector_field):: horizontal_positions
529
integer, dimension(:), pointer:: surface_node_list
530
real, dimension(vertical_normal%dim):: normal_vector
533
assert(vertical_normal%field_type==FIELD_TYPE_CONSTANT)
534
normal_vector=node_val(vertical_normal, 1)
536
call add_boundary_condition_surface_elements(positions, &
537
name=surface_name, type="verticalextrapolation", &
538
surface_element_list=surface_element_list)
539
! now get back the created surface mesh
540
! and surface_node_list a mapping between node nos in the projected mesh and node nos in the original 'positions' mesh
541
call get_boundary_condition(positions, name=surface_name, &
542
surface_mesh=surface_mesh, surface_node_list=surface_node_list)
543
call allocate( horizontal_positions, dim=positions%dim-1, mesh=surface_mesh, &
544
name=trim(surface_name)//"HorizontalCoordinate" )
546
call insert_surface_field(positions, name=surface_name, &
547
surface_field=horizontal_positions)
549
do i=1, size(surface_node_list)
550
node=surface_node_list(i)
551
call set(horizontal_positions, i, &
552
map2horizontal(node_val(positions, node), normal_vector))
555
call deallocate( horizontal_positions )
557
end subroutine create_horizontal_positions_flat
559
subroutine create_horizontal_positions_sphere(positions, surface_element_list, surface_name)
560
! adds a "boundary condition" to 'positions' with an associated vector surface field containing a dim-1
561
! horizontal coordinate field that can be used to map from the surface mesh specified
562
! by 'surface_element_list'. This "boundary condition" will be stored under the name 'surface_name'
563
! The horizontal coordinates are created by a stereographic projection from the unit sphere to the
564
! xy-plane. The northern hemisphere is projected to the unit-circle (including some extra surface
565
! elements on the southern hemisphere around the equator). The southern hemisphere (including some
566
! northern surface elements near the equator) is also projected to a unit circle but translated
567
! away from the origin to not overlap with the projected northern hemisphere. Thus the created
568
! projected horizontal positions field consists of two disjoint areas in the plane corresponding
569
! to both hemispheres, and elements in the original surface mesh (near the equator) may appear
570
! twice in the projected field.
571
type(vector_field), intent(inout):: positions
572
integer, dimension(:), intent(in):: surface_element_list
573
character(len=*), intent(in):: surface_name
575
type(vector_field), pointer:: horizontal_positions_north, horizontal_positions_south
576
type(vector_field):: horizontal_positions
577
type(mesh_type), pointer:: surface_mesh_north, surface_mesh_south
578
type(mesh_type):: surface_mesh
579
integer, dimension(:), pointer:: surface_element_list_north, surface_element_list_south
580
integer, dimension(:), allocatable:: surface_element_list_combined
581
integer:: nodes_north, elements_south, elements_north
584
! first create 2 separate horizontal coordinate fields
585
! for each hemisphere
587
call create_horizontal_positions_hemisphere(positions, &
588
surface_element_list, trim(surface_name)//'North', +1.0)
590
call create_horizontal_positions_hemisphere(positions, &
591
surface_element_list, trim(surface_name)//'South', -1.0)
593
! these are stored as bcs under the positions
594
! retrieve this information back:
596
call get_boundary_condition(positions, name=trim(surface_name)//'North', &
597
surface_mesh=surface_mesh_north, &
598
surface_element_list=surface_element_list_north)
599
call get_boundary_condition(positions, name=trim(surface_name)//'South', &
600
surface_mesh=surface_mesh_south, &
601
surface_element_list=surface_element_list_south)
603
horizontal_positions_north => extract_surface_field(positions, &
604
trim(surface_name)//'North', trim(surface_name)//"NorthHorizontalCoordinate")
605
horizontal_positions_south => extract_surface_field(positions, &
606
trim(surface_name)//'South', trim(surface_name)//"SouthHorizontalCoordinate")
608
! merge these 2 meshes
609
surface_mesh=merge_meshes( (/ surface_mesh_north, surface_mesh_south /) )
611
! and merge the positions field
612
call allocate( horizontal_positions, positions%dim-1, surface_mesh, &
613
trim(surface_name)//"HorizontalCoordinate" )
615
nodes_north=node_count(surface_mesh_north)
617
call set( horizontal_positions, i, &
618
node_val(horizontal_positions_north, i))
620
! but translate the projected southern hemisphere to the left
621
! to not overlap it with the northern hemisphere
622
do i=1, node_count(surface_mesh_south)
623
call set( horizontal_positions, nodes_north+i, &
624
node_val(horizontal_positions_south, i)+(/ HEMISPHERE_SHIFT, 0.0 /) )
627
! merge the surface element lists
628
! (note that this is different than the incoming surface_element_list
629
! as it will have some duplicate equatorial elements)
630
elements_north=size(surface_element_list_north)
631
elements_south=size(surface_element_list_south)
632
allocate(surface_element_list_combined(1:elements_north+elements_south))
633
surface_element_list_combined(1:elements_north)=surface_element_list_north
634
surface_element_list_combined(elements_north+1:)=surface_element_list_south
636
! finally the bc for the combined surface mesh:
637
! (note that this creates a new surface mesh different than
638
! the merged mesh, that we won't be using)
639
call add_boundary_condition_surface_elements( &
640
positions, name=surface_name, type="verticalextrapolation", &
641
surface_element_list=surface_element_list_combined)
643
! insert the horizontal positions under this bc
644
call insert_surface_field( positions, name=surface_name, &
645
surface_field=horizontal_positions)
647
! everything is safely stored, so we can deallocate our references
648
call deallocate(horizontal_positions)
649
call deallocate(surface_mesh)
650
deallocate(surface_element_list_combined)
652
! also we won't need the 2 hemisphere bcs anymore
653
call remove_boundary_condition( positions, trim(surface_name)//'North')
654
call remove_boundary_condition( positions, trim(surface_name)//'South')
656
end subroutine create_horizontal_positions_sphere
658
subroutine create_horizontal_positions_hemisphere(positions, &
659
surface_element_list, surface_name, hemi_sign)
660
type(vector_field), intent(inout):: positions
661
integer, dimension(:), intent(in):: surface_element_list
662
character(len=*), intent(in):: surface_name
663
real, intent(in):: hemi_sign
665
type(vector_field):: horizontal_positions
666
type(mesh_type), pointer:: surface_mesh
667
real, dimension(2):: xy
668
integer, dimension(:), pointer:: nodes, surface_node_list
669
integer:: i, j, sele, node
671
type(integer_set):: surface_element_set
673
call allocate(surface_element_set)
675
do i=1, size(surface_element_list)
676
sele=surface_element_list(i)
677
if (sele_is_on_hemisphere(sele)) then
678
call insert(surface_element_set, sele)
682
call add_boundary_condition_surface_elements(positions, &
683
name=surface_name, type="verticalextrapolation", &
684
surface_element_list=set2vector(surface_element_set))
685
call deallocate(surface_element_set)
687
call get_boundary_condition(positions, name=surface_name, &
688
surface_mesh=surface_mesh, surface_node_list=surface_node_list)
690
call allocate( horizontal_positions, dim=2, mesh=surface_mesh, &
691
name=trim(surface_name)//"HorizontalCoordinate" )
693
call insert_surface_field(positions, name=surface_name, &
694
surface_field=horizontal_positions)
696
do i=1, element_count(horizontal_positions)
697
nodes => ele_nodes(horizontal_positions, i)
699
node=surface_node_list( nodes(j) )
700
xy=map2horizontal_sphere(node_val(positions, node), hemi_sign)
701
assert(abs(xy(1))<HEMISPHERE_SHIFT/2.0)
702
call set( horizontal_positions, nodes(j), xy)
706
call deallocate(horizontal_positions)
710
logical function sele_is_on_hemisphere(sele)
711
integer, intent(in):: sele
713
real, dimension(face_loc(positions, sele)):: znodes
715
znodes=face_val(positions, 3, sele)
716
! if any of the nodes is on the hemisphere or *very* close to the equator
717
! (relative to the other nodes)
718
sele_is_on_hemisphere = any(hemi_sign * znodes > -sum(abs(znodes))*1e-6)
720
end function sele_is_on_hemisphere
722
end subroutine create_horizontal_positions_hemisphere
724
function map2horizontal(xyz, normal_vector)
725
real, dimension(:), intent(in):: xyz, normal_vector
726
real, dimension(size(xyz)-1):: map2horizontal
728
real, dimension(size(xyz)):: hxyz
729
integer:: i, c, takeout
731
! first subtract of the vertical component
732
hxyz=xyz-dot_product(xyz, normal_vector)*normal_vector
734
! then leave out the "most vertical" coordinate
735
takeout=maxloc(abs(normal_vector), dim=1)
739
if (i==takeout) cycle
740
map2horizontal(c)=hxyz(i)
744
end function map2horizontal
746
function map2horizontal_sphere(xyz, hemi_sign)
747
real, dimension(3), intent(in):: xyz
748
real, intent(in):: hemi_sign
749
real, dimension(2):: map2horizontal_sphere
754
map2horizontal_sphere=xyz(1:2)/(r+hemi_sign*xyz(3))
756
end function map2horizontal_sphere
758
function VerticalProlongationOperator(mesh, positions, vertical_normal, &
759
surface_element_list, surface_mesh)
760
!! creates a prolongation operator that prolongates values on
761
!! a surface mesh to a full mesh below using the same interpolation
762
!! as the vertical extrapolation code above. The transpose of this prolongation
763
!! operator can be used as a restriction/clustering operator
764
!! from the full mesh to the surface.
765
type(csr_matrix) :: VerticalProlongationOperator
766
!! the mesh to which to prolongate, its nodes are only considered as
767
!! a set of loose points
768
type(mesh_type), target, intent(in):: mesh
769
!! positions on the whole domain (doesn't have to be the same mesh)
770
type(vector_field), intent(inout):: positions
771
!! upward normal vector on the whole domain
772
type(vector_field), target, intent(in):: vertical_normal
773
!! the face numbers of the surface mesh
774
integer, dimension(:), intent(in):: surface_element_list
775
!! Optionally a surface mesh may be provided that represents the nodes
776
!! /from/ which to interpolate (may be of different signature than 'mesh')
777
!! Each column in the prolongator will correspond to a node in this surface_mesh.
778
!! If not provided, the "from mesh" is considered to consist of faces
779
!! given by surface_element_list (and same cont. and order as 'mesh').
780
!! In this case however empty columns, i.e. surface nodes from which no
781
!! value in the full mesh is interpolated, are removed and there is no
782
!! necessary relation between column numbering and surface node numbering.
783
type(mesh_type), intent(in), target, optional:: surface_mesh
785
type(csr_sparsity):: sparsity
786
type(mesh_type), pointer:: lsurface_mesh
787
real, dimension(:,:), allocatable:: loc_coords
788
real, dimension(:), allocatable:: mat
790
integer, dimension(:), pointer:: snodes
791
integer, dimension(:), allocatable:: seles, colm, findrm, snod2used_snod
792
integer i, j, k, rows, entries, count, snod, sele
793
! coefficient have to be at least this otherwise they're negligable in an interpolation
794
real, parameter :: COEF_EPS=1d-10
796
! only assemble the rows associted with nodes we own
797
rows=nowned_nodes(mesh)
799
! local coordinates is one more than horizontal coordinate dim
800
allocate( seles(rows), loc_coords(1:positions%dim, 1:rows) )
802
! project the positions of to_fields(1)%mesh into the horizontal plane
803
! and returns 'seles' (indices in surface_element_list) and loc_coords
804
! to tell where these projected nodes are found in the surface mesh
805
call horizontal_picker(mesh, positions, vertical_normal, &
806
surface_element_list, "TempSurfaceName", &
809
! count upper estimate for n/o entries for sparsity
812
entries=entries+face_loc(mesh, seles(i))
814
! preliminary matrix:
815
allocate( mat(1:entries), findrm(1:rows+1), &
818
if (.not. present(surface_mesh)) then
819
! We use the entire surface mesh of 'mesh'
820
lsurface_mesh => mesh%faces%surface_mesh
821
! Not all surface nodes may be used (i.e. interpolated from) - even
822
! within surface elements that /are/ in surface_element-list.
823
! We need a map between global surface node numbering and
824
! a consecutive numbering of used surface nodes.
825
! (this will be the column numbering)
826
allocate(snod2used_snod(1:node_count(lsurface_mesh)))
828
count=0 ! counts number of used surface nodes
830
lsurface_mesh => surface_mesh
833
entries=0 ! this time only count nonzero entries
836
! beginning of each row in mat
839
if (present(surface_mesh)) then
840
! element number within surface_mesh
843
! face number in 'mesh', i.e. element number within entire surface_mesh
844
sele=surface_element_list(seles(i))
847
snodes => ele_nodes(lsurface_mesh, sele)
850
coef=eval_shape(ele_shape(lsurface_mesh, sele), j, loc_coords(:,i))
852
if (abs(coef)>COEF_EPS) then
853
if (.not. present(surface_mesh)) then
854
if (snod2used_snod(snod)==0) then
855
! as of yet unused surface node
857
snod2used_snod(snod)=count
859
! this is the column index we're gonna use instead
860
snod=snod2used_snod(snod)
871
if (present(surface_mesh)) then
872
! we haven't counted used surface nodes, instead we're using all
873
! nodes of surface mesh as columns
874
count=node_count(surface_mesh)
877
call allocate(sparsity, rows, count, &
878
entries, diag=.false., name="VerticalProlongationSparsity")
879
sparsity%findrm=findrm
880
sparsity%colm=colm(1:entries)
881
! for lots of applications it's good to have sorted rows
882
call sparsity_sort(sparsity)
884
call allocate(VerticalProlongationOperator, sparsity, &
885
name="VerticalProlongationOperator")
886
call deallocate(sparsity)
888
! as the sparsity has been sorted the ordering of mat(:) no longer
889
! matches that of sparsity%colm, however it still matches the original
892
do k=findrm(i), findrm(i+1)-1
894
call set(VerticalProlongationOperator, i, j, mat(k))
898
if (.not. present(surface_mesh)) then
899
deallocate( snod2used_snod )
901
deallocate( findrm, colm, mat )
902
deallocate( seles, loc_coords )
904
call remove_boundary_condition(positions, "TempSurfaceName")
906
end function VerticalProlongationOperator
908
subroutine vertical_element_ordering(ordered_elements, face_normal_gravity, optimal_ordering)
909
!!< Calculates an element ordering such that each element is
910
!!< is preceded by all elements above it.
911
integer, dimension(:), intent(out):: ordered_elements
912
!! need to supply face_normal_gravity matrix,
913
!! created by compute_face_normal_gravity() subroutine below
914
type(csr_matrix), intent(in):: face_normal_gravity
915
!! returns .true. if an optimal ordering is found, i.e there are no
916
!! cycles, i.o.w. elements that are (indirectly) above and below each other
917
!! at the same time (deck of cards problem).
918
logical, optional, intent(out):: optimal_ordering
920
type(dynamic_bin_type) dbin
921
real, dimension(:), pointer:: inn
922
integer, dimension(:), pointer:: neigh
923
integer, dimension(:), allocatable:: bin_list
924
integer i, j, elm, bin_no
927
assert( size(ordered_elements)==size(face_normal_gravity,1) )
929
! create binlist, i.e. assign each element to a bin, according to
930
! the number of elements above it
931
allocate(bin_list(1:size(ordered_elements)))
932
do i=1, size(ordered_elements)
933
neigh => row_m_ptr(face_normal_gravity, i)
934
inn => row_val_ptr(face_normal_gravity, i)
935
! elements with no element above it go in bin 1
936
! elements with n elements above it go in bin n+1
937
! neigh>0 so we don't count exterior boundary faces
938
bin_list(i)=count( inn<-VERTICAL_INTEGRATION_EPS .and. neigh>0 )+1
941
call allocate(dbin, bin_list)
944
do i=1, size(ordered_elements)
945
! pull an element from the first non-empty bin
946
! (hopefully an element with no unprocessed elements above)
947
call pull_element(dbin, elm, bin_no)
948
ordered_elements(i)=elm
949
! if this is bin one then it is indeed an element with no unprocessed
950
! elements above, otherwise issue a warning
951
if (bin_no>1) warning=.true.
953
! update elements below:
956
neigh => row_m_ptr(face_normal_gravity, elm)
957
inn => row_val_ptr(face_normal_gravity, elm)
959
if (inn(j)>VERTICAL_INTEGRATION_EPS .and. neigh(j)>0) then
960
! element neigh(j) is below element i, therefore now has one
961
! less unprocessed element above it, so can be moved to
963
if (.not. element_pulled(dbin, neigh(j))) then
964
! but only if neigh(j) itself hasn't been selected yet
965
! (which might happen for imperfect vertical orderings)
966
call move_element(dbin, neigh(j), bin_list(neigh(j))-1)
973
! this warning may be reduced (in verbosity level) if it occurs frequently:
974
ewrite(-1,*) "Warning: vertical_element_ordering has detected a cycle."
975
ewrite(-1,*) "(deck of cards problem). This may reduce the efficiency"
976
ewrite(-1,*) "of your vertically sweeping solve."
979
if (present(optimal_ordering)) then
980
optimal_ordering=.not. warning
983
call deallocate(dbin)
985
end subroutine vertical_element_ordering
987
subroutine compute_face_normal_gravity(face_normal_gravity, &
988
positions, vertical_normal)
989
!!< Returns a matrix where A_ij is the inner product of the face normal
990
!!< and the gravity normal vector of the face between element i and j.
991
type(csr_matrix), intent(out):: face_normal_gravity
992
type(vector_field), target, intent(in):: positions, vertical_normal
994
type(mesh_type), pointer:: mesh
995
real, dimension(:), pointer:: face_normal_gravity_val
996
real, dimension(:), allocatable:: detwei_f
997
real, dimension(:,:), allocatable:: face_normal, gravity_normal
998
integer, dimension(:), pointer:: neigh, faces
1000
integer sngi, nloc, i, k
1002
mesh => positions%mesh
1003
call allocate(face_normal_gravity, mesh%faces%face_list%sparsity)
1004
call zero(face_normal_gravity)
1006
sngi=face_ngi(mesh, 1)
1007
nloc=ele_loc(mesh,1)
1008
allocate( detwei_f(1:sngi), &
1009
face_normal(1:positions%dim, 1:sngi), &
1010
gravity_normal(1:positions%dim, 1:sngi))
1012
do i=1, element_count(mesh)
1013
! elements adjacent to element i
1014
! this is a row (column indices) in the mesh%faces%face_list matrix
1015
neigh => ele_neigh(mesh, i)
1016
! the surrounding faces
1017
! this is a row (integer values) in the mesh%faces%face_list matrix
1018
faces => ele_faces(mesh, i)
1020
if (neigh(k)>i .or. neigh(k)<=0) then
1021
! only handling neigh(k)>i to ensure anti-symmetry of the matrix
1022
! (and more efficient of course)
1023
call transform_facet_to_physical(positions, faces(k), &
1024
detwei_f=detwei_f, &
1026
gravity_normal=face_val_at_quad(vertical_normal, faces(k))
1028
! inner product of face normal and vertical normal
1029
! integrated over face
1030
inn=sum(matmul(face_normal*gravity_normal, detwei_f))/area
1031
if (neigh(k)>0) then
1032
call set(face_normal_gravity, i, neigh(k), inn)
1033
call set(face_normal_gravity, neigh(k), i, -inn)
1035
! exterior surface: matrix entry does not have valid
1036
! column index, still want to store its value, so we
1038
face_normal_gravity_val => row_val_ptr(face_normal_gravity, i)
1039
face_normal_gravity_val(k)=inn
1046
end subroutine compute_face_normal_gravity
1048
subroutine vertical_integration_scalar(from_field, to_field, &
1049
positions, vertical_normal, surface_element_list, rhs)
1050
!!< See description vertical_integration_multiple
1051
type(scalar_field), intent(in):: from_field
1052
type(scalar_field), intent(in):: to_field
1053
type(vector_field), intent(in):: positions, vertical_normal
1054
integer, dimension(:), intent(in):: surface_element_list
1055
type(scalar_field), optional, intent(in):: rhs
1057
type(scalar_field) to_fields(1)
1059
to_fields=(/ to_field /)
1060
if (present(rhs)) then
1061
call vertical_integration_multiple( (/ from_field /), to_fields, &
1062
positions, vertical_normal, surface_element_list, rhs=(/ rhs /) )
1064
call vertical_integration_multiple( (/ from_field /), to_fields, &
1065
positions, vertical_normal, surface_element_list)
1068
end subroutine vertical_integration_scalar
1070
subroutine vertical_integration_vector(from_field, to_field, &
1071
positions, vertical_normal, surface_element_list, rhs)
1072
!!< See description vertical_integration_multiple
1073
type(vector_field), intent(in):: from_field
1074
type(vector_field), intent(in):: to_field
1075
type(vector_field), intent(in):: positions, vertical_normal
1076
integer, dimension(:), intent(in):: surface_element_list
1077
type(vector_field), optional, intent(in):: rhs
1079
type(scalar_field), dimension(from_field%dim):: from_field_components, &
1080
to_field_components, rhs_components
1083
assert(from_field%dim==to_field%dim)
1085
do i=1, from_field%dim
1086
from_field_components(i)=extract_scalar_field(from_field, i)
1087
to_field_components(i)=extract_scalar_field(to_field, i)
1088
if (present(rhs)) then
1089
rhs_components(i)=extract_scalar_field(rhs, i)
1093
if (present(rhs)) then
1094
call vertical_integration_multiple( from_field_components, &
1095
to_field_components, positions, vertical_normal, &
1096
surface_element_list, rhs=rhs_components)
1098
call vertical_integration_multiple( from_field_components, &
1099
to_field_components, positions, vertical_normal, &
1100
surface_element_list)
1103
end subroutine vertical_integration_vector
1105
subroutine vertical_integration_multiple(from_fields, to_fields, &
1106
positions, vertical_normal, surface_element_list, rhs)
1107
!!< This subroutine solves: dP/dz=rhs using DG
1108
!!< It can be used for vertical integration downwards (dP/dz=0) as a drop
1109
!!< in replacement of VerticalExtrapolation hence its similar interface.
1110
!!< The field P is the to_field. A boundary condition is given by the
1111
!!< from_field. Again completely similar to VerticalExtrapolation, it may
1112
!!< be defined as a surface field on surface elements given by
1113
!!< surface_element_list or it may be a field on the complete mesh
1114
!!< in which case only its values on these surface elements are used.
1115
!!< vertical_normal specifies the direction in which to integrate (usually downwards)
1116
!!< If not specified rhs is assumed zero.
1118
!!< This version accepts multiple from_fields, to_fields and rhs
1119
type(scalar_field), dimension(:), intent(in):: from_fields
1120
type(scalar_field), dimension(:), intent(inout):: to_fields
1121
type(vector_field), intent(in):: positions, vertical_normal
1122
integer, dimension(:), intent(in):: surface_element_list
1123
type(scalar_field), dimension(:), optional, intent(in):: rhs
1125
type(csr_matrix) face_normal_gravity
1126
type(element_type), pointer:: ele_shp, face_shp, x_face_shp
1127
real, dimension(:), pointer:: inn
1128
real, dimension(:,:,:), allocatable:: surface_rhs, dele_shp
1129
real, dimension(:,:), allocatable:: ele_mat, face_mat, ele_rhs
1130
real, dimension(:), allocatable:: detwei, detwei_f
1131
integer, dimension(:), pointer:: neigh, ele_nds, faces, face_lnds
1132
integer, dimension(:), allocatable:: ordered_elements, face_nds, face_nds2
1133
integer nloc, snloc, ngi, sngi
1134
integer i, j, k, f, f2, elm, it, noit
1135
logical optimal_ordering, from_surface_fields
1137
assert( size(from_fields)==size(to_fields) )
1139
! computes inner product of face normal and gravity (see above)
1140
call compute_face_normal_gravity(face_normal_gravity, &
1141
positions, vertical_normal)
1143
! determine an ordering for the elements based on this
1144
allocate( ordered_elements(1:element_count(positions)) )
1145
call vertical_element_ordering(ordered_elements, face_normal_gravity, &
1148
! General initalisation
1149
!-----------------------
1150
! various grid numbers
1151
nloc=ele_loc(to_fields(1), 1)
1152
snloc=face_loc(to_fields(1), 1)
1153
ngi=ele_ngi(positions, 1)
1154
sngi=face_ngi(positions, 1)
1156
ele_shp => ele_shape(to_fields(1), 1)
1157
face_shp => face_shape(to_fields(1), 1)
1158
x_face_shp => face_shape(positions, 1)
1159
! various allocations:
1161
surface_rhs(1:snloc, 1:size(to_fields), 1:surface_element_count(positions)), &
1162
ele_mat(1:nloc, 1:nloc), ele_rhs(1:nloc,1:size(to_fields)), &
1163
dele_shp(1:nloc, 1:ngi, 1:positions%dim), &
1164
face_mat(1:snloc, 1:snloc), face_nds(1:snloc), face_nds2(1:snloc), &
1165
detwei(1:ngi), detwei_f(1:sngi))
1167
if (element_count(from_fields(1))==size(surface_element_list)) then
1168
! from_fields are fields over the surface mesh only
1169
! so we're using all of its values:
1170
from_surface_fields=.true.
1171
! check the other fields as well:
1172
do k=2, size(from_fields)
1173
assert( element_count(from_fields(k))==size(surface_element_list) )
1176
! from_fields are on the full mesh and we only extract its values
1177
! at the specified surface_elements
1179
from_surface_fields=.false.
1183
! Compute contribution of exterior surface integral (boundary condition) to rhs
1184
!-----------------------
1185
do i=1, size(surface_element_list)
1186
f=surface_element_list(i)
1187
call transform_facet_to_physical(positions, f, detwei_f)
1188
face_mat=-shape_shape(face_shp, face_shp, detwei_f)
1189
do k=1, size(from_fields)
1190
if (from_surface_fields) then
1191
! we need to use ele_val, where i is the element number in the surface_mesh
1192
surface_rhs(:,k,f)=surface_rhs(:,k,f)+ &
1193
matmul(face_mat, ele_val(from_fields(k), i))
1195
! we can simply use face_val with face number f
1196
surface_rhs(:,k,f)=surface_rhs(:,k,f)+ &
1197
matmul(face_mat, face_val(from_fields(k), f))
1203
!-----------------------
1204
if (optimal_ordering) then
1211
do i=1, element_count(positions)
1213
elm=ordered_elements(i)
1215
! construct diagonal matrix block for this element
1216
call transform_to_physical(positions, elm, &
1217
shape=ele_shp, dshape=dele_shp, detwei=detwei)
1218
ele_mat=shape_vector_dot_dshape(ele_shp, &
1219
ele_val_at_quad(vertical_normal,elm), &
1223
if (present(rhs)) then
1224
do k=1, size(to_fields)
1225
ele_rhs(:,k)=shape_rhs(ele_shp, detwei*ele_val_at_quad(rhs(k), elm))
1232
! then add contribution of surface integrals of incoming
1233
! faces to the rhs and matrix
1234
neigh => row_m_ptr(face_normal_gravity, elm)
1235
inn => row_val_ptr(face_normal_gravity, elm)
1236
faces => ele_faces(positions, elm)
1238
if (inn(j)<-VERTICAL_INTEGRATION_EPS) then
1239
call transform_facet_to_physical(positions, faces(j), &
1241
face_mat=-shape_shape(face_shp, face_shp, detwei_f)*inn(j)
1243
face_nds=face_global_nodes(to_fields(1), faces(j))
1244
face_lnds => face_local_nodes(to_fields(1)%mesh, faces(j))
1245
ele_mat(face_lnds,face_lnds)=ele_mat(face_lnds,face_lnds)+face_mat
1247
if (neigh(j)>0) then
1248
! face of element neigh(j), facing elm:
1249
f2=ele_face(positions, neigh(j), elm)
1250
face_nds2=face_global_nodes(to_fields(1), f2)
1252
do k=1, size(to_fields)
1253
ele_rhs(face_lnds,k)=ele_rhs(face_lnds,k)+ &
1254
matmul(face_mat, node_val(to_fields(k), face_nds2))
1257
! note that we've already multiplied with face_mat above, but not with inn(j)
1258
ele_rhs(face_lnds,:)=ele_rhs(face_lnds,:)+surface_rhs(:,:,faces(j))*inn(j)
1263
call invert(ele_mat)
1265
! compute values for the to_fields:
1266
ele_nds => ele_nodes(to_fields(1), elm)
1267
do k=1, size(to_fields)
1268
call set( to_fields(k), ele_nds, matmul(ele_mat, ele_rhs(:,k)) )
1273
call deallocate(face_normal_gravity)
1275
end subroutine vertical_integration_multiple
1277
subroutine vertical_extrapolation_module_check_options
1279
if (have_option("/geometry/ocean_boundaries")) then
1280
if (.not. have_option("/physical_parameters/gravity")) then
1281
ewrite(-1,*) "If you select /geometry/ocean_boundaries, you also need to "//&
1282
&"set /physical_parameters/gravity"
1283
FLExit("Missing gravity!")
1287
end subroutine vertical_extrapolation_module_check_options
1289
end module vertical_extrapolation_module
added
removed
assemble/Vertical_Extrapolation.F90
