Added pressure, modified norm of principal stresses and made small mo…

…difs to basal stress following Till's comments

cicecore/cicedynB/analysis/ice_history.F90

@@ -5,7 +5,7 @@
 !
 ! The following variables are currently hard-wired as snapshots 
 !   (instantaneous rather than time-averages):
-!   divu, shear, sig1, sig2, trsig, mlt_onset, frz_onset, hisnap, aisnap
+!   divu, shear, sig1, sig2, sigI, trsig, mlt_onset, frz_onset, hisnap, aisnap
 !
 ! Options for histfreq: '1','h','d','m','y','x', where x means that
 !   output stream will not be used (recommended for efficiency).  
@@ -300,6 +300,7 @@ subroutine init_hist (dt)
       call broadcast_scalar (f_shear, master_task)
       call broadcast_scalar (f_sig1, master_task)
       call broadcast_scalar (f_sig2, master_task)
+      call broadcast_scalar (f_sigI, master_task)
       call broadcast_scalar (f_dvidtt, master_task)
       call broadcast_scalar (f_dvidtd, master_task)
       call broadcast_scalar (f_daidtt, master_task)
@@ -776,6 +777,11 @@ subroutine init_hist (dt)
              "norm. principal stress 2",                       &
              "sig2 is instantaneous", c1, c0,                  &
              ns1, f_sig2)
+
+         call define_hist_field(n_sigI,"sigI","1",ustr2D, ucstr, &
+             "ice pressure",                       &
+             "sigI is instantaneous", c1, c0,                  &
+             ns1, f_sigI)
 
          call define_hist_field(n_dvidtt,"dvidtt","cm/day",tstr2D, tcstr, &
              "volume tendency thermo",                                  &
@@ -1199,11 +1205,11 @@ subroutine accum_hist (dt)
           fhocn, fhocn_ai, uatm, vatm, fbot, &
           fswthru_ai, strairx, strairy, strtltx, strtlty, strintx, strinty, &
           taubx, tauby, strocnx, strocny, fm, daidtt, dvidtt, daidtd, dvidtd, fsurf, &
-          fcondtop, fsurfn, fcondtopn, flatn, fsensn, albcnt, prs_sig, &
+          fcondtop, fsurfn, fcondtopn, flatn, fsensn, albcnt, &
           stressp_1, stressm_1, stress12_1, &
           stressp_2, stressm_2, stress12_2, &
           stressp_3, stressm_3, stress12_3, &
-          stressp_4, stressm_4, stress12_4, sig1, sig2, &
+          stressp_4, stressm_4, stress12_4, sig1, sig2, sigI, &
           mlt_onset, frz_onset, dagedtt, dagedtd, fswint_ai, keffn_top, &
           snowfrac, alvdr_ai, alvdf_ai, alidr_ai, alidf_ai
       use ice_arrays_column, only: snowfracn
@@ -1922,9 +1928,10 @@ subroutine accum_hist (dt)
                                   stressp_1 (:,:,iblk), &
                                   stressm_1 (:,:,iblk), &
                                   stress12_1(:,:,iblk), &
-                                  prs_sig   (:,:,iblk), &
+                                  strength  (:,:,iblk), &
                                   sig1      (:,:,iblk), &
-                                  sig2      (:,:,iblk))
+                                  sig2      (:,:,iblk), &
+                                  sigI      (:,:,iblk))
 
            do j = jlo, jhi
            do i = ilo, ihi
@@ -1933,6 +1940,7 @@ subroutine accum_hist (dt)
                  if (n_shear    (ns) /= 0) a2D(i,j,n_shear(ns),    iblk) = spval
                  if (n_sig1     (ns) /= 0) a2D(i,j,n_sig1(ns),     iblk) = spval
                  if (n_sig2     (ns) /= 0) a2D(i,j,n_sig2(ns),     iblk) = spval
+                 if (n_sigI     (ns) /= 0) a2D(i,j,n_sigI(ns),     iblk) = spval
                  if (n_mlt_onset(ns) /= 0) a2D(i,j,n_mlt_onset(ns),iblk) = spval
                  if (n_frz_onset(ns) /= 0) a2D(i,j,n_frz_onset(ns),iblk) = spval
                  if (n_hisnap   (ns) /= 0) a2D(i,j,n_hisnap(ns),   iblk) = spval
@@ -1961,6 +1969,8 @@ subroutine accum_hist (dt)
                        sig1 (i,j,iblk)*avail_hist_fields(n_sig1(ns))%cona
                  if (n_sig2     (ns) /= 0) a2D(i,j,n_sig2(ns),iblk)      = &
                        sig2 (i,j,iblk)*avail_hist_fields(n_sig2(ns))%cona
