Porting miscellaneous changes from the various column package branche…

…s back to

the trunk.  BFB in the standard configuration, but the initializations and 
conditional changes for coszen could change the answers in other configurations. 
Also the flux calculation change in ice_therm_itd.F90 could change the answers 
in coupled simulations.  PIO changes are not tested.

Regarding logic behind the ice_therm_itd.F90 modification, there are two namelist 
flags in play here:

update_ocn_f = true if the ocean model being coupled does not assume fresh water 
and salt content of frazil ice (e.g. HadGEM, ACME); otherwise update_ocn_f = false.

ktherm = 2 includes the fresh water and salt content of frazil ice in its 
calculations, consistent with update_ocn_f = true, and therefore if 
update_ocn_f = false then these need to be removed from the fluxes. 
ktherm = 0 and 1 do not include fresh water and salt content in their 
calculations, and therefore if update_ocn_f = true then they need to be added.


M       source/ice_shortwave.F90
 - initialize scale_factor
 - change gamma to agamm (to avoid conflict with Fortran intrinsic)

 M      source/ice_firstyear.F90
 - removed executable property

M       source/ice_flux.F90
 - initialize coszen

M       source/ice_therm_itd.F90
 - alter fresh and fsalt fluxes calculation depending on update_ocn_f and ktherm

M       source/ice_history_shared.F90
 - bug fix for num_avail_hist_fields_3Dz

M       source/ice_spacecurve.F90
 - initialize res

M       source/ice_history.F90
 - remove unneeded variable declarations

M       input_templates/run_ice.Linux.LANL.conejo
M       comp_ice
 - change path to run directory at LANL

M       drivers/cice/CICE_RunMod.F90
M       drivers/cesm/CICE_RunMod.F90
M       drivers/hadgem3/CICE_RunMod.F90
 - include = in coszen conditional 

M       io_pio/ice_pio.F90
 - change dof*d from pio_offset_kind to int_kind

M       io_pio/ice_history_write.F90
 - all processors call date_and_time
 - change PIO_OFFSET_KIND to PIO_OFFSET

cice/comp_ice

@@ -7,7 +7,7 @@ setenv SITE LANL.conejo   # also works for mustang
 #setenv SITE LANL.bigsuze
 
 ### SYSTEM_USERDIR is predefined on ORNL machines
-setenv SYSTEM_USERDIR /scratch/eclare/CICE.v5.1  # conejo
+setenv SYSTEM_USERDIR /net/scratch1/eclare/CICE.v5.1  # conejo
 #setenv SYSTEM_USERDIR /Users/akt/Work/MPAS-CICE/Column_Package/cice_package/run_dirs
 
 ### Grid resolution

cice/drivers/cesm/CICE_RunMod.F90

@@ -366,7 +366,7 @@ subroutine coupling_prep (iblk)
             alidr(i,j,iblk) = alidr(i,j,iblk) &
                + alidrn(i,j,n,iblk)*aicen(i,j,n,iblk)
 
-            if (coszen(i,j,iblk) > puny) then ! sun above horizon
+            if (coszen(i,j,iblk) >= puny) then ! sun above horizon
             albice(i,j,iblk) = albice(i,j,iblk) &
                + albicen(i,j,n,iblk)*aicen(i,j,n,iblk)
             albsno(i,j,iblk) = albsno(i,j,iblk) &

cice/drivers/cice/CICE_RunMod.F90

@@ -357,7 +357,7 @@ subroutine coupling_prep (iblk)
             alidr(i,j,iblk) = alidr(i,j,iblk) &
                + alidrn(i,j,n,iblk)*aicen(i,j,n,iblk)
 
-            if (coszen(i,j,iblk) > puny) then ! sun above horizon
+            if (coszen(i,j,iblk) >= puny) then ! sun above horizon
             albice(i,j,iblk) = albice(i,j,iblk) &
                + albicen(i,j,n,iblk)*aicen(i,j,n,iblk)
             albsno(i,j,iblk) = albsno(i,j,iblk) &

cice/drivers/hadgem3/CICE_RunMod.F90

@@ -367,7 +367,7 @@ subroutine coupling_prep (iblk)
             alidr(i,j,iblk) = alidr(i,j,iblk) &
                + alidrn(i,j,n,iblk)*aicen(i,j,n,iblk)
 
-            if (coszen(i,j,iblk) > puny) then ! sun above horizon
+            if (coszen(i,j,iblk) >= puny) then ! sun above horizon
             albice(i,j,iblk) = albice(i,j,iblk) &
                + albicen(i,j,n,iblk)*aicen(i,j,n,iblk)
             albsno(i,j,iblk) = albsno(i,j,iblk) &

