<p><b>mhoffman@lanl.gov</b> 2012-07-23 11:39:48 -0600 (Mon, 23 Jul 2012)</p><p>BRANCH COMMIT -- land ice<br>
Made modifications needed to perform the ice2sea elevation feedback experiments.  Code compiles but I have not tested it.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/Registry        2012-07-23 14:00:25 UTC (rev 2040)
+++ branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/Registry        2012-07-23 17:39:48 UTC (rev 2041)
@@ -78,6 +78,7 @@
 dim R3 3
 dim FIFTEEN 15
 dim TWENTYONE 21
+dim TEN 10
 dim vertexDegree vertexDegree
 dim nVertLevels nVertLevels
 % nvertLevels+2 is only needed if attaching b.c. to a 3-d field.  We are unlikely to use this, but I am retaining it for now. MJH 04/19/12
@@ -198,9 +199,21 @@
 var persistent real    sfcAirTempTimeSeries ( nSfcAirTempTimeSlices nCells ) 0 ir sfcAirTempTimeSeries mesh - -
 var persistent real    sfcMassBalTimeSeries ( nSfcMassBalTimeSlices nCells ) 0 ir sfcMassBalTimeSeries mesh - -
 
+% variables needed for ice2sea elevation feedback experiment
+% field indicating which cells fall in the 'north' or 'south' 
+var persistent integer isNorth ( nCells ) 0 ir isNorth mesh - - 
+% the reference SMB that comes from the backwards looking average of the adjusted SMB - used for determining b
+var persistent real    SMB_reference ( nCells ) 0 r SMB_reference mesh - - 
+% the last 10 years of the adjusted SMB - used to calculate SMB_reference
+var persistent real    SMB_adjusted10 ( TEN nCells ) 0 r SMB_adjusted10 mesh - - 
+% a copy of the initial thickness field - used to calculate the thickness difference at each time step
+var persistent real    upperSurfaceIC ( nCells ) 0 r upperSurfaceIC mesh - - 
 
 % STATE VARIABLES =====================
 
+% a copy of the final SMB field that is used that is written to output.  It would be better to move sfcMassBal from grid to state, but doing it this way does not require editing the input files.
+var persistent real    SMBout ( nCells Time ) 2 o SMBout state - - 
+
 % Prognostic variables: read from input, saved in restart, and written to output
 var persistent real    thickness ( nCells Time ) 2 iro thickness state - -
 var persistent real    bedTopography ( nCells Time ) 2 ir bedTopography state - -

Modified: branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/mpas_land_ice_annual_forcing.F
===================================================================
--- branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/mpas_land_ice_annual_forcing.F        2012-07-23 14:00:25 UTC (rev 2040)
+++ branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/mpas_land_ice_annual_forcing.F        2012-07-23 17:39:48 UTC (rev 2041)
@@ -34,22 +34,37 @@
    !
    !--------------------------------------------------------------------
 
-   public :: land_ice_assign_annual_forcing
+   public :: land_ice_assign_annual_forcing, land_ice_assign_annual_forcing_init
 
    !--------------------------------------------------------------------
    !
    ! Private module variables
    !
    !--------------------------------------------------------------------
+   integer :: iYearPrevious, smb_counter
 
 
-
 !***********************************************************************
 
 contains
 
 !***********************************************************************
 !
+!  routine land_ice_assign_annual_forcing_init
+!
+!
+!-----------------------------------------------------------------------
+
+   subroutine land_ice_assign_annual_forcing_init()
+
+      iYearPrevious = -999
+      smb_counter = 0
+
+   end subroutine land_ice_assign_annual_forcing_init
+
+
+!***********************************************************************
+!
 !  routine land_ice_assign_annual_forcing
 !
 !&gt; \brief   Determines the forcing array used for the annual forcing.
@@ -61,7 +76,7 @@
 !
 !-----------------------------------------------------------------------
 