+                 if (n_sigI     (ns) /= 0) a2D(i,j,n_sigI(ns),iblk)      = &
+                       sigI (i,j,iblk)*avail_hist_fields(n_sigI(ns))%cona      
                  if (n_mlt_onset(ns) /= 0) a2D(i,j,n_mlt_onset(ns),iblk) = &
                        mlt_onset(i,j,iblk)
                  if (n_frz_onset(ns) /= 0) a2D(i,j,n_frz_onset(ns),iblk) = &
@@ -1972,7 +1982,7 @@ subroutine accum_hist (dt)
 
                  if (kdyn == 2) then  ! for EAP dynamics different time of output
                     if (n_trsig    (ns) /= 0) a2D(i,j,n_trsig(ns),iblk ) = &
-                                        prs_sig(i,j,iblk)
+                                        strength(i,j,iblk)
                  else
                     if (n_trsig    (ns) /= 0) a2D(i,j,n_trsig(ns),iblk ) = &
                                        p25*(stressp_1(i,j,iblk) &

cicecore/cicedynB/analysis/ice_history_shared.F90

@@ -5,7 +5,7 @@
 !
 ! The following variables are currently hard-wired as snapshots 
 !   (instantaneous rather than time-averages):
-!   divu, shear, sig1, sig2, trsig, mlt_onset, frz_onset, hisnap, aisnap
+!   divu, shear, sig1, sig2, sigI, trsig, mlt_onset, frz_onset, hisnap, aisnap
 !
 ! Options for histfreq: '1','h','d','m','y','x', where x means that
 !   output stream will not be used (recommended for efficiency).  
@@ -227,6 +227,7 @@ module ice_history_shared
            f_strength  = 'm', &
            f_divu      = 'm', f_shear      = 'm', &
            f_sig1      = 'm', f_sig2       = 'm', &
+           f_sigI      = 'm', &
            f_dvidtt    = 'm', f_dvidtd     = 'm', &
            f_daidtt    = 'm', f_daidtd     = 'm', &
            f_dagedtt   = 'm', f_dagedtd    = 'm', &
@@ -314,6 +315,7 @@ module ice_history_shared
            f_strength,  &
            f_divu,      f_shear    , &
            f_sig1,      f_sig2     , &
+           f_sigI,      &
            f_dvidtt,    f_dvidtd   , &
            f_daidtt,    f_daidtd   , &
            f_dagedtt,   f_dagedtd  , &
@@ -419,6 +421,7 @@ module ice_history_shared
            n_strength   , &
            n_divu       , n_shear      , &
            n_sig1       , n_sig2       , &
+           n_sigI       , &
            n_dvidtt     , n_dvidtd     , &
            n_daidtt     , n_daidtd     , &
            n_dagedtt    , n_dagedtd    , &

cicecore/cicedynB/dynamics/ice_dyn_eap.F90

@@ -93,7 +93,7 @@ subroutine eap (dt)
           cosw, sinw, denom1, uvel_init, vvel_init, arlx1i, &
           evp_prep1, evp_prep2, stepu, evp_finish, &
           basal_stress_coeff, basalstress
-      use ice_flux, only: rdg_conv, rdg_shear, prs_sig, strairxT, strairyT, &
+      use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
           strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, &
           strtltx, strtlty, strocnx, strocny, strintx, strinty, &
           strocnxT, strocnyT, strax, stray, &
@@ -176,7 +176,6 @@ subroutine eap (dt)
             rdg_shear(i,j,iblk) = c0 
             divu (i,j,iblk) = c0 
             shear(i,j,iblk) = c0 
-            prs_sig(i,j,iblk) = c0 
             e11(i,j,iblk) = c0
             e12(i,j,iblk) = c0
             e22(i,j,iblk) = c0
@@ -394,7 +393,6 @@ subroutine eap (dt)
                               yieldstress11 (:,:,iblk),                   &
                               yieldstress12 (:,:,iblk),                   &
                               yieldstress22 (:,:,iblk),                   &
-                              prs_sig   (:,:,iblk),                       &
                               rdg_conv  (:,:,iblk), rdg_shear (:,:,iblk), &
                               strtmp    (:,:,:))
 !      call ice_timer_stop(timer_tmp1) ! dynamics
@@ -486,8 +484,8 @@ subroutine eap (dt)
       if ( basalstress ) then
          !$OMP PARALLEL DO PRIVATE(iblk)
          do iblk = 1, nblocks
