Merge pull request #14 from NCAR/lipscomb/inversion2

Last round of changes before the CISM release

This PR includes many changes in various parts of the code, in preparation for the CISM.1 public release. Among the major changes:

* Added a scheme to invert for basal traction parameters and basal melt rates.  This was the original purpose of the branch. But along the way, I made some general changes that were desirable for the release, and it would have been hard to separate those changes from the inversion scheme.  So I ended up keeping all the changes in one branch. The inversion scheme itself is not yet ready for release, but I organized the inversion calculations so that it will be straightforward to redact them.

* Set new default values of several options and parameters, including:
  - which_ho_assemble_taud/beta/bfric = 1 (local assembly)
  - which_ho_flotation_function = 2 (linear)
  - beta_grounded_min = 1 Pa/(m/yr)
  - basal_mbal = 1 (include bmlt in continuity)
  - upn = 3 (smallest value for which the model will run)
  - cliff_timescale = 10 yr (gradual rather than abrupt calving of marine cliffs)
  - pmp_offset = 5 deg (colder initial bed for spin-ups)

* Changed some config files of test cases to be consistent with new defaults; also removed '[DOME-TEST]' and similar headings that are no longer needed

* Added a calving_rate field, plus tavg fields for acab_applied, bmlt_applied, calving_rate

* Cleaned up various diagnostics; added an option to write mass budget terms in Gt/yr

* Reworked eigencalving and damage-based calving; eigencalving rates are now based on stresses instead of strain rates

* Simplified the calving interfaces and converted the calving module to SI units

* Changed the local assembly of taud and beta terms, to make the model more robust with less flickering at land-terminating margins

* Modified the hybrid gradient calculation to compute gradients at edges with ice-covered land above ice-free ocean; accordingly; the lateral shelf term at these edges was removed

* Modified the f_flotation calculation for ice-free land

Following this PR, various chunks of code related to inversion, the plume model, and the grounding-line parameterization will be redacted for the release.

AUTHORS

@@ -1,20 +1,21 @@
-The following authors (listed alphabetically) have contributed to this version of Glimmer-CISM.
+The following authors (listed alphabetically) have contributed to this version of CISM.
 Affiliations shown with an asterisk (*) are no longer current.
 
 Erin Barker		Los Alamos National Laboratory (*)
 Tim Bocek               University of Montana, Missoula (*)
-Josh Campbell           University of Montana, Missoula
+Josh Campbell           University of Montana, Missoula (*)
 Katherine J. Evans      Oak Ridge National Laboratory
-Jeremy Fyke             Los Alamos National Laboratory
+Jeremy Fyke             Los Alamos National Laboratory (*)
 Glen Granzow            University of Montana, Missoula
 Magnus Hagdorn	        School of GeoSciences, University of Edinburgh
 Brian Hand              University of Montana, Missoula (*)
 Felix Hebeler           University of Zurich(*) 
 Matthew Hoffman         Los Alamos National Laboratory
 Jesse Johnson           University of Montana, Missoula
 Irina Kalashnikova      Sandia National Laboratories
+Gunter Leguy            Los Alamos National Laboratory (*), National Center for Atmospheric Research
 Jean-Francois Lemieux   New York University (*)
-William Lipscomb        Los Alamos National Laboratory
+William Lipscomb        Los Alamos National Laboratory (*), National Center for Atmospheric Research
 Daniel Martin		Lawrence Berkeley National Laboratory
 Jeffrey A. Nichols      Oak Ridge National Laboratory
 Ryan Nong               Sandia National Laboratories (*)
@@ -26,6 +27,6 @@ Ian Rutt		Dept. of Geography, Swansea University
 William Sacks           National Center for Atmospheric Research
 Andrew Salinger         Sandia National Laboratories
 James B. White III      Oak Ridge National Laboratory (*)
-Jon Wolfe               National Center for Atmospheric Research (*)
+Jon Wolfe               National Center for Atmospheric Research (*) 
 Patrick Worley          Oak Ridge National Laboratory
 Timothy Wylie           University of Montana, Missoula (*)