-   subroutine land_ice_assign_annual_forcing(timeStamp, grid, err)!{{{
+   subroutine land_ice_assign_annual_forcing(timeStamp, grid, state, err)!{{{
 
       !-----------------------------------------------------------------
       !
@@ -70,6 +85,7 @@
       !-----------------------------------------------------------------
 
       type(MPAS_Time_type), intent(in) :: timeStamp
+      type(state_type), intent(in) :: state
 
       !-----------------------------------------------------------------
       !
@@ -94,9 +110,10 @@
       ! local variables
       !
       !-----------------------------------------------------------------
-      real (kind=RKIND), dimension(:,:), pointer :: sfcMassBalTimeSeries, sfcAirTempTimeSeries, betaTimeSeries, basalHeatFluxTimeSeries
-      real (kind=RKIND), dimension(:), pointer :: sfcMassBal, sfcAirTemp, beta, basalHeatFlux
-      integer :: iYear, ierr, nCells
+      real (kind=RKIND), dimension(:,:), pointer :: sfcMassBalTimeSeries, sfcAirTempTimeSeries, betaTimeSeries, basalHeatFluxTimeSeries, SMB_adjusted10
+      real (kind=RKIND), dimension(:), pointer :: sfcMassBal, sfcAirTemp, beta, basalHeatFlux, SMB_reference
+      integer :: iYear, ierr, nCells, i
+      real b
       
 
       err = 0
@@ -112,6 +129,9 @@
       basalHeatFlux =&gt; grid % basalHeatFlux % array
 
       nCells = grid % nCells
+
+      SMB_reference =&gt; grid % SMB_reference % array
+      SMB_adjusted10 =&gt; grid % SMB_adjusted10 % array
       
       if (config_forcing_frequency .eq. 'Static') then
          iYear = 1  ! Static forcing uses the first entry in the time series at all time steps of the simulation
@@ -122,12 +142,56 @@
          err = ierr
       endif
 
-      ! Call this separately for each field - there will be cases where one field is time varying and others are constant, and this routine will deal with that properly. 
-      call assign_annual_forcing_to_1D_array(sfcMassBal, sfcMassBalTimeSeries, iYear, nCells)
-      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)
+      print *, 'Setting up SMB forcing for year ', iYear
+      if (iYear .ne. iYearPrevious) then  !only assign fields only when the year has actually changed
+        print *, 'New SMB being calculated'
+        ! Call this separately for each field - there will be cases where one field is time varying and others are constant, and this routine will deal with that properly. 
+        call assign_annual_forcing_to_1D_array(sfcMassBal, sfcMassBalTimeSeries, iYear, nCells)
+        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)
 
+        !ice2sea elevation feedback modifications
+        do i = 1, nCells
+          ! get proper b value
+          if (grid % isNorth % array(i) .eq. 1) then
+             if (SMB_reference(i) .ge. 0.0) then
+                b = 0.085
+             else
+                b = 0.543
+             endif
+          else  ! south side
+             if (SMB_reference(i) .ge. 0.0) then
+                b = 0.063
+             else
+                b = 1.890
+             endif
+          endif
+          ! scale from dimensional units of kg / m^3 / yr to 1 / yr
+          b = b / config_ice_density  
+
+          ! adjust SMB
+          sfcMassBal(i) = sfcMassBal(i) + b * (state%upperSurface%array(i) - grid % upperSurfaceIC % array(i) )
+        enddo
+
+        ! Calculate new SMB reference value
+        smb_counter = smb_counter + 1 
+        if( smb_counter &lt;= 10 )then
+            SMB_adjusted10(smb_counter,:) = sfcMassBal ! if less then first 10 yrs, keep storing 
+            SMB_reference = sum( SMB_adjusted10, dim=1) / real(smb_counter)
+        else
+            SMB_adjusted10(1:9,:) = SMB_adjusted10(2:10,:) ! move the old smb fields down in the stack and drop the oldest
+            SMB_adjusted10(10,:) = sfcMassBal  ! put the one from this time step calc on top
+            SMB_reference = sum( SMB_adjusted10, dim=1) / 10.0d0 
+        endif
+
+
+      endif ! check for a new year
+
+      iYearPrevious = iYear  ! assign this for reference next time round
+
+
+
    end subroutine land_ice_assign_annual_forcing!}}}
 
    !--------------------------------------------------------------------

Modified: branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/mpas_land_ice_mpas_core.F
===================================================================
--- branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/mpas_land_ice_mpas_core.F        2012-07-23 14:00:25 UTC (rev 2040)
+++ branches/land_ice_projects/ice2sea_elev_feedback/src/core_land_ice/mpas_land_ice_mpas_core.F        2012-07-23 17:39:48 UTC (rev 2041)
@@ -180,7 +180,11 @@
       call land_ice_setup_vertical_coords(mesh)
       currTime = mpas_get_clock_time(clock, MPAS_NOW, err_tmp)
       err = ior(err, err_tmp)
-      call land_ice_assign_annual_forcing(currTime, mesh, err_tmp)
+      ! ice2sea - init annual forcing module to do ice2sea stuff
+      call land_ice_assign_annual_forcing_init()
+      ! ice2sea initialize smb storage field to 0
+      mesh % SMB_adjusted10 % array = 0.0
+      call land_ice_assign_annual_forcing(currTime, mesh, state, err_tmp)
       err = ior(err, err_tmp)
 
       ! Init for FCT tracer advection
@@ -213,9 +217,13 @@
       ! \todo: skip this step if a velocity field was supplied in the I.C. input file
       call mpas_timer_start(&quot;initial velocity calculation&quot;)
       call land_ice_vel_solve(mesh, state, err)
-      call mpas_timer_stop(&quot;initial velocity calculation&quot;)
+      call mpas_timer_stop(&quot;initial velocity calculation&quot;)      
 
 
+      ! ice2sea - make a copy of the initial upper Surface field
+      mesh % upperSurfaceIC % array = state % upperSurface % array
+
+
       if(err == 1) then
           print *, &quot;An error has occurred in mpas_init_block. Aborting...&quot;
           call mpas_dmpar_abort(dminfo)
@@ -270,7 +278,9 @@
          call mpas_timer_start(&quot;assign forcing fields&quot;)
          block_ptr =&gt; domain % blocklist
          do while(associated(block_ptr))
-           call land_ice_assign_annual_forcing(currTime, block_ptr % mesh, ierr)
+           call land_ice_assign_annual_forcing(currTime, block_ptr % mesh, block_ptr % state % time_levs(1) % state, ierr)
+           ! assign the modified SMB field to an output field
+           block_ptr % state % time_levs(2) % state % SMBout % array = block_ptr % mesh % sfcMassBal % array
            block_ptr =&gt; block_ptr % next
          end do
          call mpas_timer_stop(&quot;assign forcing fields&quot;)

</font>
</pre>