-            taubx(:,:,iblk) = Cbu(:,:,iblk)*uvel(:,:,iblk)
-            tauby(:,:,iblk) = Cbu(:,:,iblk)*vvel(:,:,iblk)
+            taubx(:,:,iblk) = -Cbu(:,:,iblk)*uvel(:,:,iblk)
+            tauby(:,:,iblk) = -Cbu(:,:,iblk)*vvel(:,:,iblk)
          enddo
          !$OMP END PARALLEL DO
       endif
@@ -1145,7 +1143,6 @@ subroutine stress_eap  (nx_block,   ny_block,       &
                               yieldstress11,              &
                               yieldstress12,              &
                               yieldstress22,              &
-                              prs_sig,                    &
                               rdg_conv,   rdg_shear,      &
                               strtmp)
 
@@ -1196,7 +1193,6 @@ subroutine stress_eap  (nx_block,   ny_block,       &
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), & 
          intent(inout) :: &
-         prs_sig  , & ! replacement pressure, for stress calc
          shear    , & ! strain rate II component (1/s)
          divu     , & ! strain rate I component, velocity divergence (1/s)
          e11      , & ! components of strain rate tensor (1/s)
@@ -1353,8 +1349,6 @@ subroutine stress_eap  (nx_block,   ny_block,       &
          e22(i,j) = p5*p25*(divune + divunw + divuse + divusw - &
                     tensionne - tensionnw - tensionse - tensionsw) * tarear(i,j)
 
-         prs_sig(i,j) = strength(i,j) 
-
       !-----------------------------------------------------------------
       ! elastic relaxation, see Eq. A12-A14
       !-----------------------------------------------------------------

cicecore/cicedynB/dynamics/ice_dyn_evp.F90

@@ -80,7 +80,7 @@ subroutine evp (dt)
       use ice_blocks, only: block, get_block, nx_block, ny_block
       use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn
       use ice_domain_size, only: max_blocks, ncat
-      use ice_flux, only: rdg_conv, rdg_shear, prs_sig, strairxT, strairyT, &
+      use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
           strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, &
           strtltx, strtlty, strocnx, strocny, strintx, strinty, &
           strocnxT, strocnyT, strax, stray, &
@@ -177,7 +177,6 @@ subroutine evp (dt)
             rdg_shear(i,j,iblk) = c0 
             divu (i,j,iblk) = c0 
             shear(i,j,iblk) = c0 
-            prs_sig(i,j,iblk) = c0 
          enddo
          enddo
 
@@ -352,7 +351,6 @@ subroutine evp (dt)
                             stress12_1(:,:,iblk), stress12_2(:,:,iblk), & 
                             stress12_3(:,:,iblk), stress12_4(:,:,iblk), & 
                             shear     (:,:,iblk), divu      (:,:,iblk), & 
-                            prs_sig   (:,:,iblk),                       & 
                             rdg_conv  (:,:,iblk), rdg_shear (:,:,iblk), & 
                             strtmp    (:,:,:) )
 !            endif               ! yield_curve
@@ -420,8 +418,8 @@ subroutine evp (dt)
       if ( basalstress ) then
          !$OMP PARALLEL DO PRIVATE(iblk)
          do iblk = 1, nblocks
-            taubx(:,:,iblk) = Cbu(:,:,iblk)*uvel(:,:,iblk)
-            tauby(:,:,iblk) = Cbu(:,:,iblk)*vvel(:,:,iblk)
+            taubx(:,:,iblk) = -Cbu(:,:,iblk)*uvel(:,:,iblk)
+            tauby(:,:,iblk) = -Cbu(:,:,iblk)*vvel(:,:,iblk)
          enddo
          !$OMP END PARALLEL DO
       endif
@@ -551,7 +549,6 @@ subroutine stress (nx_block,   ny_block,   &
                          stress12_1, stress12_2, & 
                          stress12_3, stress12_4, & 
                          shear,      divu,       & 
-                         prs_sig,                & 
                          rdg_conv,   rdg_shear,  & 
                          str )
 
@@ -588,7 +585,6 @@ subroutine stress (nx_block,   ny_block,   &
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), & 
          intent(inout) :: &
-         prs_sig  , & ! replacement pressure, for stress calc
          shear    , & ! strain rate II component (1/s)
          divu     , & ! strain rate I component, velocity divergence (1/s)
          rdg_conv , & ! convergence term for ridging (1/s)
