Merge pull request #10 from NCAR/hgoelzer/tempunstag

Add initialization with unstaggered temperature from glide

libglide/glide_setup.F90

@@ -689,11 +689,12 @@ subroutine print_options(model)
          'constant in time      ', &
          'prognostic enthalpy   ' /)
 
-    character(len=*), dimension(0:3), parameter :: temp_init = (/ &
+    character(len=*), dimension(0:4), parameter :: temp_init = (/ &
          'set to 0 C                  ', &
          'set to surface air temp     ', &
          'linear vertical profile     ', &
-         'advective-diffusive balance ' /)
+         'advective-diffusive balance ',&
+         'temp from external file     ' /)
 
     character(len=*), dimension(0:3), parameter :: flow_law = (/ &
          'const 1e-16 Pa^-n a^-1      ', &

libglide/glide_types.F90

@@ -95,6 +95,7 @@ module glide_types
   integer, parameter :: TEMP_INIT_ARTM = 1
   integer, parameter :: TEMP_INIT_LINEAR = 2
   integer, parameter :: TEMP_INIT_ADVECTIVE_DIFFUSIVE = 3
+  integer, parameter :: TEMP_INIT_EXTERNAL = 4
 
   integer, parameter :: FLWA_CONST_FLWA = 0
   integer, parameter :: FLWA_PATERSON_BUDD_CONST_TEMP = 1
@@ -368,6 +369,7 @@ module glide_types
     !> \item[1] Initialize temperature to surface air temperature
     !> \item[2] Initialize temperature with a linear profile in each column
     !> \item[3] Initialize temperature with an advective-diffusive balance in each column
+    !> \item[4] Initialize temperature from external file
     !> \end{description}
 
     !> Method for calculating flow factor $A$:
@@ -1183,6 +1185,7 @@ module glide_types
     !       whereas bheatflx is defined as positive downward.
 
     real(dp),dimension(:,:,:),pointer :: temp => null()      !> 3D temperature field.
+    real(dp),dimension(:,:,:),pointer :: tempunstag => null()!> 3D temperature field unstaggered in the vertical.
     real(dp),dimension(:,:),  pointer :: bheatflx => null()  !> basal heat flux (W/m^2) (geothermal, positive down)
     real(dp),dimension(:,:,:),pointer :: flwa => null()      !> Glen's flow factor $A$.
     real(dp),dimension(:,:,:),pointer :: dissip => null()    !> interior heat dissipation rate, divided by rhoi*Ci (deg/s)
@@ -1934,6 +1937,9 @@ subroutine glide_allocarr(model)
        call coordsystem_allocate(model%general%ice_grid, upn, model%temper%dissip)
     else    ! glam/glissade dycore
        allocate(model%temper%temp(0:upn,1:ewn,1:nsn))
+       ! tempunstag has the same horizontal grid as the glam/glissade temp, but a
+       ! vertical axis like the glide temp
+       allocate(model%temper%tempunstag(upn,1:ewn,1:nsn))
        call coordsystem_allocate(model%general%ice_grid, upn-1, model%temper%flwa)
        call coordsystem_allocate(model%general%ice_grid, upn-1, model%temper%dissip)
     endif
@@ -1943,6 +1949,8 @@ subroutine glide_allocarr(model)
     model%temper%temp(:,:,:) = unphys_val  ! large negative number
     model%temper%flwa(:,:,:) = unphys_val
     model%temper%dissip(:,:,:) = 0.d0
+    if (associated(model%temper%tempunstag)) &
+         model%temper%tempunstag(:,:,:) = unphys_val
 
     call coordsystem_allocate(model%general%ice_grid,  model%temper%bheatflx)
     call coordsystem_allocate(model%general%ice_grid,  model%temper%bwat)
@@ -2264,6 +2272,8 @@ subroutine glide_deallocarr(model)
 
     if (associated(model%temper%temp)) &
         deallocate(model%temper%temp)
+    if (associated(model%temper%tempunstag)) &
+        deallocate(model%temper%tempunstag)
     if (associated(model%temper%bheatflx)) &
         deallocate(model%temper%bheatflx)
     if (associated(model%temper%bwat)) &

libglide/glide_vars.def

@@ -1216,6 +1216,15 @@ standard_name: land_ice_temperature_stag
 load:          1
 coordinates:   lon lat
 
+[tempunstag]
+dimensions:    time, level, y1, x1
+units:         degree_Celsius
+long_name:     ice temperature on unstaggered vertical levels
+data:          data%temper%tempunstag(up,1:data%general%ewn,1:data%general%nsn)
+standard_name: land_ice_temperature_unstag
+load:          1
+coordinates:   lon lat
+
 [bpmp]
 dimensions:    time, y1, x1
 units:         degree_Celsius

libglissade/glissade.F90

@@ -364,7 +364,8 @@ subroutine glissade_initialise(model, evolve_ice)
                              model%climate%acab*thk0/tim0,                          & ! m/s
                              model%temper%bheatflx,                                 & ! W/m^2, positive down
                              model%temper%pmp_offset,                               & ! deg C
-                             model%temper%temp)                                       ! deg C
+                             model%temper%temp,                                     & ! deg C
+                             model%temper%tempunstag)                                 ! deg C
 
     if (model%options%gthf == GTHF_COMPUTE) then
        call not_parallel(__FILE__,__LINE__)
