Merge pull request #61 from ESCOMP/lipscomb/basal_physics4

Basal hydrology scheme, new basal friction and inversion options

cism_driver/cism_front_end.F90

@@ -226,7 +226,7 @@ subroutine cism_init_dycore(model)
 
   ! --- Output the initial state -------------
 
-  if (model%options%is_restart == RESTART_FALSE) then
+  if (model%options%is_restart == RESTART_FALSE .or. model%options%forcewrite_restart) then
      call t_startf('initial_io_writeall')
      call glide_io_writeall(model, model, time=time)          ! MJH The optional time argument needs to be supplied 
                                                               !     since we have not yet set model%numerics%time

cism_driver/eismint_forcing.F90

@@ -273,8 +273,12 @@ subroutine eismint_printconfig(eismint_climate)
        call write_log(message)
     end select
 
-    if ( (eismint_climate%eismint_type > 0) .and. (tasks > 1) ) then
-       call write_log('EISMINT tests are not supported for more than one processor', GM_FATAL)
+    ! Note: The EISMINT-1 tests might work with a parallel Glissade solver, but this has not been tested.
+    !       Only EISMINT-2 has been tested in parallel as of August 2021.
+    if (eismint_climate%eismint_type == 1 .or. eismint_climate%eismint_type == 2) then
+       if (tasks > 1) then
+          call write_log('EISMINT-1 tests are not supported for more than one processor', GM_WARNING)
+       endif
     end if
 
     call write_log('')
@@ -285,27 +289,27 @@ subroutine eismint_massbalance(eismint_climate,model,time)
 
     ! calculate eismint mass balance
 
-!TODO - Remove acc0
-
     use glimmer_global, only : dp
     use glide_types
     use glimmer_paramets, only : len0, acc0, scyr
     use glimmer_physcon, only : pi
     use glimmer_scales, only : scale_acab
+    use cism_parallel, only: parallel_globalindex
+
     implicit none
 
     type(eismint_climate_type) :: eismint_climate         ! structure holding climate info
     type(glide_global_type) :: model        ! model instance
     real(dp), intent(in) :: time            ! current time
 
-    !WHL - Changed 'periodic_bc' to 'periodic' to avoid a name conflict with parallel modules
     ! local variables
-    integer  :: ns,ew
+    integer  :: ew, ns, ew_global, ns_global
     real(dp) :: dist, ewct, nsct, grid, rel
-    real(dp) :: periodic = 1.d0  !TODO - Make this an integer?
+    real(dp) :: periodic = 1.d0
+
+    ewct = (model%parallel%global_ewn + 1.d0) / 2.d0
+    nsct = (model%parallel%global_nsn + 1.d0) / 2.d0
 
-    ewct = (real(model%general%ewn,dp) + 1.d0) / 2.d0
-    nsct = (real(model%general%nsn,dp) + 1.d0) / 2.d0
     grid = real(model%numerics%dew,dp) * len0
 
     if (model%options%periodic_ew) then
@@ -322,7 +326,6 @@ subroutine eismint_massbalance(eismint_climate,model,time)
        if (eismint_climate%period .ne. 0.d0) then
           model%climate%acab(:,:) = model%climate%acab(:,:) + eismint_climate%mb_amplitude * &
                sin(2.d0*pi*time/eismint_climate%period)/ (acc0 * scyr)
-!          model%climate%acab(:,:) = model%climate%acab(:,:) + climate%mb_amplitude * sin(2.d0*pi*time/climate%period) / scale_acab
        end if
 
     case(2)
@@ -335,7 +338,8 @@ subroutine eismint_massbalance(eismint_climate,model,time)
 
        do ns = 1,model%general%nsn
           do ew = 1,model%general%ewn
-             dist = grid * sqrt(periodic*(real(ew,kind=dp) - ewct)**2 + (real(ns,kind=dp) - nsct)**2)
+             call parallel_globalindex(ew, ns, ew_global, ns_global, model%parallel)
+             dist = grid * sqrt((real(ew_global,kind=dp) - ewct)**2 + (real(ns_global,kind=dp) - nsct)**2)
              model%climate%acab(ew,ns) = min(eismint_climate%nmsb(1), eismint_climate%nmsb(2) * (rel - dist))
           end do
        end do
@@ -346,7 +350,8 @@ subroutine eismint_massbalance(eismint_climate,model,time)
 
        do ns = 1,model%general%nsn
           do ew = 1,model%general%ewn
-             dist = grid * sqrt(periodic*(real(ew,kind=dp) - ewct)**2 + (real(ns,kind=dp) - nsct)**2)
+             call parallel_globalindex(ew, ns, ew_global, ns_global, model%parallel)
+             dist = grid * sqrt((real(ew_global,kind=dp) - ewct)**2 + (real(ns_global,kind=dp) - nsct)**2)
              model%climate%acab(ew,ns) = min(eismint_climate%nmsb(1), eismint_climate%nmsb(2) * (rel - dist))
           end do
        end do
