<p><b>mhoffman@lanl.gov</b> 2013-01-30 20:40:57 -0700 (Wed, 30 Jan 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Some first modifications necessary for modeling ice shelves.<br>
<br>
1.  Added a marine basal mass balance field.  This is input as a time series (but can only have 1 time level) and is defined across the entire domain.  Each time step the mask is used to 0 this field for grounded ice cells.  Therefore this field can be defined for areas that will unground in the future.  For now this is simply an input field, but eventually it can come from a coupler or ocean model.  Note: grounded ice basal mass balance is not yet defined because a column diffusion temperature solve does not yet exist.<br>
<br>
2. Set beta to 0 for floating ice.  LifeV may already do this, but doing it here for completeness.<br>
<br>
3. Adjusting Registry to output lower and upper surface which are useful for visualizing floating ice.<br>
<br>
(also added messages to stderr when trying to run a lifev velocity solver without a lifev build.  such a message would have saved me some time recently.)<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry        2013-01-30 22:34:57 UTC (rev 2399)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry        2013-01-31 03:40:57 UTC (rev 2400)
@@ -87,8 +87,8 @@
 dim nBetaTimeSlices nBetaTimeSlices
 dim nSfcAirTempTimeSlices nSfcAirTempTimeSlices
 dim nBasalHeatFluxTimeSlices nBasalHeatFluxTimeSlices
+dim nMarineBasalMassBalTimeSlices nMarineBasalMassBalTimeSlices
 
-
 % VARIABLES =====================
 %
 % var persistence type  name_in_file  ( dims )  time_levs iro-  name_in_code struct super-array array_class
@@ -191,14 +191,15 @@
 var persistent real    basalHeatFlux ( nCells ) 0 - basalHeatFlux mesh - -
 var persistent real    sfcAirTemp ( nCells ) 0 - sfcAirTemp mesh - -
 var persistent real    sfcMassBal ( nCells ) 0 - sfcMassBal mesh - -
+var persistent real    marineBasalMassBal ( nCells ) 0 - marineBasalMassBal mesh - -
 
 % Time-varying forcing for B.C.
 var persistent real    betaTimeSeries ( nBetaTimeSlices nCells ) 0 ir betaTimeSeries mesh - -
 var persistent real    basalHeatFluxTimeSeries ( nBasalHeatFluxTimeSlices nCells ) 0 ir basalHeatFluxTimeSeries mesh - -
 var persistent real    sfcAirTempTimeSeries ( nSfcAirTempTimeSlices nCells ) 0 ir sfcAirTempTimeSeries mesh - -
 var persistent real    sfcMassBalTimeSeries ( nSfcMassBalTimeSlices nCells ) 0 ir sfcMassBalTimeSeries mesh - -
+var persistent real    marineBasalMassBalTimeSeries ( nMarineBasalMassBalTimeSlices nCells ) 0 ir marineBasalMassBalTimeSeries mesh - -
 
-
 % STATE VARIABLES =====================
 
 % Prognostic variables: read from input, saved in restart, and written to output
@@ -222,8 +223,8 @@
 var persistent real    uReconstructMeridional ( nVertLevels nCells Time ) 2 - uReconstructMeridional state - -
 var persistent real    layerThickness ( nVertLevels nCells Time ) 2 - layerThickness state - -
 var persistent real    layerThicknessEdge ( nVertLevels nEdges Time ) 2 - layerThicknessEdge state - -
-var persistent real    lowerSurface ( nCells Time ) 2 - lowerSurface state - -
-var persistent real    upperSurface ( nCells Time ) 2 - upperSurface state - -
+var persistent real    lowerSurface ( nCells Time ) 2 o lowerSurface state - -
+var persistent real    upperSurface ( nCells Time ) 2 o upperSurface state - -
 var persistent real    temperatureUpperBC ( nCells Time ) 2 - temperatureUpperBC state - -
 var persistent real    temperatureLowerBC ( nCells Time ) 2 - temperatureLowerBC state - -
 % Other diagnostic variables: neither read nor written to any files

Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_annual_forcing.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_annual_forcing.F        2013-01-30 22:34:57 UTC (rev 2399)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_annual_forcing.F        2013-01-31 03:40:57 UTC (rev 2400)
@@ -94,8 +94,8 @@
       ! local variables
       !
       !-----------------------------------------------------------------
-      real (kind=RKIND), dimension(:,:), pointer :: sfcMassBalTimeSeries, sfcAirTempTimeSeries, betaTimeSeries, basalHeatFluxTimeSeries
-      real (kind=RKIND), dimension(:), pointer :: sfcMassBal, sfcAirTemp, beta, basalHeatFlux
+      real (kind=RKIND), dimension(:,:), pointer :: sfcMassBalTimeSeries, sfcAirTempTimeSeries, betaTimeSeries, basalHeatFluxTimeSeries, marineBasalMassBalTimeSeries
+      real (kind=RKIND), dimension(:), pointer :: sfcMassBal, sfcAirTemp, beta, basalHeatFlux, marineBasalMassBal
       integer :: iYear, ierr, nCells
       
 
@@ -105,11 +105,13 @@
       sfcAirTempTimeSeries =&gt; grid % sfcAirTempTimeSeries % array
       betaTimeSeries =&gt; grid % betaTimeSeries % array
       basalHeatFluxTimeSeries =&gt; grid % basalHeatFluxTimeSeries % array
+      marineBasalMassBalTimeSeries =&gt; grid % marineBasalMassBalTimeSeries % array
 
       sfcMassBal =&gt; grid % sfcMassBal % array
       sfcAirTemp =&gt; grid % sfcAirTemp % array
       beta =&gt; grid % beta % array
       basalHeatFlux =&gt; grid % basalHeatFlux % array
+      marineBasalMassBal =&gt; grid % marineBasalMassBal % array
 
       nCells = grid % nCells
       
@@ -127,6 +129,7 @@
       call assign_annual_forcing_to_1D_array(sfcAirTemp, sfcAirTempTimeSeries, iYear, nCells)
       call assign_annual_forcing_to_1D_array(beta, betaTimeSeries, iYear, nCells)
       call assign_annual_forcing_to_1D_array(basalHeatFlux, basalHeatFluxTimeSeries, iYear, nCells)
+      call assign_annual_forcing_to_1D_array(marineBasalMassBal, marineBasalMassBalTimeSeries, iYear, nCells)
 
    end subroutine land_ice_assign_annual_forcing!}}}
 