cice/input_templates/run_ice.Linux.LANL.conejo

@@ -10,7 +10,7 @@ limit stacksize unlimited
 set  stamp   = `date -u '+%y%m%d%H%M%S'`
 setenv ICE_LOG_FILE ice.log.${stamp}
 
-cd /scratch/eclare/CICE.v5.1/rundir
+cd /net/scratch1/eclare/CICE.v5.1/rundir
 mpirun -np 4 ./cice  >&! $ICE_LOG_FILE
 
 exit 0

cice/io_pio/ice_history_write.F90

@@ -718,11 +718,7 @@ subroutine ice_write_hist (ns)
         status =  &
              pio_put_att(File,pio_global,'conventions',trim(title))
 
-        if (my_task == master_task) then
-           call date_and_time(date=current_date, time=current_time)
-        endif
-        call broadcast_scalar(current_date, master_task)
-        call broadcast_scalar(current_time, master_task)
+        call date_and_time(date=current_date, time=current_time)
         write(start_time,1000) current_date(1:4), current_date(5:6), &
                                current_date(7:8), current_time(1:2), &
                                current_time(3:4)
@@ -894,7 +890,7 @@ subroutine ice_write_hist (ns)
             if (status /= pio_noerr) call abort_ice( &
                'ice: Error getting varid for '//avail_hist_fields(n)%vname)
             workr2(:,:,:) = a2D(:,:,n,1:nblocks)
-            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET_KIND))
+            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET))
             call pio_write_darray(File, varid, iodesc2d,&
                                   workr2, status, fillval=spval_dbl)
          endif
@@ -915,7 +911,7 @@ subroutine ice_write_hist (ns)
                workr3(:,:,j,i) = a3Dc(:,:,i,nn,j)
             enddo
             enddo
-            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET_KIND))
+            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET))
             call pio_write_darray(File, varid, iodesc3dc,&
                                   workr3, status, fillval=spval_dbl)
          endif
@@ -935,7 +931,7 @@ subroutine ice_write_hist (ns)
                workr3(:,:,j,i) = a3Dz(:,:,i,nn,j)
             enddo
             enddo
-            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET_KIND))
+            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET))
             call pio_write_darray(File, varid, iodesc3di,&
                                   workr3, status, fillval=spval_dbl)
          endif
@@ -955,7 +951,7 @@ subroutine ice_write_hist (ns)
                workr3(:,:,j,i) = a3Db(:,:,i,nn,j)
             enddo
             enddo
-            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET_KIND))
+            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET))
             call pio_write_darray(File, varid, iodesc3db,&
                                   workr3, status, fillval=spval_dbl)
          endif
@@ -977,7 +973,7 @@ subroutine ice_write_hist (ns)
             enddo ! k
             enddo ! i
             enddo ! j
-            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET_KIND))
+            call pio_setframe(File, varid, int(1,kind=PIO_OFFSET))
             call pio_write_darray(File, varid, iodesc4di,&
                                   workr4, status, fillval=spval_dbl)
          endif

cice/io_pio/ice_pio.F90

@@ -136,7 +136,7 @@ subroutine ice_pio_initdecomp_2d(iodesc)
 
       type(block) :: this_block 
 
-      integer(kind=pio_offset_kind), pointer :: dof2d(:)
+      integer(kind=int_kind), pointer :: dof2d(:)
 
       allocate(dof2d(nx_block*ny_block*nblocks))
 
@@ -183,7 +183,7 @@ subroutine ice_pio_initdecomp_3d (ndim3, iodesc, remap)
 
       type(block) :: this_block 
       logical :: lremap
-      integer(kind=pio_offset_kind), pointer :: dof3d(:)
+      integer(kind=int_kind), pointer :: dof3d(:)
 
       allocate(dof3d(nx_block*ny_block*nblocks*ndim3))
       lremap=.false.
@@ -261,7 +261,7 @@ subroutine ice_pio_initdecomp_3d_inner(ndim3, inner_dim, iodesc)
 
       type(block) :: this_block 
 
-      integer(kind=pio_offset_kind), pointer :: dof3d(:)
+      integer(kind=int_kind), pointer :: dof3d(:)
 
       allocate(dof3d(nx_block*ny_block*nblocks*ndim3))
 
@@ -308,7 +308,7 @@ subroutine ice_pio_initdecomp_4d (ndim3, ndim4, iodesc)
 
       type(block) :: this_block 
 
-      integer(kind=pio_offset_kind), pointer :: dof4d(:)
+      integer(kind=int_kind), pointer :: dof4d(:)
 
       allocate(dof4d(nx_block*ny_block*nblocks*ndim3*ndim4))
 