@@ -367,19 +372,21 @@ subroutine eismint_surftemp(eismint_climate,model,time)
     use glimmer_global, only: dp
     use glimmer_paramets, only : len0
     use glimmer_physcon, only : pi
+    use cism_parallel, only: parallel_globalindex
+
     implicit none
 
     type(eismint_climate_type) :: eismint_climate         ! structure holding climate info
     type(glide_global_type) :: model        ! model instance
     real(dp), intent(in) :: time            ! current time
 
     ! local variables
-    integer  :: ns,ew
+    integer  :: ew, ns, ew_global, ns_global
     real(dp) :: dist, ewct, nsct, grid
     real(dp) :: periodic = 1.d0
 
-    ewct = (real(model%general%ewn,dp)+1.d0) / 2.d0
-    nsct = (real(model%general%nsn,dp)+1.d0) / 2.d0
+    ewct = (model%parallel%global_ewn + 1.d0) / 2.d0
+    nsct = (model%parallel%global_nsn + 1.d0) / 2.d0
     grid = real(model%numerics%dew,dp) * len0
 
     if (model%options%periodic_ew) then
@@ -394,7 +401,8 @@ subroutine eismint_surftemp(eismint_climate,model,time)
        ! EISMINT-1 fixed margin
        do ns = 1,model%general%nsn
           do ew = 1,model%general%ewn
-             dist = grid * max(periodic*abs(real(ew,kind=dp) - ewct),abs(real(ns,kind=dp) - nsct))*1d-3
+             call parallel_globalindex(ew, ns, ew_global, ns_global, model%parallel)
+             dist = grid * max(abs(real(ew_global,kind=dp) - ewct),abs(real(ns_global,kind=dp) - nsct))*1d-3
              model%climate%artm(ew,ns) = eismint_climate%airt(1) + eismint_climate%airt(2) * dist*dist*dist
           end do
        end do
@@ -413,7 +421,8 @@ subroutine eismint_surftemp(eismint_climate,model,time)
        ! EISMINT-2
        do ns = 1,model%general%nsn
           do ew = 1,model%general%ewn
-             dist = grid * sqrt(periodic*(real(ew,kind=dp) - ewct)**2 + (real(ns,kind=dp) - nsct)**2)
+             call parallel_globalindex(ew, ns, ew_global, ns_global, model%parallel)
+             dist = grid * sqrt((real(ew_global,kind=dp) - ewct)**2 + (real(ns_global,kind=dp) - nsct)**2)
              model%climate%artm(ew,ns) = eismint_climate%airt(1)+eismint_climate%airt(2) * dist
           end do
        end do
@@ -428,6 +437,10 @@ end subroutine eismint_surftemp
   !which_call - eismint_surftemp(0)/eismint_massbalance(1)/both(2)
   !which_test - test f(0)/test g(1)/exact(2)
 
+  !Note: This subroutine is called only for eismint_climate%eismint_type = 4,
+  !      which is not currently supported.
+  !      It would need some changes to work for parallel runs; see eismint_surftemp above.
+
   subroutine exact_surfmass(eismint_climate,model,time,which_call,which_test)
 
     use glide_types
@@ -450,7 +463,7 @@ subroutine exact_surfmass(eismint_climate,model,time,which_call,which_test)
 
     !TODO - Change which_call to an integer?
     !       Modify for Glissade? (dissip has smaller vertical dimension)
-    if (which_call .eq. 0.d0 .or. which_call .eq. 2.d0) then
+    if (which_call == 0 .or. which_call == 2) then
 
         !point by point call to the function 
         do ns = 1,model%general%nsn
@@ -464,9 +477,9 @@ subroutine exact_surfmass(eismint_climate,model,time,which_call,which_test)
                     if(r>0.d0 .and. r<L) then
                         if (which_test .eq. 0.d0) then
                             call testF(r,z,H,TT,U,w,Sig,M,Sigc)
-                        else if (which_test .eq. 1.d0) then
+                        else if (which_test == 1) then
                             call testG(time,r,z,H,TT,U,w,Sig,M,Sigc)
-                        else if (which_test .eq. 2.d0) then
+                        else if (which_test == 2) then
                             !H_0 = H0
                             H_0 = model%geometry%thck(center,center)
                             !TT = 0.d0
@@ -482,7 +495,7 @@ subroutine exact_surfmass(eismint_climate,model,time,which_call,which_test)
             end do
         end do
 
-    else if (which_call .eq. 1.d0 .or. which_call .eq. 2.d0) then
+    else if (which_call == 1 .or. which_call == 2) then
 
         do ns = 1,model%general%nsn
             do ew = 1,model%general%ewn