@@ -197,7 +200,6 @@
    end subroutine assign_annual_forcing_to_1D_array
 
 
-
 !***********************************************************************
 
 end module land_ice_annual_forcing

Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F        2013-01-30 22:34:57 UTC (rev 2399)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F        2013-01-31 03:40:57 UTC (rev 2400)
@@ -104,7 +104,8 @@
       !call lifeV and set the grid of the velocity solver
       call velocity_solver_init_mpi(domain % dminfo % comm)
 #else
-      write(*,*) &quot;LifeV needed to run L1L2, FO and Stokes models&quot;
+      write(*,*) &quot;LifeV needed to run L1L2, FO and Stokes models&quot;  ! stdout
+      write(0,*) &quot;LifeV needed to run L1L2, FO and Stokes models&quot;  ! stderr
       err = 1
       return
 #endif
@@ -163,6 +164,7 @@
       call phg_init(domain % dminfo % comm)
 #else
       write(*,*) &quot;PHG needed to run Stokes model&quot;
+      write(0,*) &quot;PHG needed to run Stokes model&quot;
       err = 1
       return
 #endif
@@ -341,6 +343,7 @@
       call mpas_timer_stop(&quot;velocity_solver_set_grid_data&quot;)
 #else
       write(*,*) &quot;LifeV needed to run L1L2, FO and Stokes models&quot;
+      write(0,*) &quot;LifeV needed to run L1L2, FO and Stokes models&quot;
       err = 1
       return
 #endif

Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F        2013-01-30 22:34:57 UTC (rev 2399)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F        2013-01-31 03:40:57 UTC (rev 2400)
@@ -112,7 +112,10 @@
       ! local variables
       !
       !-----------------------------------------------------------------
