<p><b>mhoffman@lanl.gov</b> 2013-03-27 09:31:01 -0600 (Wed, 27 Mar 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Setting up a simple dome ice sheet on a spherical mesh. This directory has instructions and tools for doing this.<br>
</p><hr noshade><pre><font color="gray">Added: branches/land_ice_projects/test_cases/dome_sphere/README
===================================================================
--- branches/land_ice_projects/test_cases/dome_sphere/README (rev 0)
+++ branches/land_ice_projects/test_cases/dome_sphere/README 2013-03-27 15:31:01 UTC (rev 2672)
@@ -0,0 +1,52 @@
+Instructions for making a spherical mesh for the land ice core. Right now, this is pretty simple - it uses the serial grid generator, it does not cull out cells, and it does not provide variable resolution.
+
+------------------
+A. Run global_scvt
+------------------
+** Compile the code at: trunk/grid_gen/global_scvt
+Instructions for compiling and running can be found in the Appendix of the nhyd_atm users guide at: /trunk/documents/users_guide/mpas_nhyd_atm_users_guide.tex
+Note: for using gfortran, I had to remove the compiler flags -fsecond-underscore from the Makefile to get it to link properly.
+
+** Run global_scvt 3 consecutive times using the default locs.dat file:
+First increase n_scvt_iterations in the namelist file to something large (e.g. 10000).
+Run it:
+> ./grid_gen
+
+Use the locs.dat.out that it generates to seed a new run:
+> cp locs.dat.out locs.dat
+Find the number of generating points available
+> wc -l locs.dat
+(you want the number printed minus one since there is one header line in the file)
+Edit np in namelist.input to be set to the new number of generating points (642 in this case).
+Start the new run:
+> ./grid_gen
+
+Repeat the process one more time. This is the minimum resolution you might want to try running the ice sheet model at. You can run it a few more times if you want higher resolution. Use MpasDraw to visualize the spherical mesh.
+
+
+------------
+B. Run basin
+------------
+The ocean core tool 'basin' will expand the sphere from a unit sphere (that global_scvt produces) to an earth-sized sphere. The tool can also cull out part of the mesh, but we won't bother with that yet.
+It can be found here: branches/ocean_projects/basin
+Compile basin following instructions in that directory
+Copy the grid.nc file from step A into the directory that has basin in it.
+Run basin: > ./map
+You don't need to modify the namelist for basin. It is doing a bunch of stuff that we don't need, but the only thing that matters is that these config settings are set:
+ on_a_sphere = 'YES'
+ expand_from_unit_sphere = .true.
+ sphere_radius = 6.37122e6
+
+-------------
+C. Setup the land ice grid and add IC
+-------------
+Copy ocean.nc from step B into your work directory for the land ice run
+
+Use the python script for setting up a land ice grid (works on both planar and spherical meshes).
+land_ice_projects/grid_tools/add_land_ice_variables_to_mpas_grid.py ocean.nc
+It produces land_ice_grid.nc.
+
+Add the I.C. to the grid:
+Run: python setup_spherical_dome_initial_conditions.py
+
+You can use MpasDraw to visualize the grid, specifically the thickness variable. You can also modify the python script to modify the I.C.
Added: branches/land_ice_projects/test_cases/dome_sphere/namelist.basin
===================================================================
--- branches/land_ice_projects/test_cases/dome_sphere/namelist.basin (rev 0)
+++ branches/land_ice_projects/test_cases/dome_sphere/namelist.basin 2013-03-27 15:31:01 UTC (rev 2672)
@@ -0,0 +1,33 @@
+&basin
+ nVertLevelsMOD = 10
+ on_a_sphere = 'YES'
+ expand_from_unit_sphere = .true.
+ sphere_radius = 6.37122e6
+ zLevel_thickness = 'POP_40_zLevel'
+ bottom_topography = 'flat_bottom'
+ initial_conditions = 'uniform_TS'
+ eliminate_inland_seas=.false.
+ load_woce_IC = .false.
+ write_OpenDX_flag = .false.
+ monthly_forcing = .false.
+ check_mesh = .true.
+ cut_domain_from_sphere = .false.
+ solid_boundary_in_y = .false.
+ solid_boundary_in_x = .false.
+ top_layers_without_land = 3
+
+ ! These variables may be used for acc wind amplification
+ amplify_acc_wind = .false.
+ amp_wind_factor = 2.0
+ amp_wind_center_lat = -35.0
+ amp_wind_spread_lat = 3.0
+
+ ! These variables are not needed for realistic global topography:
+
+ ! h_total_max = 2000.0
+ ! u_src_max = 0.1
+ ! f0 = -1.1e-4
+ ! beta = 1.4e-11
+ ! omega = 7.29212e-5
+ ! Lx = 3200.0e3
+/
Added: branches/land_ice_projects/test_cases/dome_sphere/namelist.input.global_scvt
===================================================================
--- branches/land_ice_projects/test_cases/dome_sphere/namelist.input.global_scvt (rev 0)
+++ branches/land_ice_projects/test_cases/dome_sphere/namelist.input.global_scvt 2013-03-27 15:31:01 UTC (rev 2672)
@@ -0,0 +1,9 @@
+&domains
+ np = 2562
+ locs_as_xyz = .true.
+ n_scvt_iterations = 100
+ eps = 0.000000001
+ l2_conv = .true.
+ inf_conv = .false.
+ min_dx = 120000.0
+/
Added: branches/land_ice_projects/test_cases/dome_sphere/setup_spherical_dome_initial_conditions.py
===================================================================
--- branches/land_ice_projects/test_cases/dome_sphere/setup_spherical_dome_initial_conditions.py (rev 0)
+++ branches/land_ice_projects/test_cases/dome_sphere/setup_spherical_dome_initial_conditions.py 2013-03-27 15:31:01 UTC (rev 2672)
@@ -0,0 +1,107 @@
+#!/usr/bin/python
+# Generate initial conditions for dome land ice test case on a spherical mesh
+
+import sys, numpy
+try:
+ from Scientific.IO.NetCDF import NetCDFFile
+ netCDF_module = 'Scientific.IO.NetCDF'
+except ImportError:
+ try:
+ from netCDF4 import Dataset as NetCDFFile
+ netCDF_module = 'netCDF4'
+ except ImportError:
+ print 'Unable to import any of the following python modules:'
+ print ' Scientific.IO.NetCDF </font>
<font color="blue"> netcdf4 '
+ print 'One of them must be installed.'
+ raise ImportError('No netCDF module found')
+from math import sqrt
+
+# Check to see if a grid file was specified on the command line.
+# If not, land_ice_grid.nc is used.
+if len(sys.argv) > 1:
+ if sys.argv[1][0] == '-': # The filename can't begin with a hyphen
+ print '</font>
<font color="black">Usage: python setup_dome_initial_conditions.py [GRID.NC]</font>
<font color="blue">If no filename is supplied, land_ice_grid.nc will be used.'
+ sys.exit(0)
+ else:
+ gridfilename = sys.argv[1]
+else:
+ gridfilename = 'land_ice_grid.nc'
+
+
+
+# Open the file, get needed dimensions
+try:
+ gridfile = NetCDFFile(gridfilename,'r+')
+ if (netCDF_module == 'Scientific.IO.NetCDF'):
+ nVertLevels = gridfile.dimensions['nVertLevels']
+ else:
+ nVertLevels = len(gridfile.dimensions['nVertLevels'])
+ if nVertLevels != 9:
+ print 'nVerLevels in the supplied file was ', nVertLevels, '. Were you expecting 9?'
+ # Get variables
+ latCell = gridfile.variables['latCell'][:]
+ lonCell = gridfile.variables['lonCell'][:]
+ thickness = gridfile.variables['thickness'][:]
+ bedTopography = gridfile.variables['bedTopography'][:]
+ normalVelocity = gridfile.variables['normalVelocity'][:]
+ layerThicknessFractions = gridfile.variables['layerThicknessFractions'][:]
+ temperature = gridfile.variables['temperature'][:]
+ # Get b.c. variables
+ beta = gridfile.variables['betaTimeSeries'][:]
+ SMB = gridfile.variables['sfcMassBalTimeSeries'][:]
+ Tsfc = gridfile.variables['sfcAirTempTimeSeries'][:]
+ G = gridfile.variables['basalHeatFluxTimeSeries'][:]
+ BMB = gridfile.variables['marineBasalMassBalTimeSeries'][:]
+except:
+ sys.exit('Error: The grid file specified is either missing or lacking needed dimensions/variables.')
+
+
+# Assign variable values for dome
+
+# center of dome at South Pole
+
+
+
+# Set default value for non-dome cells
+thickness[:] = 0.0
+
+h0 = 5000.0 # thickness in me of center of dome
+r0 = -90.0 * numpy.pi/180.0
+r1 = -60.0 * numpy.pi/180.0
+ind = (latCell < r1)
+
+# Calculate the dome thickness for cells within the desired radius (thickness will be NaN otherwise)
+#thickness[0, ind] = h0 * ( 1.0 - (latCell[ind] / r0 )**2 )**0.5
+thickness[0, ind] = h0 * ( (latCell[ind] - r1) / (r0 - r1) )**0.5
+# zero velocity everywhere
+normalVelocity[:] = 0.0
+# flat bed at sea level
+bedTopography[:] = 0.0
+# constant, arbitrary temperature, degrees C
+temperature[:] = -20.0
+# Setup layerThicknessFractions
+layerThicknessFractions[:] = 1.0 / nVertLevels
+
+# boundary conditions
+# Sample values to use, or comment these out for them to be 0.
+beta[:] = 10000.
+#SMB[:] = 2.0/1000.0 * (thickness[:] + bedTopography[:]) - 1.0 # units: m/yr, lapse rate of 1 m/yr with 0 at 500 m
+Tsfc[:,0] = 0.0
+G = 0.05
+#BMB[:] = -20.0 # units: m/yr
+
+# Reassign the modified numpy array values back into netCDF variable objects
+gridfile.variables['thickness'][:] = thickness
+gridfile.variables['normalVelocity'][:] = normalVelocity
+gridfile.variables['bedTopography'][:] = bedTopography
+gridfile.variables['temperature'][:] = temperature
+gridfile.variables['layerThicknessFractions'][:] = layerThicknessFractions
+gridfile.variables['betaTimeSeries'][:] = beta
+gridfile.variables['sfcMassBalTimeSeries'][:] = SMB
+gridfile.variables['sfcAirTempTimeSeries'][:] = Tsfc
+gridfile.variables['basalHeatFluxTimeSeries'][:] = G
+gridfile.variables['marineBasalMassBalTimeSeries'][:] = BMB
+
+gridfile.close()
+print 'Successfully added dome initial conditions to ', gridfilename
+
</font>
</pre>