Merge branch 'NOAA-EMC:emc/develop' into feature/endofrun_restart

.github/workflows/test-cice.yml

@@ -46,7 +46,6 @@ jobs:
         run: |
           sudo xcode-select -r
           sudo xcode-select -s /Library/Developer/CommandLineTools
-          sudo ln -s /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/* /usr/local/include/
           echo "xcrun --show-sdk-path: $(xcrun --show-sdk-path)"
           echo "xcode-select -p: $(xcode-select -p)"
       - name: system info
@@ -55,7 +54,7 @@ jobs:
           type wget
           type curl
           type csh
-          echo "readlink \$(which csh): $(python -c 'import os, sys; print os.path.realpath(sys.argv[1])' $(which csh))"
+          echo "readlink \$(which csh): $(python -c 'import os, sys; print(os.path.realpath(sys.argv[1]))' $(which csh))"
           echo "csh --version: $(csh --version)"
           echo "uname -a: $(uname -a)"
           echo "sw_vers: $(sw_vers)"

cicecore/cicedyn/analysis/ice_history.F90

@@ -362,6 +362,7 @@ subroutine init_hist (dt)
          f_sidmasslat = 'mxxxx'
          f_sndmasssnf = 'mxxxx'
          f_sndmassmelt = 'mxxxx'
+         f_sndmassdyn = 'mxxxx'
          f_siflswdtop = 'mxxxx'
          f_siflswutop = 'mxxxx'
          f_siflswdbot = 'mxxxx'
@@ -402,6 +403,11 @@ subroutine init_hist (dt)
          f_siu = f_CMIP
          f_siv = f_CMIP
          f_sispeed = f_CMIP
+         f_sndmasssubl = f_CMIP
+         f_sndmasssnf = f_CMIP
+         f_sndmassmelt = f_CMIP
+         f_sndmassdyn = f_CMIP
+         f_sidmasssi = f_CMIP
       endif
 
       if (grid_ice == 'CD' .or. grid_ice == 'C') then
@@ -447,6 +453,14 @@ subroutine init_hist (dt)
       if (f_Tsnz (1:1) /= 'x')                          f_VGRDs = .true.
       if (tr_fsd)                                       f_NFSD  = .true.
 
+      call broadcast_scalar (f_tlon, master_task)
+      call broadcast_scalar (f_tlat, master_task)
+      call broadcast_scalar (f_ulon, master_task)
+      call broadcast_scalar (f_ulat, master_task)
+      call broadcast_scalar (f_nlon, master_task)
+      call broadcast_scalar (f_nlat, master_task)
+      call broadcast_scalar (f_elon, master_task)
+      call broadcast_scalar (f_elat, master_task)
       call broadcast_scalar (f_tmask, master_task)
       call broadcast_scalar (f_umask, master_task)
       call broadcast_scalar (f_nmask, master_task)
@@ -646,6 +660,7 @@ subroutine init_hist (dt)
       call broadcast_scalar (f_sidmasslat, master_task)
       call broadcast_scalar (f_sndmasssnf, master_task)
       call broadcast_scalar (f_sndmassmelt, master_task)
+      call broadcast_scalar (f_sndmassdyn, master_task)
       call broadcast_scalar (f_siflswdtop, master_task)
       call broadcast_scalar (f_siflswutop, master_task)
       call broadcast_scalar (f_siflswdbot, master_task)
@@ -1640,7 +1655,7 @@ subroutine init_hist (dt)
              "none", c1, c0,         &
              ns1, f_sidmassevapsubl)
 
-         call define_hist_field(n_sndmasssubl,"sndmassubl","kg m-2 s-1",tstr2D, tcstr,  &
+         call define_hist_field(n_sndmasssubl,"sndmasssubl","kg m-2 s-1",tstr2D, tcstr,  &
              "snow mass change from evaporation and sublimation", &
              "none", c1, c0,         &
              ns1, f_sndmasssubl)
@@ -1670,6 +1685,11 @@ subroutine init_hist (dt)
              "none", c1, c0,         &
              ns1, f_sndmassmelt)
 
+         call define_hist_field(n_sndmassdyn,"sndmassdyn","kg m-2 s-1",tstr2D, tcstr,  &
+             "snow mass change from dynamics ridging",                      &
+             "none", c1, c0,         &
+             ns1, f_sndmassdyn)
+
          call define_hist_field(n_siflswdtop,"siflswdtop","W/m2",tstr2D, tcstr, &
              "down shortwave flux over sea ice", &
              "positive downward", c1, c0,                            &
@@ -1973,6 +1993,21 @@ subroutine init_hist (dt)
       ! floe size distribution
        call init_hist_fsd_4Df
 
+      !-----------------------------------------------------------------
+      ! fill icoord array with namelist values
+      !-----------------------------------------------------------------
+
+       icoord=.true.
+
+       icoord(n_tlon   ) = f_tlon
+       icoord(n_tlat   ) = f_tlat
+       icoord(n_ulon   ) = f_ulon
+       icoord(n_ulat   ) = f_ulat
+       icoord(n_nlon   ) = f_nlon
+       icoord(n_nlat   ) = f_nlat
+       icoord(n_elon   ) = f_elon
+       icoord(n_elat   ) = f_elat
+
       !-----------------------------------------------------------------
       ! fill igrd array with namelist values
       !-----------------------------------------------------------------
@@ -2137,7 +2172,7 @@ subroutine accum_hist (dt)
           taubxN, taubyN, strocnxN, strocnyN, &
           strairxE, strairyE, strtltxE, strtltyE, strintxE, strintyE, &
           taubxE, taubyE, strocnxE, strocnyE, &
-          fmU, fmN, fmE, daidtt, dvidtt, daidtd, dvidtd, fsurf, &
+          fmU, fmN, fmE, daidtt, dvidtt, daidtd, dvidtd, dvsdtd, fsurf, &
           fcondtop, fcondbot, fsurfn, fcondtopn, flatn, fsensn, albcnt, snwcnt, &
           stressp_1, stressm_1, stress12_1, &
           stresspT, stressmT, stress12T, &
@@ -3045,7 +3080,7 @@ subroutine accum_hist (dt)
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = evaps(i,j,iblk)*rhos
+                 worka(i,j) = evaps(i,j,iblk)
               endif
            enddo
            enddo
@@ -3057,7 +3092,7 @@ subroutine accum_hist (dt)
            do j = jlo, jhi
            do i = ilo, ihi
               if (aice(i,j,iblk) > puny) then
-                 worka(i,j) = aice(i,j,iblk)*fsnow(i,j,iblk)*rhos
+                 worka(i,j) = aice(i,j,iblk)*fsnow(i,j,iblk)
               endif
            enddo
            enddo
@@ -3076,6 +3111,18 @@ subroutine accum_hist (dt)
            call accum_hist_field(n_sndmassmelt, iblk, worka(:,:), a2D)
          endif
 
+         if (f_sndmassdyn(1:1) /= 'x') then
+           worka(:,:) = c0
+           do j = jlo, jhi
+           do i = ilo, ihi
+              if (aice(i,j,iblk) > puny) then
+                 worka(i,j) = dvsdtd(i,j,iblk)*rhos
+              endif
+           enddo
+           enddo
+           call accum_hist_field(n_sndmassdyn, iblk, worka(:,:), a2D)
+         endif
+
          if (f_siflswdtop(1:1) /= 'x') then
            worka(:,:) = c0
            do j = jlo, jhi

cicecore/cicedyn/analysis/ice_history_shared.F90

@@ -57,7 +57,7 @@ module ice_history_shared
          history_rearranger      ! history file rearranger, box or subset for pio
 
       character (len=char_len), public :: &
-         hist_suffix(max_nstrm)  ! appended to 'h' in filename when not 'x'
+         hist_suffix(max_nstrm)  ! appended to history_file in filename
 
       integer (kind=int_kind), public :: &
          history_iotasks     , & ! iotasks, root, stride defines io pes for pio
@@ -131,6 +131,7 @@ module ice_history_shared
          avail_hist_fields(max_avail_hist_fields)
 
       integer (kind=int_kind), parameter, public :: &
+         ncoord   = 8           , & ! number of coordinate variables: TLON, TLAT, ULON, ULAT, NLON, NLAT, ELON, ELAT
          nvar_grd = 21          , & ! number of grid fields that can be written
                                     !   excluding grid vertices
          nvar_grdz = 6              ! number of category/vertical grid fields written
@@ -165,6 +166,7 @@ module ice_history_shared
          avgct(max_nstrm)   ! average sample counter
 
       logical (kind=log_kind), public :: &
+         icoord(ncoord) , &    ! true if coord field is written to output file
          igrd (nvar_grd), &    ! true if grid field is written to output file
          igrdz(nvar_grdz)      ! true if category/vertical grid field is written
 
@@ -194,6 +196,10 @@ module ice_history_shared
       !---------------------------------------------------------------
 
       logical (kind=log_kind), public :: &
+           f_tlon      = .true., f_tlat       = .true., &
+           f_ulon      = .true., f_ulat       = .true., &
+           f_nlon      = .true., f_nlat       = .true., &
+           f_elon      = .true., f_elat       = .true., &
            f_tmask     = .true., f_umask      = .true., &
            f_nmask     = .true., f_emask      = .true., &
            f_blkmask   = .true., &
@@ -308,6 +314,7 @@ module ice_history_shared
            f_sidmasslat = 'x', &
            f_sndmasssnf = 'x', &
            f_sndmassmelt = 'x', &
+           f_sndmassdyn = 'x', &
            f_siflswdtop = 'x', &
            f_siflswutop = 'x', &
            f_siflswdbot = 'x', &
@@ -362,6 +369,10 @@ module ice_history_shared
       !---------------------------------------------------------------
 
       namelist / icefields_nml /     &
+           f_tlon     , f_tlat     , &
+           f_ulon     , f_ulat     , &
+           f_nlon     , f_nlat     , &
+           f_elon     , f_elat     , &
            f_tmask    , f_umask    , &
            f_nmask    , f_emask    , &
            f_blkmask  , &
@@ -475,6 +486,7 @@ module ice_history_shared
            f_sidmasslat, &
            f_sndmasssnf, &
            f_sndmassmelt, &
+           f_sndmassdyn, &
            f_siflswdtop, &
            f_siflswutop, &
            f_siflswdbot, &
@@ -529,6 +541,15 @@ module ice_history_shared
       !---------------------------------------------------------------
 
       integer (kind=int_kind), parameter, public :: &
+           n_tlon       = 1,  &
+           n_tlat       = 2,  &
+           n_ulon       = 3,  &
+           n_ulat       = 4,  &
+           n_nlon       = 5,  &
+           n_nlat       = 6,  &
+           n_elon       = 7,  &
+           n_elat       = 8,  &
+
            n_tmask      = 1,  &
            n_umask      = 2,  &
            n_nmask      = 3,  &
@@ -665,6 +686,7 @@ module ice_history_shared
            n_sidmasslat,  &
            n_sndmasssnf,  &
            n_sndmassmelt,  &
+           n_sndmassdyn,  &
            n_siflswdtop,  &
            n_siflswutop,  &
            n_siflswdbot,  &
@@ -735,18 +757,22 @@ subroutine construct_filename(ncfile,suffix,ns)
                               dt
       use ice_restart_shared, only: lenstr
 
-      character (char_len_long), intent(inout) :: ncfile
-      character (len=2), intent(in) :: suffix
+      character (len=*), intent(inout) :: ncfile
+      character (len=*), intent(in) :: suffix
       integer (kind=int_kind), intent(in) :: ns
 
       integer (kind=int_kind) :: iyear, imonth, iday, isec
-      character (len=1) :: cstream
+      integer (kind=int_kind) :: n
+      character (len=char_len) :: cstream
+      character (len=char_len_long), save :: ncfile_last(max_nstrm) = 'UnDefineD'
       character(len=*), parameter :: subname = '(construct_filename)'
 
         iyear = myear
         imonth = mmonth
         iday = mday
         isec = int(msec - dt,int_kind)
+        cstream = ''
+        if (hist_suffix(ns) /= 'x') cstream = hist_suffix(ns)
 
         ! construct filename
         if (write_ic) then
@@ -771,9 +797,6 @@ subroutine construct_filename(ncfile,suffix,ns)
               endif
            endif
 
-           cstream = ''
-           if (hist_suffix(ns) /= 'x') cstream = hist_suffix(ns)
-
            if (hist_avg(ns)) then    ! write averaged data
               if (histfreq(ns) == '1' .and. histfreq_n(ns) == 1)  then ! timestep
                  write(ncfile,'(a,a,i4.4,a,i2.2,a,i2.2,a,i5.5,a,a)')  &
@@ -809,6 +832,25 @@ subroutine construct_filename(ncfile,suffix,ns)
 
         endif
 
+        ! Check whether the filename is already in use.
+        ! Same filename in multiple history streams leads to files being overwritten (not good).
+        ! The current filename convention means we just have to check latest filename, 
+        ! not all filenames ever generated because of use of current model date/time in filename.
+
+        ! write(nu_diag,'(2a,i2,1x,a)') subname, 'debug ncfile= ',ns,trim(ncfile)
+        do n = 1,max_nstrm
+           ! write(nu_diag,'(2a,i2,1x,a)') subname, 'debug nfile_last= ',n,trim(ncfile_last(n))
+           if (ncfile == ncfile_last(n)) then
+              write(nu_diag,*) subname,' history stream = ',ns
+              write(nu_diag,*) subname,' history filename = ',trim(ncfile)
+              write(nu_diag,*) subname,' filename in use for stream ',n
+              write(nu_diag,*) subname,' filename for stream ',trim(ncfile_last(n))
+              write(nu_diag,*) subname,' Use namelist hist_suffix so history filenames are unique'
+              call abort_ice(subname//' ERROR: history filename already used for another history stream '//trim(ncfile))
+           endif
+        enddo
+        ncfile_last(ns) = ncfile
+
       end subroutine construct_filename
 
 !=======================================================================
@@ -869,7 +911,7 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
       if(present(mask_ice_free_points)) l_mask_ice_free_points = mask_ice_free_points
 
       if (histfreq(ns) == 'x') then
-         call abort_ice(subname//'ERROR: define_hist_fields has histfreq x')
+         call abort_ice(subname//' ERROR: define_hist_fields has histfreq x')
       endif
 
       if (ns == 1) id(:) = 0
@@ -879,7 +921,7 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
          if (vhistfreq(ns1:ns1) == histfreq(ns)) then
 
             if (ns1 > 1 .and. index(vhistfreq(1:ns1-1),'x') /= 0) then
-               call abort_ice(subname//'ERROR: history frequency variable f_' // vname // ' can''t contain ''x'' along with active frequencies')
+               call abort_ice(subname//' ERROR: history frequency variable f_' // vname // ' can''t contain ''x'' along with active frequencies')
             endif
 
             num_avail_hist_fields_tot = num_avail_hist_fields_tot + 1
@@ -909,7 +951,7 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
                   write(nu_diag,*) subname,' num_avail_hist_fields_tot = ',num_avail_hist_fields_tot
                   write(nu_diag,*) subname,' max_avail_hist_fields     = ',max_avail_hist_fields
                endif
-               call abort_ice(subname//'ERROR: Need in computation of max_avail_hist_fields')
+               call abort_ice(subname//' ERROR: Need in computation of max_avail_hist_fields')
             endif
 
             if (num_avail_hist_fields_tot /= &
@@ -925,7 +967,7 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
                if (my_task == master_task) then
                   write(nu_diag,*) subname,' num_avail_hist_fields_tot = ',num_avail_hist_fields_tot
                endif
-               call abort_ice(subname//'ERROR: in num_avail_hist_fields')
+               call abort_ice(subname//' ERROR: in num_avail_hist_fields')
             endif
 
             id(ns) = num_avail_hist_fields_tot

cicecore/cicedyn/general/ice_flux.F90

@@ -107,6 +107,7 @@ module ice_flux
          strintyE, & ! divergence of internal ice stress, y at E points (N/m^2)
          daidtd  , & ! ice area tendency due to transport   (1/s)
          dvidtd  , & ! ice volume tendency due to transport (m/s)
+         dvsdtd  , & ! snow volume tendency due to transport (m/s)
          dagedtd , & ! ice age tendency due to transport (s/s)
          dardg1dt, & ! rate of area loss by ridging ice (1/s)
          dardg2dt, & ! rate of area gain by new ridges (1/s)
@@ -319,6 +320,7 @@ module ice_flux
          dsnow,  & ! change in snow thickness (m/step-->cm/day)
          daidtt, & ! ice area tendency thermo.   (s^-1)
          dvidtt, & ! ice volume tendency thermo. (m/s)
+         dvsdtt, & ! snow volume tendency thermo. (m/s)
          dagedtt,& ! ice age tendency thermo.    (s/s)
          mlt_onset, &! day of year that sfc melting begins
          frz_onset, &! day of year that freezing begins (congel or frazil)
@@ -419,6 +421,7 @@ subroutine alloc_flux
          strintyU   (nx_block,ny_block,max_blocks), & ! divergence of internal ice stress, y (N/m^2)
          daidtd     (nx_block,ny_block,max_blocks), & ! ice area tendency due to transport   (1/s)
          dvidtd     (nx_block,ny_block,max_blocks), & ! ice volume tendency due to transport (m/s)
+         dvsdtd     (nx_block,ny_block,max_blocks), & ! snow volume tendency due to transport (m/s)
          dagedtd    (nx_block,ny_block,max_blocks), & ! ice age tendency due to transport (s/s)
          dardg1dt   (nx_block,ny_block,max_blocks), & ! rate of area loss by ridging ice (1/s)
          dardg2dt   (nx_block,ny_block,max_blocks), & ! rate of area gain by new ridges (1/s)
@@ -530,6 +533,7 @@ subroutine alloc_flux
          dsnow      (nx_block,ny_block,max_blocks), & ! change in snow thickness (m/step-->cm/day)
          daidtt     (nx_block,ny_block,max_blocks), & ! ice area tendency thermo.   (s^-1)
          dvidtt     (nx_block,ny_block,max_blocks), & ! ice volume tendency thermo. (m/s)
+         dvsdtt     (nx_block,ny_block,max_blocks), & ! snow volume tendency thermo. (m/s)
          dagedtt    (nx_block,ny_block,max_blocks), & ! ice age tendency thermo.    (s/s)
          mlt_onset  (nx_block,ny_block,max_blocks), & ! day of year that sfc melting begins
          frz_onset  (nx_block,ny_block,max_blocks), & ! day of year that freezing begins (congel or frazil)
@@ -918,7 +922,7 @@ end subroutine init_flux_ocn
 
       subroutine init_history_therm
 
-      use ice_state, only: aice, vice, trcr
+      use ice_state, only: aice, vice, vsno, trcr
       use ice_arrays_column, only: &
           hfreebd, hdraft, hridge, distrdg, hkeel, dkeel, lfloe, dfloe, &
           Cdn_atm_skin, Cdn_atm_floe, Cdn_atm_pond, Cdn_atm_rdg, &
@@ -965,6 +969,7 @@ subroutine init_history_therm
       meltl  (:,:,:) = c0
       daidtt (:,:,:) = aice(:,:,:) ! temporary initial area
       dvidtt (:,:,:) = vice(:,:,:) ! temporary initial volume
+      dvsdtt (:,:,:) = vsno(:,:,:) ! temporary initial volume
       if (tr_iage) then
          dagedtt(:,:,:) = trcr(:,:,nt_iage,:) ! temporary initial age
       else
@@ -1022,7 +1027,7 @@ end subroutine init_history_therm
 
       subroutine init_history_dyn
 
-      use ice_state, only: aice, vice, trcr, strength, divu, shear, vort
+      use ice_state, only: aice, vice, vsno, trcr, strength, divu, shear, vort
       use ice_grid,  only: grid_ice
 
       logical (kind=log_kind) :: &
@@ -1061,6 +1066,7 @@ subroutine init_history_dyn
       opening (:,:,:) = c0
       daidtd  (:,:,:) = aice(:,:,:) ! temporary initial area
       dvidtd  (:,:,:) = vice(:,:,:) ! temporary initial volume
+      dvsdtd  (:,:,:) = vsno(:,:,:) ! temporary initial volume
       if (tr_iage) &
          dagedtd (:,:,:) = trcr(:,:,nt_iage,:) ! temporary initial age
       fmU     (:,:,:) = c0