Lipscomb/hma glaciers4

cism_driver/cism_front_end.F90

@@ -137,8 +137,12 @@ subroutine cism_init_dycore(model)
   call eismint_surftemp(model%eismint_climate,model,time)
 
   ! read forcing time slice if needed - this will overwrite values from IC file if there is a conflict.
+  ! Note: The first 'model' is passed to the argument 'data', which is filled by calling glide_read.
   call glide_read_forcing(model, model)
 
+  ! Optionally, read all the time slices at once from selected forcing files.
+  call glide_read_forcing_once(model, model)
+
   call spinup_lithot(model)
 
   call t_stopf('initialization')
@@ -225,8 +229,10 @@ subroutine cism_init_dycore(model)
 
 
   ! --- Output the initial state -------------
+  ! Note: For a standard restart, the initial state is not output, because this state
+  !       should already have been written to the output file when the previous run ended.
 
-  if (model%options%is_restart == RESTART_FALSE .or. model%options%forcewrite_restart) then
+  if (model%options%is_restart == NO_RESTART .or. model%options%is_restart == HYBRID_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
@@ -283,16 +289,16 @@ subroutine cism_run_dycore(model)
     do while(time + time_eps < model%numerics%tend)
 
       !!! SFP moved block of code for applying time dependent forcing read in from netCDF here,
-      !!! as opposed to at the end of the time step (commented it out in it's original location for now)
+      !!! as opposed to at the end of the time step (commented it out in its original location for now)
       !!! This is a short-term fix. See additional discussion as part of issue #19 (in cism-piscees github repo).
 
       ! Forcing from a 'forcing' data file - will read time slice if needed
-      ! Note: Forcing is read from the appropriate time slice after every dynamic time step.
-      !       This is not strictly necessary if there are multiple time steps per forcing time slice.
-      !       We would need additional logic if we wanted to read a new time slice only when needed
-      !        to replace the current data.  TODO: Add this logic?
       call t_startf('read_forcing')
       call glide_read_forcing(model, model)
+
+      ! If any forcing data have been read once into Fortran arrays at initialization,
+      !  simply copy the data based on the current forcing time.
+      call glide_retrieve_forcing(model, model)
       call t_stopf('read_forcing')
 
       ! Increment time step

libglide/glide.F90

@@ -218,6 +218,11 @@ subroutine glide_initialise(model)
     ! read first time slice
     call glide_io_readall(model,model)
 
+    ! Compute grid cell areas
+    ! Note: cell_area is used for diagnostics only. It is set to dew*dns by default but can be corrected below.
+    !       For the purposes of CISM dynamics, all grid cells are rectangles of dimension dew*dns.
+    model%geometry%cell_area = model%numerics%dew*model%numerics%dns
+
     ! Compute area scale factors for stereographic map projection.
     ! This should be done after reading the input file, in case the input file contains mapping info.
     ! Note: Not yet enabled for other map projections.
@@ -231,6 +236,14 @@ subroutine glide_initialise(model)
                                       model%numerics%dew*len0, &
                                       model%numerics%dns*len0)
 
+       ! Given the stereographic area correction factors, correct the diagnostic grid cell areas.
+       ! Note: area_factor is actually a length correction factor k; must divide by k^2 to adjust areas.
+       ! TODO: Change the name of area_factor
+       where (model%projection%stere%area_factor > 0.0d0)
+          model%geometry%cell_area = &
+               model%geometry%cell_area / model%projection%stere%area_factor**2
+       endwhere
+
     endif
 
     ! write projection info to log
@@ -292,9 +305,6 @@ subroutine glide_initialise(model)
 !    print*, 'Created Glide variables'
 !    print*, 'max, min bheatflx (W/m2)=', maxval(model%temper%bheatflx), minval(model%temper%bheatflx)
 
-    ! Compute the cell areas of the grid
-    model%geometry%cell_area = model%numerics%dew*model%numerics%dns
-
     ! If a 2D bheatflx field is present in the input file, it will have been written 
     !  to model%temper%bheatflx.  For the case model%options%gthf = 0, we want to use
     !  a uniform heat flux instead.
@@ -453,7 +463,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice)
        l_evolve_ice = .true.
     end if
 
-    if (model%options%is_restart == RESTART_TRUE) then
+    if (model%options%is_restart == STANDARD_RESTART .or. &
+        model%options%is_restart == HYBRID_RESTART) then
        ! On a restart, just assign the basal velocity from uvel/vvel (which are restart variables)
        ! to ubas/vbas which are used by the temperature solver to calculate basal heating.
        ! During time stepping ubas/vbas are calculated by slipvelo during thickness evolution or below on a cold start.