Add saltflux_option

cicecore/cicedynB/analysis/ice_diagnostics.F90

@@ -86,6 +86,8 @@ module ice_diagnostics
          totms            , & ! total ice/snow water mass (sh)
          totmin           , & ! total ice water mass (nh)
          totmis           , & ! total ice water mass (sh)
+         totsn            , & ! total salt mass (nh)
+         totss            , & ! total salt mass (sh)
          toten            , & ! total ice/snow energy (J)
          totes                ! total ice/snow energy (J)
 
@@ -154,14 +156,16 @@ subroutine runtime_diags (dt)
          rhofresh, lfresh, lvap, ice_ref_salinity, Tffresh
 
       character (len=char_len) :: &
-         snwredist
+         snwredist, saltflux_option
 
       ! hemispheric state quantities
       real (kind=dbl_kind) :: &
          umaxn,   hmaxn,   shmaxn,    arean,   snwmxn, extentn, shmaxnt, &
          umaxs,   hmaxs,   shmaxs,    areas,   snwmxs, extents, shmaxst, &
          etotn,   mtotn,   micen,     msnwn,   pmaxn,  ketotn, &
          etots,   mtots,   mices,     msnws,   pmaxs,  ketots, &
+         stotn, &
+         stots, &
          urmsn,   albtotn, arean_alb, mpndn,   ptotn,  spondn, &
          urmss,   albtots, areas_alb, mpnds,   ptots,  sponds 
 
@@ -217,7 +221,7 @@ subroutine runtime_diags (dt)
            awtvdr_out=awtvdr, awtidr_out=awtidr, awtvdf_out=awtvdf, awtidf_out=awtidf, &
            rhofresh_out=rhofresh, lfresh_out=lfresh, lvap_out=lvap, &
            ice_ref_salinity_out=ice_ref_salinity,snwredist_out=snwredist, &
-           snwgrain_out=snwgrain)
+           snwgrain_out=snwgrain, saltflux_option_out=saltflux_option)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
          file=__FILE__, line=__LINE__)
@@ -443,6 +447,15 @@ subroutine runtime_diags (dt)
          etots = global_sum(work1, distrb_info, &
                             field_loc_center, tareas)
 
+         ! total salt volume
+         call total_salt   (work2)
+
+         stotn = global_sum(work2, distrb_info, &
+                            field_loc_center, tarean)
+         stots = global_sum(work2, distrb_info, &
+                            field_loc_center, tareas)
+
+
       !-----------------------------------------------------------------
       ! various fluxes
       !-----------------------------------------------------------------
@@ -693,12 +706,22 @@ subroutine runtime_diags (dt)
          swerrs = (fswnets - fswdns) / (fswnets - c1)
 
          ! salt mass
-         msltn = micen*ice_ref_salinity*p001
-         mslts = mices*ice_ref_salinity*p001
+         if (saltflux_option == 'prognostic') then
+            ! compute the total salt mass
+            msltn = stotn*rhoi*p001
+            mslts = stots*rhoi*p001
+
+            ! change in salt mass
+            delmsltn = rhoi*(stotn-totsn)*p001
+            delmslts = rhoi*(stots-totss)*p001
+         else
+            msltn = micen*ice_ref_salinity*p001
+            mslts = mices*ice_ref_salinity*p001
 
-         ! change in salt mass
-         delmsltn = delmxn*ice_ref_salinity*p001
-         delmslts = delmxs*ice_ref_salinity*p001
+            ! change in salt mass
+            delmsltn = delmxn*ice_ref_salinity*p001
+            delmslts = delmxs*ice_ref_salinity*p001
+         endif
 
          ! salt error
          serrn = (sfsaltn + delmsltn) / (msltn + c1)
@@ -1163,7 +1186,7 @@ subroutine init_mass_diags
          rhoi, rhos, rhofresh
 
       real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
-         work1
+         work1, work2
 
       character(len=*), parameter :: subname = '(init_mass_diags)'
 
@@ -1198,6 +1221,12 @@ subroutine init_mass_diags
       toten = global_sum(work1, distrb_info, field_loc_center, tarean)
       totes = global_sum(work1, distrb_info, field_loc_center, tareas)
 
+      ! north/south salt
+      call total_salt (work2)
+      totsn = global_sum(work2, distrb_info, field_loc_center, tarean)
+      totss = global_sum(work2, distrb_info, field_loc_center, tareas)
+
+
       if (print_points) then
          do n = 1, npnt
             if (my_task == pmloc(n)) then

cicecore/cicedynB/analysis/ice_history.F90

@@ -1790,12 +1790,13 @@ subroutine accum_hist (dt)
 
       real (kind=dbl_kind) :: awtvdr, awtidr, awtvdf, awtidf, puny, secday, rad_to_deg
       real (kind=dbl_kind) :: Tffresh, rhoi, rhos, rhow, ice_ref_salinity
-      real (kind=dbl_kind) :: rho_ice, rho_ocn, Tice, Sbr, phi, rhob, dfresh, dfsalt
+      real (kind=dbl_kind) :: rho_ice, rho_ocn, Tice, Sbr, phi, rhob, dfresh, dfsalt, sicen
       logical (kind=log_kind) :: formdrag, skl_bgc
       logical (kind=log_kind) :: tr_pond, tr_aero, tr_brine, tr_snow
       integer (kind=int_kind) :: ktherm
       integer (kind=int_kind) :: nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY, nt_Tsfc, &
                                  nt_alvl, nt_vlvl
