Merge pull request #36 from ESCOMP/lipscomb/ais_coupled

Changes to support CESM ocean coupling

libglad/glad_main.F90

@@ -125,10 +125,10 @@ module glad_main
   ! When coupled to CESM, Glad receives several fields from the coupler on the ice sheet grid:
   !   qsmb = surface mass balance (kg/m^2/s)
   !   tsfc = surface ground temperature (deg C)
-  !   salinity1..7 = ocean salinity at levels 0,10,19,26,30,33,35 (g/kg)
-  !   tocn1..7     = ocean temperatures at levels 0,10,19,26,30,33,35 (deg K)
+  !   salinity = ocean salinity at one or more ocean levels (g/kg)
+  !   tocn = ocean temperature at one or more levels (deg K)
   ! Both qsmb and tsfc are computed in the CESM land model.
-  ! Both set of fields salinity1..7 and tocan1..7 are computed in POP.
+  ! Both salinity and tocn are computed in the CESM ocean model.
   ! Seven fields are returned to CESM on the ice sheet grid:
   !   ice_covered = whether a grid cell is ice-covered [0,1]
   !   topo = surface elevation (m)
@@ -584,13 +584,17 @@ subroutine glad_gcm(params,         instance_index, time,                 &
 
     ! Main Glad subroutine for GCM coupling.
     !
-    ! It does all necessary temporal averaging, 
-    ! and calls the dynamic ice sheet model when required. 
+    ! It does all temporal averaging and calls the dynamic ice sheet model when required.
     !
     ! Input fields should be taken as means over the period since the last call.
     ! See the user documentation for more information.
     !
-    ! Input fields are assumed to NOT have halo cells
+    ! Input fields are assumed to NOT have halo cells.
+    !
+    ! Glad accumulates and averages ocean thermal_forcing based on salinity and tocn,
+    ! but this thermal forcing is used by CISM only if model%options%ocean_data_domain = 2.
+    ! Otherwise, CISM will compute thermal_forcing internally (ocean_data_domain = 0)
+    ! or read it from an external file (ocean_data_domain = 1).
 
     use glimmer_utils
     use glad_timestep, only : glad_i_tstep_gcm
@@ -711,13 +715,14 @@ subroutine glad_gcm(params,         instance_index, time,                 &
                thermal_forcing_haloed(k,:,:))
        enddo
 
+
        if (verbose_glad .and. this_rank == rtest) then
           i = itest
           j = jtest
           print*, 'r, i, j =', this_rank, i, j
           print*, 'k, zocn, temperature, salinity, thermal forcing:'
           do k = 1, nzocn
-             write(6,'(i4, 4f10.3)') k, zocn(k), &
+             write(6,'(i4, 4f11.3)') k, zocn(k), &
                   tocn_haloed(k,i,j), salinity_haloed(k,i,j), thermal_forcing_haloed(k,i,j)
           enddo
        endif
@@ -799,7 +804,7 @@ subroutine glad_gcm(params,         instance_index, time,                 &
              j = jtest
              print*, 'Before calling glad_i_tstep_gcm, k, zocn, average thermal forcing:'
              do k = 1, nzocn
-                write(6,'(i4, 2f10.3)') k, zocn(k), params%instances(instance_index)%thermal_forcing(k,i,j)
+                write(6,'(i4, 2f11.3)') k, zocn(k), params%instances(instance_index)%thermal_forcing(k,i,j)
              enddo
           endif
 

libglad/glad_mbal_coupling.F90

@@ -104,9 +104,7 @@ end subroutine glad_mbc_init
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-  subroutine glad_accumulate_input_gcm(params, time, acab, artm, thermal_forcing) !, &
-!                                       thermal_forcing2, thermal_forcing3, thermal_forcing4,  &
-!                                       thermal_forcing5, thermal_forcing6, thermal_forcing7)
+  subroutine glad_accumulate_input_gcm(params, time, acab, artm, thermal_forcing)
 
     ! In glint, this was done in glint_downscale.F90
 
@@ -115,7 +113,7 @@ subroutine glad_accumulate_input_gcm(params, time, acab, artm, thermal_forcing)
 
     real(dp),dimension(:,:),intent(in) :: acab   ! Surface mass balance (m)
     real(dp),dimension(:,:),intent(in) :: artm   ! Mean air temperature (degC)
-    real(dp),dimension(:,:,:),intent(in) :: thermal_forcing   ! Mean thermal_forcing at level 0 (degK)
+    real(dp),dimension(:,:,:),intent(in) :: thermal_forcing   ! Ocean thermal_forcing (degK)
 
     ! Things to do the first time
 
@@ -161,7 +159,7 @@ subroutine glad_average_input_gcm(params, dt, acab, artm, thermal_forcing)
     integer,                intent(in)   :: dt     !> mbal accumulation time (hours)
     real(dp),dimension(:,:),intent(out)  :: artm   !> Mean air temperature (degC)
     real(dp),dimension(:,:),intent(out)  :: acab   !> Mass-balance (m/yr)
-    real(dp),dimension(:,:,:),intent(out)  :: thermal_forcing   ! Mean thermal_forcing at level 0 (degK)
+    real(dp),dimension(:,:,:),intent(out)  :: thermal_forcing   ! Ocean thermal_forcing (degK)
 
     if (.not. params%new_accum) then
        params%artm_save = params%artm_save / real(params%av_count,dp)

libglad/glad_timestep.F90

@@ -66,6 +66,7 @@ subroutine glad_i_tstep_gcm(time,            instance,       &
     use glide
     use glissade
     use glide_io
+    use glide_types
     use glad_mbal_coupling, only : glad_accumulate_input_gcm, glad_average_input_gcm
     use glad_io
     use glad_mbal_io
@@ -200,16 +201,21 @@ subroutine glad_i_tstep_gcm(time,            instance,       &
           call glide_set_acab(instance%model, instance%acab * rhow/rhoi)
           call glide_set_artm(instance%model, instance%artm)
 
+          ! Note: CISM has several thermal forcing options, as determined by ocean_data_domain:
+          !        compute internally, read from external file, or receive from Glad.
+          !       Only if receiving from Glad do we set instance%model%ocean_data%thermal_forcing here.
+          !
           ! Note: The ocean thermal forcing is reset only on the first ice dynamics timestep within this loop.
-          ! The reason is that CISM has the option of extrapolating thermal forcing from open ocean
-          !  (i.e., the overlap region between the ocean and ice sheet domains) into sub-shelf cavities.
-          ! It is more efficient to do this extrapolation only once within the mass_balance timestep
-          !  (with minor subsequent adjustments if CISM ice shelves retreat) than to extrapolate
-          !  from the open ocean every ice dynamics time step.
+          !       The reason is that CISM has the option of extrapolating thermal forcing from open ocean
+          !        (i.e., the overlap region between the ocean and ice sheet domains) into sub-shelf cavities.
+          !       It is more efficient to do this extrapolation only once within the mass_balance timestep
+          !        (with minor subsequent adjustments if CISM ice shelves retreat) than to extrapolate
+          !        from the open ocean every ice dynamics time step.
 
           if (iter == 1) then
-             if ( associated(instance%model%ocean_data%thermal_forcing) ) then
-                ! GL: At this point, glide_set does not work for 3D variables.
+             if ( associated(instance%model%ocean_data%thermal_forcing) .and. &
+                  instance%model%options%ocean_data_domain == OCEAN_DATA_GLAD) then
+                ! GRL: At this point, glide_set does not work for 3D variables.
                 instance%model%ocean_data%thermal_forcing = instance%thermal_forcing
              endif
           endif

libglide/glide_setup.F90

@@ -722,6 +722,7 @@ subroutine handle_options(section, model)
     call GetValue(section,'bmlt_float_thermal_forcing_param',model%options%bmlt_float_thermal_forcing_param)
     call GetValue(section,'bmlt_float_ismip6_magnitude',model%options%bmlt_float_ismip6_magnitude)
     call GetValue(section,'ocean_data_domain',model%options%ocean_data_domain)
+    call GetValue(section,'ocean_data_extrapolate',model%options%ocean_data_extrapolate)
     call GetValue(section,'enable_bmlt_anomaly',model%options%enable_bmlt_anomaly)
     call GetValue(section,'basal_mass_balance',model%options%basal_mbal)
     call GetValue(section,'smb_input',model%options%smb_input)
@@ -936,9 +937,13 @@ subroutine print_options(model)
          'highest forcing magnitude '  /)
 
     character(len=*), dimension(0:2), parameter :: ocean_data_domain = (/ &
-         'data from external ocean model domain; extrapolate to cavities', &
-         'ocean data already extrapolated to ice shelf cavities         ', &
-         'apply data in CISM ice-free ocean; extrapolate to cavities    ' /)
+         'ocean data computed internally by CISM', &
+         'ocean data read from external file    ', &
+         'ocean data from coupler via Glad      ' /)
+
+    character(len=*), dimension(0:1), parameter :: ocean_data_extrapolate = (/ &
+         'ocean data not extrapolated to cavities', &
+         'ocean data extrapolated to cavities    ' /)
 
     character(len=*), dimension(0:1), parameter :: smb_input = (/ &
          'SMB input in units of m/yr ice  ', &
@@ -1522,6 +1527,9 @@ subroutine print_options(model)
        write(message,*) 'ocean data domain       : ', model%options%ocean_data_domain, &
             ocean_data_domain(model%options%ocean_data_domain)
        call write_log(message)
+       write(message,*) 'ocean data extrapolate  : ', model%options%ocean_data_extrapolate, &
+            ocean_data_extrapolate(model%options%ocean_data_extrapolate)
+       call write_log(message)
     endif
 
     if (model%options%basal_mbal < 0 .or. model%options%basal_mbal >= size(b_mbal)) then

libglide/glide_types.F90

@@ -137,9 +137,12 @@ module glide_types
   integer, parameter :: BMLT_FLOAT_ISMIP6_MEDIAN = 1
   integer, parameter :: BMLT_FLOAT_ISMIP6_PCT95 = 2
 
-  integer, parameter :: DATA_OCEAN_ONLY = 0
-  integer, parameter :: DATA_OCEAN_ICE = 1
-  integer, parameter :: DATA_CISM_OCEAN_MASK = 2
+  integer, parameter :: OCEAN_DATA_INTERNAL = 0
+  integer, parameter :: OCEAN_DATA_EXTERNAL = 1
+  integer, parameter :: OCEAN_DATA_GLAD = 2
+
+  integer, parameter :: OCEAN_DATA_EXTRAPOLATE_FALSE = 0
+  integer, parameter :: OCEAN_DATA_EXTRAPOLATE_TRUE = 1
 
   integer, parameter :: BASAL_MBAL_NO_CONTINUITY = 0
   integer, parameter :: BASAL_MBAL_CONTINUITY = 1
@@ -522,12 +525,18 @@ module glide_types
     !> \item[2] High level of forcing (e.g., pct95)
     !> \end{description}
 
-    integer :: ocean_data_domain = 0
+    integer :: ocean_data_domain = 1
+
+    !> \begin{description}
+    !> \item[0] ocean data computed internally by CISM
+    !> \item[1] ocean data read from external file
+    !> \item[2] ocean data received from coupler via Glad
+    !> \end{description}
 
+    integer :: ocean_data_extrapolate = 0
     !> \begin{description}
-    !> \item[0] ocean data on ocean domain only; extrapolate data to shelf cavities
-    !> \item[1] ocean data available everywhere; already extrapolated to shelf cavities
-    !> \item[2] ocean data applied where CISM has ice-free ocean; extrapolated to shelf cavities
+    !> \item[0] ocean data not extrapolated to shelf cavities
+    !> \item[1] ocean data extrapolated to shelf cavities
     !> \end{description}
 
     logical :: enable_bmlt_anomaly = .false.

libglimmer/parallel_mpi.F90

@@ -248,6 +248,11 @@ module cism_parallel
      module procedure distributed_scatter_var_col_real8_2d
   end interface
 
+  interface parallel_boundary_value
+     module procedure parallel_boundary_value_real8_2d
+     module procedure parallel_boundary_value_real8_3d
+  end interface parallel_boundary_value
+
   interface parallel_convert_haloed_to_nonhaloed
      module procedure parallel_convert_haloed_to_nonhaloed_real4_2d
      module procedure parallel_convert_haloed_to_nonhaloed_real8_2d
@@ -5020,7 +5025,7 @@ end subroutine parallel_barrier
 
 !=======================================================================
 
-  function parallel_boundary(ew, ewn, ns, nsn, parallel)
+  function parallel_boundary(ew, ns, parallel)
 
     implicit none
     logical :: parallel_boundary
@@ -5029,6 +5034,8 @@ function parallel_boundary(ew, ewn, ns, nsn, parallel)
 
     ! begin
     associate(  &
+         local_ewn   => parallel%local_ewn,    &
+         local_nsn   => parallel%local_nsn,    &
          global_ewn  => parallel%global_ewn,   &
          global_nsn  => parallel%global_nsn,   &
          ewlb        => parallel%ewlb,         &
@@ -5038,14 +5045,83 @@ function parallel_boundary(ew, ewn, ns, nsn, parallel)
          )
 
     parallel_boundary = (ewlb<1 .and. ew==1+lhalo) .or. &
-             (ewub>global_ewn .and. ew==ewn-uhalo) .or. &
+               (ewub>global_ewn .and. ew==local_ewn-uhalo) .or. &
                         (nslb<1 .and. ns==1+lhalo) .or. &
-             (nsub>global_nsn .and. ns==nsn-uhalo)
+               (nsub>global_nsn .and. ns==local_nsn-uhalo)
 
     end associate
 
   end function parallel_boundary
 
+!=======================================================================
+
+  ! subroutines belonging to the parallel_boundary_value interface
+
+  subroutine parallel_boundary_value_real8_2d(&
+       field,          &
+       boundary_value, &
+       parallel)
+
+    ! Insert a specified value into cells on the global boundary.
+    ! Typically called with value = 0.0 or an special value.
+
+    real(dp), dimension(:,:), intent(inout) :: field
+    real(dp), intent(in) :: boundary_value
+    type(parallel_type) :: parallel
+
+    integer :: ew, ns
+
+    ! begin
+    associate(  &
+         local_ewn   => parallel%local_ewn,    &
+         local_nsn   => parallel%local_nsn     &
+         )
+
+    do ns = 1, local_nsn
+       do ew = 1, local_ewn
+          if (parallel_boundary(ew, ns, parallel)) then
+             field(ew,ns) = boundary_value
+          endif
+       enddo
+    enddo
+
+    end associate
+
+  end subroutine parallel_boundary_value_real8_2d
+
+
+  subroutine parallel_boundary_value_real8_3d(&
+       field,          &
+       boundary_value, &
+       parallel)
+
+    ! Insert a specified value into cells on the global boundary.
+    ! Typically called with value = 0.0 or an special value.
+
+    real(dp), dimension(:,:,:), intent(inout) :: field
+    real(dp), intent(in) :: boundary_value
+    type(parallel_type) :: parallel
+
+    integer :: ew, ns
+
+    ! begin
+    associate(  &
+         local_ewn   => parallel%local_ewn,    &
+         local_nsn   => parallel%local_nsn     &
+         )
+
+    do ns = 1, local_nsn
+       do ew = 1, local_ewn
+          if (parallel_boundary(ew, ns, parallel)) then
+             field(:,ew,ns) = boundary_value
+          endif
+       enddo
+    enddo
+
+    end associate
+
+  end subroutine parallel_boundary_value_real8_3d
+
 !=======================================================================
 
   function parallel_close(ncid)