@@ -492,11 +505,11 @@ subroutine exact_surfmass(eismint_climate,model,time,which_call,which_test)
                 z = model%geometry%thck(ew,ns)
                 !the function only returns values within the radius
                 if(r>0.d0 .and. r<L) then
-                    if (which_test .eq. 0.d0) then
+                    if (which_test == 0) then
                         call testF(r,z,H,TT,U,w,Sig,M,Sigc)
-                    else if (which_test .eq. 1.d0) then
+                    else if (which_test == 1) then
                         call testG(time,r,z,H,TT,U,w,Sig,M,Sigc)
-                    else if (which_test .eq. 2.d0) then
+                    else if (which_test == 2) then
                         H_0 = H0 !H_0 = model%geometry%thck(center,center)
                         TT = 0.d0 !TT = model%temper%dissip(lev,ew,ns)
                         call model_exact(time,r,z,model%geometry%thck(ew,ns),H_0,TT,U,w,Sig,M,Sigc)

libglide/glide.F90

@@ -85,8 +85,6 @@ subroutine glide_config(model,config,fileunit)
     type(ConfigSection), pointer :: ncconfig
     integer :: unit
 
-    integer :: k   !WHL - debug
-
     unit = 99
     if (present(fileunit)) then
        unit = fileunit
@@ -594,8 +592,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice)
           call calcbwat(model, &
                         model%options%whichbwat, &
                         model%basal_melt%bmlt, &
-                        model%temper%bwat, &
-                        model%temper%bwatflx, &
+                        model%basal_hydro%bwat, &
+                        model%basal_hydro%bwatflx, &
                         model%geometry%thck, &
                         model%geometry%topg, &
                         model%temper%temp(model%general%upn,:,:), &
@@ -604,8 +602,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice)
 
 
           ! This call is redundant for now, but is needed if the call to calcbwat is removed