+      character (len=char_len) :: saltflux_option
 
       type (block) :: &
          this_block           ! block information for current block
@@ -1807,6 +1808,7 @@ subroutine accum_hist (dt)
       call icepack_query_parameters(Tffresh_out=Tffresh, rhoi_out=rhoi, rhos_out=rhos, &
            rhow_out=rhow, ice_ref_salinity_out=ice_ref_salinity)
       call icepack_query_parameters(formdrag_out=formdrag, skl_bgc_out=skl_bgc, ktherm_out=ktherm)
+      call icepack_query_parameters(saltflux_option_out=saltflux_option)
       call icepack_query_tracer_flags(tr_pond_out=tr_pond, tr_aero_out=tr_aero, &
            tr_brine_out=tr_brine, tr_snow_out=tr_snow)
       call icepack_query_tracer_indices(nt_sice_out=nt_sice, nt_qice_out=nt_qice, &
@@ -1892,7 +1894,7 @@ subroutine accum_hist (dt)
 
 ! MHRI: CHECK THIS OMP ... Maybe ok after "dfresh,dfsalt" added
       !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
-      !$OMP             k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob,dfresh,dfsalt, &
+      !$OMP             k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob,dfresh,dfsalt,sicen, &
       !$OMP             worka,workb,worka3,Tinz4d,Sinz4d,Tsnz4d)
 
       do iblk = 1, nblocks
@@ -2750,7 +2752,16 @@ subroutine accum_hist (dt)
                     dfresh = -rhoi*frazil(i,j,iblk)/dt 
                  endif
                  endif
-                 dfsalt = ice_ref_salinity*p001*dfresh
+                 if (saltflux_option == 'prognostic') then
+                    sicen = c0
+                    do k = 1, nzilyr
+                       sicen = sicen + trcr(i,j,nt_sice+k-1,iblk)*vice(i,j,iblk) &
+                                     / real(nzilyr,kind=dbl_kind)
+                    enddo
+                    dfsalt = sicen*p001*dfresh
+                 else
+                    dfsalt = ice_ref_salinity*p001*dfresh
+                 endif
                  worka(i,j) = aice(i,j,iblk)*(fsalt(i,j,iblk)+dfsalt)
               endif
            enddo

cicecore/cicedynB/general/ice_init.F90

@@ -141,7 +141,7 @@ subroutine input_data
         kitd, kcatbound, ktransport
 
       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
-        tfrz_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
+        tfrz_option, saltflux_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
 
       logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, wave_spec, &
                                  sw_redist, calc_dragio, use_smliq_pnd, snwgrain
@@ -247,6 +247,7 @@ subroutine input_data
         highfreq,       natmiter,        atmiter_conv,  calc_dragio,    &
         ustar_min,      emissivity,      iceruf,        iceruf_ocn,     &
         fbot_xfer_type, update_ocn_f,    l_mpond_fresh, tfrz_option,    &
+        saltflux_option, &
         oceanmixed_ice, restore_ice,     restore_ocn,   trestore,       &
         precip_units,   default_season,  wave_spec_type,nfreq,          &
         atm_data_type,  ocn_data_type,   bgc_data_type, fe_data_type,   &
@@ -466,6 +467,7 @@ subroutine input_data
       precip_units    = 'mks'     ! 'mm_per_month' or
                                   ! 'mm_per_sec' = 'mks' = kg/m^2 s
       tfrz_option     = 'mushy'   ! freezing temp formulation
+      saltflux_option = '4psu'    ! saltflux calculation
       oceanmixed_ice  = .false.   ! if true, use internal ocean mixed layer
       wave_spec_type  = 'none'    ! type of wave spectrum forcing
       nfreq           = 25        ! number of wave frequencies
@@ -929,6 +931,7 @@ subroutine input_data
       call broadcast_scalar(wave_spec_file,       master_task)
       call broadcast_scalar(nfreq,                master_task)
       call broadcast_scalar(tfrz_option,          master_task)
+      call broadcast_scalar(saltflux_option,      master_task)
       call broadcast_scalar(ocn_data_format,      master_task)
       call broadcast_scalar(bgc_data_type,        master_task)
       call broadcast_scalar(fe_data_type,         master_task)
@@ -2164,6 +2167,7 @@ subroutine input_data
          wave_spec_type_in = wave_spec_type, &
          wave_spec_in=wave_spec, nfreq_in=nfreq, &
          tfrz_option_in=tfrz_option, kalg_in=kalg, fbot_xfer_type_in=fbot_xfer_type, &
+         saltflux_option_in=saltflux_option, &
          Pstar_in=Pstar, Cstar_in=Cstar, iceruf_in=iceruf, iceruf_ocn_in=iceruf_ocn, calc_dragio_in=calc_dragio, &
          windmin_in=windmin, drhosdwind_in=drhosdwind, &
          rsnw_fall_in=rsnw_fall, rsnw_tmax_in=rsnw_tmax, rhosnew_in=rhosnew, &