+      integer nVertLevels
 
+     nVertLevels = mesh % nVertLevels
+
      ! 0 tendency
      layerThickness_tend = 0.0
 
@@ -156,11 +159,23 @@
         ! Do nothing - don't add the MB
      case default
         ! Add surface mass balance to tendency
-        ! NOTE: Need to decide where to put this for layer by layer thickness tendency, and how to deal with negative SMB that eliminates top layer or all ice (check for negative thickness?)
-        layerThickness_tend(1,:) = layerThickness_tend(1,:) + mesh % sfcMassBal % array  ! (tendency in meters per year)
-        ! THIS MIGHT RESULT IN  NEGATIVE LAYER THICKNESS!
+        ! TODO: Need to decide where to put this for layer by layer thickness tendency, and how to deal with negative SMB that eliminates top layer or all ice (check for negative thickness?)
+        layerThickness_tend(1,:) = layerThickness_tend(1,:) + mesh % sfcMassBal % array   ! (tendency in meters per year)
+        ! TODO THIS MIGHT RESULT IN  NEGATIVE LAYER THICKNESS!
 
-        ! Do basal mass balance too?
+        ! Add basal mass balance to tendency
+        !    Apply marineBasalMassBal to floating ice only.  Apply groundedBasalMassBal to grounded ice only.
+        !    TODO: more complicated treatment at GL?
+        !    Assign 0.0 values to the marineBasalMassBal where the ice is grounded.
+        !    It's ok to overwrite the values with 0's here, because each time step
+        !    we get a fresh copy of the array from the annual_forcing subroutine.
+        where ( MASK_IS_GROUNDED(state % cellMask % array) )
+           mesh % marineBasalMassBal % array = 0.0
+        end where
+        layerThickness_tend(nVertLevels,:) = layerThickness_tend(nVertLevels,:) &amp;
+                                    + mesh % marineBasalMassBal % array    ! (tendency in meters per year)
+        ! TODO Add in grounded ice basal mass balance once temperature diffusion is calculated
+        ! TODO THIS MIGHT RESULT IN  NEGATIVE LAYER THICKNESS!
 
      end select
 
@@ -346,7 +361,7 @@
       !-----------------------------------------------------------------
       integer :: iCell, iLevel, nCells, nVertLevels
       real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &amp;
-        lowerSurface, bedTopography, layerThicknessFractions
+        lowerSurface, bedTopography, layerThicknessFractions, beta
       integer, dimension(:), pointer :: cellMask
       real (kind=RKIND), dimension(:,:), pointer :: layerThickness
 
@@ -359,15 +374,23 @@
       lowerSurface =&gt; state % lowerSurface % array
       bedTopography =&gt; state % bedTopography % array
       layerThickness =&gt; state % layerThickness % array
+      beta =&gt; mesh % beta % array
       cellMask =&gt; state % cellMask % array
 
 
       ! Calculate masks - needs to happen before calculating lower surface so we know where the ice is floating
       call land_ice_calculate_mask(mesh, state, err)
 
-      ! no ice shelves -- lower surface same as bed topography
-      ! Eventually this should consider floatation!
+      ! Update beta before the velocity solve occurs, now that we have the new state and its mask.
+      !    It's ok to overwrite the beta values with 0's here, because each time step
+      !    we get a fresh copy of the array from the annual_forcing subroutine.
+      !    Note: some velocity solvers may do this on their own, but we are doing it here for completeness.
       where ( MASK_IS_FLOATING(cellMask) )
+         beta = 0.0
+      end where
+
+      ! Lower surface is based on floatation for floating ice.  For grounded ice it is the bed.
+      where ( MASK_IS_FLOATING(cellMask) )
          lowerSurface = config_sealevel - thickness * (config_ice_density / config_ocean_density)
       elsewhere
          lowerSurface = bedTopography

</font>
</pre>