@@ -698,7 +694,6 @@ subroutine stress (nx_block,   ny_block,   &
          c0nw = strength(i,j)/max(Deltanw,tinyarea(i,j))
          c0sw = strength(i,j)/max(Deltasw,tinyarea(i,j))
          c0se = strength(i,j)/max(Deltase,tinyarea(i,j))
-         prs_sig(i,j) = c0ne*Deltane ! northeast
 
          c1ne = c0ne*arlx1i
          c1nw = c0nw*arlx1i

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

@@ -851,6 +851,8 @@ subroutine basal_stress_coeff (nx_block, ny_block, icellu, &
                                      uold,     vold,             &
                                      Cbu)
 
+      use ice_constants, only: c0, c1                                     
+
       integer (kind=int_kind), intent(in) :: &
          nx_block, ny_block, &  ! block dimensions
          icellu                 ! no. of cells where icetmask = 1
@@ -878,7 +880,7 @@ subroutine basal_stress_coeff (nx_block, ny_block, icellu, &
          hu,  & ! volume per unit area of ice at u location (mean thickness)
          hwu, & ! water depth at u location
          hcu, & ! critical thickness at u location
-         k1 = 8.0_dbl_kind , &  ! first free parameter for landfast parametrization 
+         k1 = 80.0_dbl_kind , &  ! first free parameter for landfast parametrization 
          k2 = 15.0_dbl_kind, &  ! second free parameter (Nm^-3) for landfast parametrization 
          u0 = 5e-5_dbl_kind, &  ! residual velocity (m/s)
          CC = 20.0_dbl_kind     ! CC=Cb factor in Lemieux et al 2015
@@ -900,11 +902,9 @@ subroutine basal_stress_coeff (nx_block, ny_block, icellu, &
          hcu = au * hwu / k1
 
          ! 2- calculate stress factor
-         if (au > p01 .and. hu > hcu ) then
-   !       endif
-           Cbu(i,j) = ( k2 / (sqrt(uold(i,j)**2 + vold(i,j)**2) + u0) ) &
-                      * (hu - hcu) * exp(-CC * (1 - au))
-         endif
+
+         Cbu(i,j) = ( k2 / (sqrt(uold(i,j)**2 + vold(i,j)**2) + u0) ) &
+                    * max(c0,(hu - hcu)) * exp(-CC * (c1 - au))
 
       enddo                     ! ij
 
@@ -919,8 +919,9 @@ end subroutine basal_stress_coeff
 
       subroutine principal_stress(nx_block,   ny_block,  &
                                   stressp_1,  stressm_1, &
-                                  stress12_1, prs_sig,   &
-                                  sig1,       sig2)
+                                  stress12_1, strength,  &
+                                  sig1,       sig2,      &
+                                  sigI)
 
       use ice_constants, only: spval_dbl, p5, c4
 
@@ -931,11 +932,12 @@ subroutine principal_stress(nx_block,   ny_block,  &
          stressp_1 , & ! sigma11 + sigma22
          stressm_1 , & ! sigma11 - sigma22
          stress12_1, & ! sigma12
-         prs_sig       ! replacement pressure, for stress calc
+         strength      ! for normalization of sig1 and sig2
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: &
-         sig1    , & ! principal stress component
-         sig2        ! principal stress component
+         sig1    , & ! normalized principal stress component
+         sig2    , & ! normalized principal stress component
+         sigI        ! internal ice pressure (N/m)
 
       ! local variables
 
@@ -950,16 +952,20 @@ subroutine principal_stress(nx_block,   ny_block,  &
 
       do j = 1, ny_block
       do i = 1, nx_block
-         if (prs_sig(i,j) > puny) then
+         if (strength(i,j) > puny) then
+            ! not normalized yet
             sig1(i,j) = (p5*(stressp_1(i,j) &
-                      + sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2))) &
-                      / prs_sig(i,j)
+                      + sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2)))
             sig2(i,j) = (p5*(stressp_1(i,j) &
-                      - sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2))) &
-                      / prs_sig(i,j)
+                      - sqrt(stressm_1(i,j)**2+c4*stress12_1(i,j)**2)))
+            sigI(i,j) = -p5*( sig1(i,j) + sig2(i,j) ) 
+            ! normalization of sig1 and sig2
+            sig1(i,j) = sig1(i,j) / strength(i,j)
+            sig2(i,j) = sig2(i,j) / strength(i,j)
          else
             sig1(i,j) = spval_dbl
             sig2(i,j) = spval_dbl
+            sigI(i,j) = spval_dbl
          endif
       enddo
       enddo

cicecore/cicedynB/general/ice_flux.F90

@@ -58,8 +58,9 @@ module ice_flux
        ! diagnostic
 
       real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), public :: &