cism_driver/cism_front_end.F90

@@ -216,8 +216,7 @@ subroutine cism_init_dycore(model)
 
     call t_startf('initial_write_diagnostics')
     call glide_write_diagnostics(model,        time,       &
-                                 tstep_count = model%numerics%tstep_count,  &
-                                 minthick_in = model%numerics%thklim*thk0)  ! m
+                                 tstep_count = model%numerics%tstep_count)
     call t_stopf('initial_write_diagnostics')
 
   end if ! whichdycore .ne. DYCORE_BISICLES
@@ -346,8 +345,7 @@ subroutine cism_run_dycore(model)
 
       call t_startf('write_diagnostics')
       call glide_write_diagnostics(model,        time,       &
-                                   tstep_count = model%numerics%tstep_count,  &
-                                   minthick_in = model%numerics%thklim*thk0)  ! m
+                                   tstep_count = model%numerics%tstep_count)
       call t_stopf('write_diagnostics')
 
       ! update time from dycore advance

cism_driver/eismint_forcing.F90

@@ -99,9 +99,6 @@ subroutine eismint_initialise(eismint_climate,config)
        eismint_climate%nmsb(1) = eismint_climate%nmsb(1) / acab_scale
        eismint_climate%nmsb(2) = eismint_climate%nmsb(2) / acab_scale
 
-    case(4)   ! MISMIP-1
-       eismint_climate%nmsb(1) = eismint_climate%nmsb(1) / acab_scale
-
     end select
 
   end subroutine eismint_initialise
@@ -188,35 +185,12 @@ subroutine eismint_readconfig(eismint_climate, config)
        return
     end if
 
-    !mismip tests
-
-    !TODO - Assign reasonable default values if not present in config file
-
-    call GetSection(config,section,'MISMIP-1')
-    if (associated(section)) then
-       eismint_climate%eismint_type = 4
-       dummy=>NULL()
-       call GetValue(section,'temperature',dummy,2)
-       if (associated(dummy)) then
-           eismint_climate%airt = dummy
-           deallocate(dummy)
-           dummy=>NULL()
-       end if
-       call GetValue(section,'massbalance',dummy,3)
-       if (associated(dummy)) then
-           eismint_climate%nmsb = dummy
-           deallocate(dummy)
-           dummy=>NULL()
-       end if
-       return
-    end if
-
     !exact verification
-    !TODO - Is this test currently supported?
+    !TODO - Is this test still supported?
 
     call GetSection(config,section,'EXACT')
     if (associated(section)) then
-       eismint_climate%eismint_type = 5
+       eismint_climate%eismint_type = 4
        dummy=>NULL()
        call GetValue(section,'temperature',dummy,2)
        if (associated(dummy)) then
@@ -227,53 +201,8 @@ subroutine eismint_readconfig(eismint_climate, config)
        return
     end if
 
-    ! Standard higher-order tests
-    ! These do not require EISMINT-type input parameters.
-
-    call GetSection(config,section,'DOME-TEST')
-    if (associated(section)) then
-        return
-    end if
-
-    call GetSection(config,section,'ISMIP-HOM-TEST')
-    if (associated(section)) then
-        return
-    end if
-
-    call GetSection(config,section,'SHELF-TEST')
-    if (associated(section)) then 
-        return
-    end if
-
-    call GetSection(config,section,'STREAM-TEST')
-    if (associated(section)) then 
-        return
-    end if
-
-    call GetSection(config,section,'ROSS-TEST')
-    if (associated(section)) then 
-        return
-    end if
-
-    call GetSection(config,section,'GIS-TEST')
-    if (associated(section)) then 
-        return
-    end if
-
-    call GetSection(config,section,'MISMIP+')
-    if (associated(section)) then
-        return
-    end if
-
-    call GetSection(config,section,'MISMIP')
-    if (associated(section)) then
-        return
-    end if
-
-    !TODO - Any other allowed tests to add here?
-
-    ! Abort if one of the above cases has not been specified.
-    call write_log('No EISMINT forcing selected',GM_FATAL)
+    ! Other tests (DOME-TEST, ISMIP-HOM-TEST, etc.) do not require
+    ! EISMINT-type input parameters.
 
   end subroutine eismint_readconfig
 