cice/source/ice_flux.F90

@@ -485,6 +485,7 @@ subroutine init_coupler_flux
       ! derived or computed fields
       !-----------------------------------------------------------------
 
+      coszen  (:,:,:) = c0            ! cosine of the zenith angle
       fsw     (:,:,:) = c0            ! shortwave radiation (W/m^2)
       scale_factor(:,:,:) = c1        ! shortwave scaling factor 
       wind    (:,:,:) = sqrt(uatm(:,:,:)**2 &

cice/source/ice_history.F90

@@ -1151,7 +1151,7 @@ subroutine accum_hist (dt)
       use ice_blocks, only: block, get_block, nx_block, ny_block
       use ice_fileunits, only: nu_diag
       use ice_constants, only: c0, c1, p25, puny, secday, depressT, &
-          awtvdr, awtidr, awtvdf, awtidf, Lfresh, rhos, cp_ice, spval, hs_min
+          awtvdr, awtidr, awtvdf, awtidf, Lfresh, rhos, cp_ice, spval
       use ice_domain, only: blocks_ice, nblocks
       use ice_grid, only: tmask, lmask_n, lmask_s
       use ice_calendar, only: new_year, write_history, &
@@ -1176,7 +1176,6 @@ subroutine accum_hist (dt)
           mlt_onset, frz_onset, dagedtt, dagedtd, fswint_ai, keffn_top, &
           snowfrac, alvdr_ai, alvdf_ai, alidr_ai, alidf_ai
       use ice_atmo, only: formdrag
-      use ice_meltpond_cesm, only: hs0
       use ice_history_shared ! almost everything
       use ice_history_write, only: ice_write_hist
       use ice_history_bgc, only: accum_hist_bgc

cice/source/ice_history_shared.F90

@@ -595,7 +595,7 @@ subroutine define_hist_field(id, vname, vunit, vcoord, vcellmeas, &
                num_avail_hist_fields_2D  = num_avail_hist_fields_2D + 1
             elseif (vcoord(11:14) == 'NCAT' .and. vcoord(16:19) == 'time') then
                num_avail_hist_fields_3Dc = num_avail_hist_fields_3Dc + 1
-            elseif (vcoord(11:14) == 'VGRDi' .and. vcoord(16:19) == 'time') then
+            elseif (vcoord(11:15) == 'VGRDi' .and. vcoord(17:20) == 'time') then
                num_avail_hist_fields_3Dz = num_avail_hist_fields_3Dz + 1
             elseif (vcoord(11:15) == 'VGRDb' .and. vcoord(17:20) == 'time') then
                num_avail_hist_fields_3Db = num_avail_hist_fields_3Db + 1

cice/source/ice_shortwave.F90

@@ -142,14 +142,15 @@ subroutine init_shortwave
       use ice_domain, only: nblocks, blocks_ice
       use ice_flux, only: alvdf, alidf, alvdr, alidr, &
                           alvdr_ai, alidr_ai, alvdf_ai, alidf_ai, &
-                          swvdr, swvdf, swidr, swidf, &
+                          swvdr, swvdf, swidr, swidf, scale_factor, &
                           albice, albsno, albpnd, albcnt, coszen, fsnow, &
                           apeff_ai, snowfrac
       use ice_orbital, only: init_orbit
       use ice_state, only: aicen, vicen, vsnon, trcrn, nt_Tsfc
       use ice_blocks, only: block, get_block
       use ice_grid, only: tmask, tlat, tlon
       use ice_meltpond_lvl, only: dhsn, ffracn
+      use ice_restart_shared, only: restart, runtype
 
       integer (kind=int_kind) :: &
          icells          ! number of cells with aicen > puny
@@ -345,6 +346,21 @@ subroutine init_shortwave
          enddo
          enddo
 
+      !----------------------------------------------------------------
+      ! Save net shortwave for scaling factor in scale_factor
+      !----------------------------------------------------------------
+         if (runtype == 'initial' .and. .not. restart) then
+            do j = 1, ny_block
+            do i = 1, nx_block
+               scale_factor(i,j,iblk) = &
+	 	      swvdr(i,j,iblk)*(c1 - alvdr_ai(i,j,iblk)) &
+	 	    + swvdf(i,j,iblk)*(c1 - alvdf_ai(i,j,iblk)) &
+ 	            + swidr(i,j,iblk)*(c1 - alidr_ai(i,j,iblk)) &
+	 	    + swidf(i,j,iblk)*(c1 - alidf_ai(i,j,iblk))
+            enddo
+            enddo
+         endif
+
       enddo     ! nblocks
       !$OMP END PARALLEL DO
 
@@ -3265,8 +3281,8 @@ subroutine solution_dEdd                                 &
 
       real (kind=dbl_kind) :: &
          alpha    , & ! term in direct reflectivity and transmissivity
-         gamma    , & ! term in direct reflectivity and transmissivity
-         el       , & ! term in alpha,gamma,n,u
+         agamm    , & ! term in direct reflectivity and transmissivity
+         el       , & ! term in alpha,agamm,n,u
          taus     , & ! scaled extinction optical depth
          omgs     , & ! scaled single particle scattering albedo
          asys     , & ! scaled asymmetry parameter
@@ -3287,7 +3303,7 @@ subroutine solution_dEdd                                 &
 
       real (kind=dbl_kind) :: &
          alp      , & ! temporary for alpha
-         gam      , & ! temporary for gamma
+         gam      , & ! temporary for agamm
          ue       , & ! temporary for u
          extins   , & ! extinction
          amg      , & ! alp - gam
@@ -3322,7 +3338,7 @@ subroutine solution_dEdd                                 &
 
       ! Delta-Eddington solution expressions
       alpha(w,uu,gg,e) = p75*w*uu*((c1 + gg*(c1-w))/(c1 - e*e*uu*uu))
-      gamma(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu) &
+      agamm(w,uu,gg,e) = p5*w*((c1 + c3*gg*(c1-w)*uu*uu) &
                         / (c1-e*e*uu*uu))
       n(uu,et)         = ((uu+c1)*(uu+c1)/et ) - ((uu-c1)*(uu-c1)*et)
       u(w,gg,e)        = c1p5*(c1 - w*gg)/e
@@ -3430,7 +3446,7 @@ subroutine solution_dEdd                                 &
            ! evaluate rdir,tdir for direct beam
            trnlay(k) = max(exp_min, exp(-ts/mu0n))
            alp = alpha(ws,mu0n,gs,lm)
-           gam = gamma(ws,mu0n,gs,lm)
+           gam = agamm(ws,mu0n,gs,lm)
            apg = alp + gam
            amg = alp - gam
            rdir(k) = apg*rdif_a(k) +  amg*(tdif_a(k)*trnlay(k) - c1)
@@ -3451,7 +3467,7 @@ subroutine solution_dEdd                                 &
              swt = swt + mu*gwt
              trn = max(exp_min, exp(-ts/mu))
              alp = alpha(ws,mu,gs,lm)
-             gam = gamma(ws,mu,gs,lm)
+             gam = agamm(ws,mu,gs,lm)
              apg = alp + gam
              amg = alp - gam
              rdr = apg*R1 + amg*T1*trn - amg

cice/source/ice_spacecurve.F90

@@ -1139,6 +1139,7 @@ function FirstFactor(fac) result(res)
        integer (int_kind) :: i
 
        found = .false.
+       res = -1
        i=1
        do while (i<=fac%numfact .and. (.not. found))
           if(fac%used(i) == 0) then

cice/source/ice_therm_itd.F90

@@ -1365,6 +1365,7 @@ subroutine add_new_ice (nx_block,  ny_block,   &
          rnilyr       , & ! real(nilyr)
          dfresh       , & ! change in fresh
          dfsalt       , & ! change in fsalt
+         vi0tmp       , & ! frzmlt part of frazil
          Ti               ! frazil temperature
 
       real (kind=dbl_kind), dimension (icells) :: &
@@ -1524,11 +1525,22 @@ subroutine add_new_ice (nx_block,  ny_block,   &
       !-----------------------------------------------------------------
 
          if (update_ocn_f) then
-            dfresh = -rhoi*vi0new(ij)/dt 
-            dfsalt = ice_ref_salinity*p001*dfresh
-
-            fresh(i,j)      = fresh(i,j)      + dfresh
-            fsalt(i,j)      = fsalt(i,j)      + dfsalt
+            if (ktherm <= 1) then
+               dfresh = -rhoi*vi0new(ij)/dt 
+               dfsalt = ice_ref_salinity*p001*dfresh
+               fresh(i,j) = fresh(i,j) + dfresh
+               fsalt(i,j) = fsalt(i,j) + dfsalt
+            ! elseif (ktherm == 2) the fluxes are added elsewhere
+            endif
+         else ! update_ocn_f = false
+            if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
+               vi0tmp = fnew*dt / (rhoi*Lfresh)
+               dfresh = -rhoi*(vi0new(ij) - vi0tmp)/dt
+               dfsalt = ice_ref_salinity*p001*dfresh
+               fresh(i,j) = fresh(i,j) + dfresh
+               fsalt(i,j) = fsalt(i,j) + dfsalt
+            ! elseif ktherm==1 do nothing
+            endif
          endif
 
       !-----------------------------------------------------------------