-         sig1    , & ! principal stress component
-         sig2    , & ! principal stress component
+         sig1    , & ! normalized principal stress component
+         sig2    , & ! normalized principal stress component
+         sigI    , & ! internal ice pressure (N/m)
          taubx   , & ! basal stress (x) (N/m^2)
          tauby   , & ! basal stress (y) (N/m^2)
          strairx , & ! stress on ice by air, x-direction
@@ -108,7 +109,6 @@ module ice_flux
        ! internal
 
       real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), public :: &
-         prs_sig  , & ! replacement pressure, for stress calc
          fm       , & ! Coriolis param. * mass in U-cell (kg/s)
          Cbu          ! coefficient for basal stress (landfast ice)
 
@@ -721,7 +721,6 @@ subroutine init_history_dyn
       if (tr_iage) &
          dagedtd (:,:,:) = trcr(:,:,nt_iage,:) ! temporary initial age
       fm      (:,:,:) = c0
-      prs_sig (:,:,:) = c0
       ardgn   (:,:,:,:) = c0
       vrdgn   (:,:,:,:) = c0
       krdgn   (:,:,:,:) = c1

cicecore/cicedynB/infrastructure/io/io_binary/ice_history_write.F90

@@ -147,7 +147,7 @@ subroutine ice_write_hist(ns)
             if (histfreq(ns) == '1' .or. .not. hist_avg         &
                 .or. n==n_divu(ns)      .or. n==n_shear(ns)     &  ! snapshots
                 .or. n==n_sig1(ns)      .or. n==n_sig2(ns)      &
-                .or. n==n_trsig(ns)                             &
+                .or. n==n_sigI(ns)      .or. n==n_trsig(ns)     &
                 .or. n==n_mlt_onset(ns) .or. n==n_frz_onset(ns) &
                 .or. n==n_hisnap(ns)    .or. n==n_aisnap(ns)) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &

cicecore/cicedynB/infrastructure/io/io_netcdf/ice_history_write.F90

@@ -505,7 +505,8 @@ subroutine ice_write_hist (ns)
       !-----------------------------------------------------------------
             if (hist_avg) then
               if (TRIM(avail_hist_fields(n)%vname)/='sig1' &
-              .or.TRIM(avail_hist_fields(n)%vname)/='sig2') then
+              .or.TRIM(avail_hist_fields(n)%vname)/='sig2' & 
+              .or.TRIM(avail_hist_fields(n)%vname)/='sigI') then
                 status = nf90_put_att(ncid,varid,'cell_methods','time: mean')
                 if (status /= nf90_noerr) call abort_ice( &
                  'Error defining cell methods for '//avail_hist_fields(n)%vname)
@@ -515,7 +516,7 @@ subroutine ice_write_hist (ns)
             if (histfreq(ns) == '1' .or. .not. hist_avg         &
                 .or. n==n_divu(ns)      .or. n==n_shear(ns)     &  ! snapshots
                 .or. n==n_sig1(ns)      .or. n==n_sig2(ns)      &
-                .or. n==n_trsig(ns)                             &
+                .or. n==n_sigI(ns)      .or. n==n_trsig(ns)     &
                 .or. n==n_mlt_onset(ns) .or. n==n_frz_onset(ns) &
                 .or. n==n_hisnap(ns)    .or. n==n_aisnap(ns)) then
                status = nf90_put_att(ncid,varid,'time_rep','instantaneous')

cicecore/cicedynB/infrastructure/io/io_pio/ice_history_write.F90

@@ -425,15 +425,16 @@ subroutine ice_write_hist (ns)
             ! Add cell_methods attribute to variables if averaged
             if (hist_avg .and. histfreq(ns) /= '1') then
               if (TRIM(avail_hist_fields(n)%vname)/='sig1' &
-              .or.TRIM(avail_hist_fields(n)%vname)/='sig2') then
+              .or.TRIM(avail_hist_fields(n)%vname)/='sig2' &
+              .or.TRIM(avail_hist_fields(n)%vname)/='sigI') then
                 status = pio_put_att(File,varid,'cell_methods','time: mean')
               endif
             endif
 
             if (histfreq(ns) == '1' .or. .not. hist_avg         &
                 .or. n==n_divu(ns)      .or. n==n_shear(ns)     &  ! snapshots
                 .or. n==n_sig1(ns)      .or. n==n_sig2(ns)      &
-                .or. n==n_trsig(ns)                             &
+                .or. n==n_sigI(ns)      .or. n==n_trsig(ns)     &
                 .or. n==n_mlt_onset(ns) .or. n==n_frz_onset(ns) &
                 .or. n==n_hisnap(ns)    .or. n==n_aisnap(ns)) then
                status = pio_put_att(File,varid,'time_rep','instantaneous')