@@ -985,6 +986,8 @@ subroutine glissade_thermal_solve(model, dt)
        bmlt_ground_unscaled,   & ! basal melt rate for grounded ice (m/s)
        bwat_unscaled             ! basal water thickness (m)
 
+    integer :: up
+
     call t_startf('glissade_thermal_solve')
 
     if (main_task .and. verbose_glissade) print*, 'Call glissade_therm_driver'
@@ -1034,6 +1037,17 @@ subroutine glissade_thermal_solve(model, dt)
     model%basal_melt%bmlt_ground(:,:) = bmlt_ground_unscaled(:,:) * tim0/thk0
     model%temper%bwat(:,:) = bwat_unscaled(:,:) / thk0
 
+    ! Update tempunstag as sigma weighted interpolation from temp to layer interfaces
+    do up = 2, model%general%upn-1
+      model%temper%tempunstag(up,:,:) = model%temper%temp(up-1,:,:) +       &
+           (model%temper%temp(up,:,:) - model%temper%temp(up-1,:,:)) *      &
+           (model%numerics%sigma(up) - model%numerics%stagsigma(up-1)) /    &
+           (model%numerics%stagsigma(up) - model%numerics%stagsigma(up-1))
+    end do
+    ! boundary conditions are identical on both grids, but temp starts at index 0
+    model%temper%tempunstag(1,:,:) = model%temper%temp(0,:,:)
+    model%temper%tempunstag(model%general%upn,:,:) = model%temper%temp(model%general%upn,:,:)
+
     !------------------------------------------------------------------------ 
     ! Halo updates
     !------------------------------------------------------------------------ 

libglissade/glissade_therm.F90

@@ -102,7 +102,8 @@ subroutine glissade_init_therm (temp_init,  is_restart,        &
                                   artm,                          &
                                   acab,                          &
                                   bheatflx,                      &
-                                  pmp_offset, temp)
+                                  pmp_offset,                    &
+                                  temp,       tempunstag)
 
     ! initialization subroutine for higher-order dycores, where temperature is defined at
     ! the midpoint of each layer plus the upper and lower surfaces