-          call stagvarb(model%temper%bwat, &
-                        model%temper%stagbwat ,&
+          call stagvarb(model%basal_hydro%bwat, &
+                        model%basal_hydro%stagbwat ,&
                         model%general%ewn, &
                         model%general%nsn)
 
@@ -859,8 +857,8 @@ subroutine glide_tstep_p1(model,time)
        call calcbwat(model, &
                      model%options%whichbwat, &
                      model%basal_melt%bmlt, &
-                     model%temper%bwat, &
-                     model%temper%bwatflx, &
+                     model%basal_hydro%bwat, &
+                     model%basal_hydro%bwatflx, &
                      model%geometry%thck, &
                      model%geometry%topg, &
                      model%temper%temp(model%general%upn,:,:), &
@@ -911,9 +909,6 @@ subroutine glide_tstep_p2(model)
 
     type(glide_global_type), intent(inout) :: model    ! model instance
 
-    !debug
-    integer :: j
-
     ! ------------------------------------------------------------------------ 
     ! Calculate flow evolution by various different methods
     ! ------------------------------------------------------------------------ 

libglide/glide_bwater.F90

@@ -195,8 +195,8 @@ subroutine calcbwat(model, which, bmlt, bwat, bwatflx, thck, topg, btem, floater
     end select
 
     ! now also calculate basal water in velocity (staggered) coord system
-    call stagvarb(model%temper%bwat, &
-                  model%temper%stagbwat ,&
+    call stagvarb(model%basal_hydro%bwat, &
+                  model%basal_hydro%stagbwat ,&
                   model%general%ewn, &
                   model%general%nsn)
 
@@ -429,7 +429,7 @@ subroutine pressure_wphi(thck,topg,N,wphi,thicklim,floater)
   !> Compute the pressure wphi at the base of the ice sheet according to
   !> ice overburden plus bed height minus effective pressure.
   !>
-  !> whpi/(rhow*g) = topg + bwat * rhoi / rhow * thick - N / (rhow * g)
+  !> wphi/(rhow*g) = topg + bwat * rhoi / rhow * thick - N / (rhow * g)
 
     use glimmer_physcon, only : rhoi,rhow,grav
     implicit none

libglide/glide_diagnostics.F90

@@ -200,9 +200,12 @@ subroutine glide_write_diag (model,       time)
          max_thck, max_thck_global,     &    ! max ice thickness (m)
          max_temp, max_temp_global,     &    ! max ice temperature (deg C)
          min_temp, min_temp_global,     &    ! min ice temperature (deg C)
-         max_bmlt, max_bmlt_global,     &    ! max basal melt rate (m/yr)
-         max_spd_sfc, max_spd_sfc_global,   &    ! max surface ice speed (m/yr)
-         max_spd_bas, max_spd_bas_global,   &    ! max basal ice speed (m/yr)
+         max_bmlt,                      &    ! max basal melt rate (m/yr)
+         max_bmlt_global,               &
+         max_bmlt_ground,               &    ! max basal melt rate, grounded ice (m/yr)
+         max_bmlt_ground_global,        &
+         max_spd_sfc, max_spd_sfc_global, &  ! max surface ice speed (m/yr)
+         max_spd_bas, max_spd_bas_global, &  ! max basal ice speed (m/yr)
          spd,                           &    ! speed
          thck_diag, usrf_diag,          &    ! local column diagnostics
          topg_diag, relx_diag,          &    
@@ -768,7 +771,8 @@ subroutine glide_write_diag (model,       time)
                     min_temp_global, imin_global, jmin_global, kmin_global
     call write_log(trim(message), type = GM_DIAGNOSTIC)
 
-    ! max basal melt rate
+    ! max applied basal melt rate
+    ! Usually, this will be for floating ice, if floating ice is present
     imax = 0
     jmax = 0
     max_bmlt = unphys_val
@@ -791,11 +795,39 @@ subroutine glide_write_diag (model,       time)
     ! Write to diagnostics only if nonzero
 
     if (abs(max_bmlt_global*thk0*scyr/tim0) > eps) then
-       write(message,'(a25,f24.16,2i6)') 'Max bmlt (m/yr), i, j    ',   &
+       write(message,'(a25,f24.16,2i6)') 'Max bmlt (m/y), i, j     ',   &
             max_bmlt_global*thk0*scyr/tim0, imax_global, jmax_global
        call write_log(trim(message), type = GM_DIAGNOSTIC)
     endif
 
+    ! max basal melt rate for grounded ice
+    imax = 0
+    jmax = 0
+    max_bmlt_ground = unphys_val
+
+    do j = lhalo+1, nsn-uhalo
+       do i = lhalo+1, ewn-uhalo
+          if (model%basal_melt%bmlt_ground(i,j) > max_bmlt_ground) then
+             max_bmlt_ground = model%basal_melt%bmlt_ground(i,j)
+             imax = i
+             jmax = j
+          endif
+       enddo
+    enddo
+
+    call parallel_reduce_maxloc(xin=max_bmlt_ground, xout=max_bmlt_ground_global, xprocout=procnum)
+    call parallel_globalindex(imax, jmax, imax_global, jmax_global, parallel)
+    call broadcast(imax_global, proc = procnum)
+    call broadcast(jmax_global, proc = procnum)
+
+    ! Write to diagnostics only if nonzero
+
+    if (abs(max_bmlt_global*thk0*scyr/tim0) > eps) then
+       write(message,'(a25,f24.16,2i6)') 'Max bmlt grnd (m/y), i, j',   &
+            max_bmlt_ground_global*thk0*scyr/tim0, imax_global, jmax_global
+       call write_log(trim(message), type = GM_DIAGNOSTIC)
+    endif
+
     ! max surface speed
     imax = 0
     jmax = 0
@@ -818,7 +850,7 @@ subroutine glide_write_diag (model,       time)
     call broadcast(imax_global, proc = procnum)
     call broadcast(jmax_global, proc = procnum)
 
-    write(message,'(a25,f24.16,2i6)') 'Max sfc spd (m/yr), i, j ',   &
+    write(message,'(a25,f24.16,2i6)') 'Max sfc spd (m/y), i, j  ',   &
                     max_spd_sfc_global*vel0*scyr, imax_global, jmax_global
     call write_log(trim(message), type = GM_DIAGNOSTIC)
 
@@ -843,7 +875,7 @@ subroutine glide_write_diag (model,       time)
     call parallel_globalindex(imax, jmax, imax_global, jmax_global, parallel)
     call broadcast(imax_global, proc = procnum)
     call broadcast(jmax_global, proc = procnum)
-    write(message,'(a25,f24.16,2i6)') 'Max base spd (m/yr), i, j',   &
+    write(message,'(a25,f24.16,2i6)') 'Max base spd (m/y), i, j ',   &
                     max_spd_bas_global*vel0*scyr, imax_global, jmax_global
     call write_log(trim(message), type = GM_DIAGNOSTIC)
 
@@ -884,7 +916,7 @@ subroutine glide_write_diag (model,       time)
           artm_diag = model%climate%artm_corrected(i,j)  ! artm_corrected = artm + artm_anomaly
           acab_diag = model%climate%acab_applied(i,j) * thk0*scyr/tim0
           bmlt_diag = model%basal_melt%bmlt_applied(i,j) * thk0*scyr/tim0
-          bwat_diag = model%temper%bwat(i,j) * thk0
+          bwat_diag = model%basal_hydro%bwat(i,j) * thk0
           bheatflx_diag = model%temper%bheatflx(i,j)
 
           temp_diag(:) = model%temper%temp(1:upn,i,j)