@@ -422,10 +351,6 @@ subroutine eismint_massbalance(eismint_climate,model,time)
        end do
 
     case(4)
-       !mismip 1
-       model%climate%acab = eismint_climate%nmsb(1)
-
-    case(5)
        !verification 
        call exact_surfmass(eismint_climate,model,time,1.d0,eismint_climate%airt(2))
 
@@ -493,9 +418,6 @@ subroutine eismint_surftemp(eismint_climate,model,time)
        end do
 
     case(4)
-       model%climate%artm = eismint_climate%airt(1)
-
-    case(5)
        !call both massbalance and surftemp at the same time to save computing time. 
        call exact_surfmass(eismint_climate,model,time,0.d0,eismint_climate%airt(2))
     end select

libglad/glad_initialise.F90

@@ -263,8 +263,7 @@ subroutine glad_i_initialise_gcm(config,           instance,         &
 
     call glide_write_diagnostics(instance%model,                  &
                                  instance%model%numerics%time,    &
-                                 tstep_count = instance%model%numerics%tstep_count,  &
-                                 minthick_in = instance%model%numerics%thklim*thk0)  ! m 
+                                 tstep_count = instance%model%numerics%tstep_count)
 
     ! Write netCDF output for this instance
 

libglad/glad_timestep.F90

@@ -255,8 +255,7 @@ subroutine glad_i_tstep_gcm(time,            instance,       &
 
           call glide_write_diagnostics(instance%model,                  &
                                        instance%model%numerics%time,    &
-                                       tstep_count = instance%model%numerics%tstep_count,  &
-                                       minthick_in = instance%model%numerics%thklim*thk0)  ! m 
+                                       tstep_count = instance%model%numerics%tstep_count)
 
           ! write netCDF output
 

libglide/glam_strs2.F90

@@ -3386,7 +3386,7 @@ subroutine bodyset(ew,  ns,  up,           &
            elseif( whichbabc == HO_BABC_BETA_CONSTANT .or. whichbabc == HO_BABC_YIELD_PICARD     .or. &
                    whichbabc == HO_BABC_BETA_LARGE    .or. whichbabc == HO_BABC_BETA_EXTERNAL    .or. &
                    whichbabc == HO_BABC_POWERLAW      .or. whichbabc == HO_BABC_COULOMB_FRICTION .or. &
-                   whichbabc == HO_BABC_COULOMB_CONST_BASAL_FLWA .or. whichbabc == HO_BABC_SIMPLE) then
+                   whichbabc == HO_BABC_COULOMB_POWERLAW_SCHOOF .or. whichbabc == HO_BABC_SIMPLE) then
 
                 bcflag = (/1,1/)              ! flag for specififed stress at bed: Tau_zx = beta * u_bed,
                                               ! where beta is MacAyeal-type traction parameter
@@ -3601,7 +3601,7 @@ subroutine bodyset(ew,  ns,  up,           &
            elseif( whichbabc == HO_BABC_BETA_CONSTANT .or. whichbabc == HO_BABC_YIELD_PICARD     .or. &
                    whichbabc == HO_BABC_BETA_LARGE    .or. whichbabc == HO_BABC_BETA_EXTERNAL    .or. &
                    whichbabc == HO_BABC_POWERLAW      .or. whichbabc == HO_BABC_COULOMB_FRICTION .or. &
-                   whichbabc == HO_BABC_COULOMB_CONST_BASAL_FLWA .or. whichbabc == HO_BABC_SIMPLE) then
+                   whichbabc == HO_BABC_COULOMB_POWERLAW_SCHOOF .or. whichbabc == HO_BABC_SIMPLE) then
 
                 bcflag = (/1,1/)              ! flag for specififed stress at bed: Tau_zx = beta * u_bed,
                                               ! where beta is MacAyeal-type traction parameter