@@ -137,6 +138,12 @@ subroutine glissade_init_therm (temp_init,  is_restart,        &
          temp              ! ice temperature
                            ! intent(inout) because it might have been read already from an input file,
                            !  but otherwise is set in this subroutine
+
+    real(dp), dimension(:,:,:), intent(inout) ::  &
+         tempunstag        ! ice temperature on unstaggered grid
+                           ! may be used to initialize from an unstaggered glide temperature field
+                           ! intent(inout) because it might have been read from an input file,
+                           ! but otherwise it is set (diagnostically) in this subroutine
 
     ! Local variables
 
@@ -146,12 +153,6 @@ subroutine glissade_init_therm (temp_init,  is_restart,        &
 
     logical :: verbose_column    ! if true, then write diagnostic info for the column
 
-    !WHL - option to overwrite initial temperature in input file
-    !      Can be useful is working with an input file that includes initial temperatures,
-    !       but we want to set a linear temperature profile instead
-    !TODO - Make this a config option?
-    logical, parameter :: overwrite_input_temps = .false.
-
     ! Precompute some grid quantities used in the vertical temperature solve
 
     allocate(dups(upn+1,2))   !TODO - upn-1 instead?
@@ -171,61 +172,93 @@ subroutine glissade_init_therm (temp_init,  is_restart,        &
     up = upn-1
     dups(up,2) = 1.d0/((sigma(up+1) - sigma(up)) * (sigma(up+1) - stagsigma(up)) )
 
-    ! Check for a possible input error.  If the user supplies a file with the 'temp' field, which has
-    ! vertical dimension (1:upn), then the temperature in layers 1:upn may appear correct (though
-    ! staggered incorrectly), but the temperature in layer 0 will remain at an unphysical value.
-    ! Let the user know if this has happened.
-    !WHL - Nov. 2014 - I verified that the code aborts here if temp (rather than tempstag) is in the restart file.
-
-    if (minval(temp(0,:,:)) < (-1.d0*trpt) .and. minval(temp(1:upn,:,:)) > (-1.d0*trpt)) then
-       call write_log('Error, temperature field has been read incorrectly. Note that the '  &
-                   // 'Glissade dycore must be initialized with tempstag, not temp.', GM_FATAL)
-    endif
 
     !==== Initialize ice temperature.============
-    ! Six possibilities:
-    ! (1) Set ice temperature to 0 C everywhere in column (TEMP_INIT_ZERO).
-    ! (2) Set ice temperature to surface air temperature everywhere in column (TEMP_INIT_ARTM).
-    ! (3) Set up a linear temperature profile, with T = artm at the surface and T <= Tpmp
+    ! Following possibilities:
+    ! (0) Set ice temperature to 0 C everywhere in column (TEMP_INIT_ZERO).
+    ! (1) Set ice temperature to surface air temperature everywhere in column (TEMP_INIT_ARTM).
+    ! (2) Set up a linear temperature profile, with T = artm at the surface and T <= Tpmp
     !     at the bed (TEMP_INIT_LINEAR).
     !     A parameter (pmpt_offset) controls how far below Tpmp the initial bed temp is set.
-    ! (4) Set up a temperature profile based on advective-diffusive balance, with T = artm
+    ! (3) Set up a temperature profile based on advective-diffusive balance, with T = artm
     !     at the surface and dT/dz = -F_geo/k at the bed (TEMP_INIT_ADVECTIVE_DIFFUSIVE).
     !     The temperature at each level is capped at the value computed by method (3).
-    ! (5) Read ice temperature from an initial input file.
-    ! (6) Read ice temperature from a restart file.
+    ! (4) Read ice temperature from external file (TEMP_INIT_EXTERNAL).
+    !   (4a) If variable tempstag is present: Set ice temperature to tempstag from input file.
+    !   (4b) If variable tempunstag is present (and tempstag is not): interpolate tempunstag
+    !        to the staggered ice temperature needed by the Glissade dycore.
     !
-    ! The default is (2).
-    ! Method (4) may be optimal for reducing spinup time in the interior of large ice sheets.
-    ! If not restarting and the temperature field is present in the input file, we do (5).
-    ! If restarting, we always do (6).
-    ! If (5) or (6), then the temperature field should already have been read from a file,
-    !  and the rest of this subroutine will do nothing.
-    ! Otherwise, the initial temperature is controlled by model%options%temp_init,
-    !  which can be read from the config file.
+    ! The default is (1).
+    ! Methods (0-3) overwrite any temperature given in input files.
+    ! Method (3) may be optimal for reducing spinup time in the interior of large ice sheets.
+    ! Option (4) requires that temperature is present in the input file.
+
+    if (temp_init == TEMP_INIT_EXTERNAL) then
+
+       ! Temperature from external file
+
+       ! Check for a possible input error.  If the user supplies a file with the 'temp' field, which has
+       ! vertical dimension (1:upn), then the temperature in layers 1:upn may appear correct (though
+       ! staggered incorrectly), but the temperature in layer 0 will remain at an unphysical value.
+       ! Let the user know if this has happened.
+       !WHL - Nov. 2014 - I verified that the code aborts here if temp (rather than tempstag) is in the restart file.
+
+       if (minval(temp(0,:,:)) < (-1.d0*trpt) .and. minval(temp(1:upn,:,:)) > (-1.d0*trpt)) then
+          call write_log('Error, temperature field has been read incorrectly. Note that the '  &
+             // 'Glissade dycore must be initialized with tempstag, not temp.' &
+             // 'You can rename temp to tempunstag if you want it to be interpolated ' &
+             // 'to the vertically staggered tempstag (loss of detail).', GM_FATAL)
+       endif
 
-    if (is_restart == RESTART_TRUE) then
 
-       ! Temperature has already been initialized from a restart file. 
-       ! (Temperature is always a restart variable.)
+       if ( minval(temp) > (-1.0d0 * trpt)) then  ! temperature has been read from an input file
+                                                      ! Note: trpt = 273.15 K
 
-       call write_log('Initializing ice temperature from the restart file')
+          ! Temperature has already been initialized from a restart or input file.
+          ! We know this because the default initial temps of unphys_val
+          ! (a large negative number) have been overwritten.
 
-    !WHL - debug - option to overwrite input file
-!!    elseif ( minval(temp) > (-1.0d0 * trpt) ) then  ! temperature has been read from an input file
-    elseif ( minval(temp) > (-1.0d0 * trpt) .and. .not.overwrite_input_temps) then  ! temperature has been read from an input file
-                                                      ! Note: trpt = 273.15 K
+          ! Initialise tempunstag, the temperature on unstaggered vertical grid
+          ! (not used for calculations)
+          ! Use sigma weighted interpolation from temp to layer interfaces
+          do up = 2, upn-1
+            tempunstag(up,:,:) = temp(up-1,:,:) + (temp(up,:,:) - temp(up-1,:,:)) *   &
+                 (sigma(up) - stagsigma(up-1)) / (stagsigma(up) - stagsigma(up-1))
+          end do
+          ! boundary conditions are identical on both grids, but temp starts at index 0
+          tempunstag(1,:,:) = temp(0,:,:)
+          tempunstag(upn,:,:) = temp(upn,:,:)
 
-       ! Temperature has already been initialized from an input file.
-       ! We know this because the default initial temps of unphys_val (a large negative number) have been overwritten.
+          call write_log('Initializing ice temperature from a restart/input file')
+
+
+       elseif ( maxval(tempunstag(:,:,:)) > (-1.0d0 * trpt)) then
+
+          ! Test if any temperature in physical range
+          ! If yes, we have data from restart or input for unstaggered tempunstag that we use to initialise temp
+
+          call write_log('Initializing ice temperature by interpolating tempunstag from restart/input file')
+
+          ! Set temp to linear interpolation from tempunstag at layer interfaces
+          do up = 1, upn-1
+            temp(up,:,:) = (tempunstag(up,:,:) + tempunstag(up+1,:,:)) * 0.5d0
+          end do
+          ! boundary conditions are identical on both grids, but temp starts at index 0
+          temp(0,:,:) = tempunstag(1,:,:)
+          temp(upn,:,:) = tempunstag(upn,:,:)
+
+       else
+          ! fatal error, neither tempstag nor tempunstag are present in the input files
+          call write_log('Error, temp_init = 4 requires temperature variable tempstag or ' &
+             // 'tempunstag specified in the restart or input files. ', GM_FATAL)
+       endif
 
-       call write_log('Initializing ice temperature from an input file')
 
     else   ! not reading temperature from restart or input file
-           ! initialize it here based on temp_init
+           ! initialize it here based on temp_init = (0-3)
 
        ! initialize T = 0 C everywhere
-       temp(:,:,:) = 0.0d0                                              
+       temp(:,:,:) = 0.0d0
 
        ! set temperature in each column based on the value of temp_init