From f5e093f5148554674079d5c7fc0702a41b81f744 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Tue, 22 Aug 2023 10:11:49 -0700 Subject: [PATCH 1/2] Deprecate zsalinity (#448) * Deprecate zsalinity mostly with ifdef and comments for testing * Deprecate zsalinity, remove code * Update zsalinity warning messages Update interface documentation --- columnphysics/icepack_brine.F90 | 261 +--- columnphysics/icepack_intfc.F90 | 2 +- columnphysics/icepack_itd.F90 | 48 +- columnphysics/icepack_mechred.F90 | 7 +- columnphysics/icepack_parameters.F90 | 6 +- columnphysics/icepack_therm_bl99.F90 | 4 +- columnphysics/icepack_therm_itd.F90 | 32 +- columnphysics/icepack_tracers.F90 | 8 +- columnphysics/icepack_zbgc.F90 | 217 +-- columnphysics/icepack_zbgc_shared.F90 | 7 - columnphysics/icepack_zsalinity.F90 | 1158 ----------------- configuration/driver/icedrv_InitMod.F90 | 4 +- configuration/driver/icedrv_RunMod.F90 | 16 +- configuration/driver/icedrv_arrays_column.F90 | 16 - .../driver/icedrv_diagnostics_bgc.F90 | 178 +-- configuration/driver/icedrv_domain_size.F90 | 4 +- configuration/driver/icedrv_flux.F90 | 8 +- configuration/driver/icedrv_init_column.F90 | 98 +- configuration/driver/icedrv_restart.F90 | 7 +- configuration/driver/icedrv_restart_bgc.F90 | 65 +- configuration/driver/icedrv_step.F90 | 19 +- configuration/scripts/icepack.build | 2 +- configuration/scripts/icepack.settings | 2 - .../scripts/options/set_env.bgcispol | 2 - configuration/scripts/options/set_env.bgcnice | 2 - .../scripts/options/set_env.bgcsklnice | 2 - configuration/scripts/options/set_env.zsal | 10 - configuration/scripts/options/set_nml.zsal | 5 - .../developer_guide/dg_adding_diagnostics.rst | 2 +- doc/source/developer_guide/dg_col_phys.rst | 1 - doc/source/user_guide/interfaces.include | 37 +- doc/source/user_guide/lg_protocols.rst | 2 - doc/source/user_guide/lg_sequence.rst | 5 +- doc/source/user_guide/ug_case_settings.rst | 5 +- 34 files changed, 139 insertions(+), 2103 deletions(-) delete mode 100644 columnphysics/icepack_zsalinity.F90 delete mode 100644 configuration/scripts/options/set_env.zsal delete mode 100644 configuration/scripts/options/set_nml.zsal diff --git a/columnphysics/icepack_brine.F90 b/columnphysics/icepack_brine.F90 index 952b326f5..17a4bcec6 100644 --- a/columnphysics/icepack_brine.F90 +++ b/columnphysics/icepack_brine.F90 @@ -11,7 +11,7 @@ module icepack_brine use icepack_parameters, only: gravit, rhoi, rhow, rhos, depressT use icepack_parameters, only: salt_loss, min_salin, rhosi use icepack_parameters, only: dts_b, l_sk - use icepack_tracers, only: ntrcr, nt_qice, nt_sice, nt_bgc_S + use icepack_tracers, only: ntrcr, nt_qice, nt_sice use icepack_tracers, only: nt_Tsfc use icepack_zbgc_shared, only: k_o, exp_h, Dm, Ra_c, viscos_dynamic, thinS use icepack_zbgc_shared, only: remap_zbgc @@ -27,10 +27,9 @@ module icepack_brine public :: preflushing_changes, & compute_microS_mushy, & update_hbrine, & - compute_microS, & calculate_drho, & icepack_init_hbrine, & - icepack_init_zsalinity + icepack_init_zsalinity ! deprecated real (kind=dbl_kind), parameter :: & maxhbr = 1.25_dbl_kind , & ! brine overflows if hbr > maxhbr*hin @@ -131,7 +130,6 @@ subroutine preflushing_changes (aicen, vicen, vsnon, & end subroutine preflushing_changes !======================================================================= - ! Computes ice microstructural properties for updating hbrine ! ! NOTE: This subroutine uses thermosaline_vertical output to compute @@ -558,239 +556,6 @@ subroutine update_hbrine (meltt, & end subroutine update_hbrine -!======================================================================= -! -! Computes ice microstructural properties for zbgc and zsalinity -! -! NOTE: In this subroutine, trcrn(nt_fbri) is the volume fraction of ice with -! dynamic salinity or the height ratio == hbr/vicen*aicen, where hbr is the -! height of the brine surface relative to the bottom of the ice. -! This volume fraction -! may be > 1 in which case there is brine above the ice surface (meltponds). -! - subroutine compute_microS (n_cat, nilyr, nblyr, & - bgrid, cgrid, igrid, & - trcrn, hice_old, & - hbr_old, sss, sst, & - bTin, iTin, bphin, & - kperm, bphi_min, & - Rayleigh_criteria, firstice, & - bSin, brine_sal, & - brine_rho, iphin, ibrine_rho, & - ibrine_sal, sice_rho, sloss) - - integer (kind=int_kind), intent(in) :: & - n_cat , & ! ice category - nilyr , & ! number of ice layers - nblyr ! number of bio layers - - real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & - bgrid ! biology nondimensional vertical grid points - - real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & - igrid ! biology vertical interface points - - real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: & - cgrid ! CICE vertical coordinate - - real (kind=dbl_kind), intent(in) :: & - hice_old , & ! previous timestep ice height (m) - sss , & ! ocean salinity (ppt) - sst ! ocean temperature (oC) - - real (kind=dbl_kind), dimension(ntrcr), intent(inout) :: & - trcrn - - real (kind=dbl_kind), intent(inout) :: & - hbr_old , & ! old brine height - sice_rho ! average ice density - - real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: & - bTin , & ! Temperature of ice layers on bio grid for history file (^oC) - bphin , & ! Porosity of layers - brine_sal , & ! equilibrium brine salinity (ppt) - brine_rho ! Internal brine density (kg/m^3) - - real (kind=dbl_kind), dimension (nblyr+2), intent(out) :: & - bSin - - real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: & - iTin ! Temperature on the interface grid - - real (kind=dbl_kind), intent(out) :: & - bphi_min , & ! surface porosity - kperm , & ! average ice permeability (m^2) - sloss ! (g/m^2) salt from brine runoff for hbr > maxhbr*hin - - logical (kind=log_kind), intent(inout) :: & - Rayleigh_criteria ! .true. if ice exceeded a minimum thickness hin >= Ra_c - - logical (kind=log_kind), intent(in) :: & - firstice ! .true. if ice is newly formed - - real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: & - iphin , & ! porosity on the igrid - ibrine_rho , & ! brine rho on interface - ibrine_sal ! brine sal on interface - - ! local variables - - integer (kind=int_kind) :: & - k ! vertical biology layer index - - real (kind=dbl_kind) :: & - surface_S , & ! salinity of ice above hin > hbr - hinc_old ! ice thickness (cell quantity) before current melt/growth (m) - -! logical (kind=log_kind) :: & -! Rayleigh ! .true. if ice exceeded a minimum thickness hin >= Ra_c - - real (kind=dbl_kind), dimension (ntrcr+2) :: & - trtmp0 , & ! temporary, remapped tracers - trtmp ! temporary, remapped tracers - - real (kind=dbl_kind) :: & - Tmlts ! melting temperature - - character(len=*),parameter :: subname='(compute_microS)' - - !----------------------------------------------------------------- - ! Initialize - !----------------------------------------------------------------- - - sloss = c0 - bTin(:) = c0 - bSin(:) = c0 - - trtmp(:) = c0 - surface_S = min_salin - - hinc_old = hice_old - - !----------------------------------------------------------------- - ! Rayleigh condition for salinity and bgc: - ! Implemented as a minimum thickness criteria for category 1 ice only. - ! When hin >= Ra_c (m), pressure flow is allowed. - ! Turn off by putting Ra_c = 0 in ice_in namelist. - !----------------------------------------------------------------- - -! Rayleigh = .true. -! if (n_cat == 1 .AND. hbr_old < Ra_c) then -! Rayleigh = Rayleigh_criteria ! only category 1 ice can be false -! endif - - !----------------------------------------------------------------- - ! Define ice salinity on Sin - !----------------------------------------------------------------- - - if (firstice) then - - do k = 1, nilyr - trcrn(nt_sice+k-1) = sss*salt_loss - enddo - - call remap_zbgc(nilyr, & - nt_sice, & - trcrn, trtmp, & - 0, nblyr, & - hinc_old, hinc_old, & - cgrid(2:nilyr+1), & - bgrid(2:nblyr+1), surface_S ) - if (icepack_warnings_aborted(subname)) return - - do k = 1, nblyr - trcrn(nt_bgc_S+k-1) = max(min_salin,trtmp(nt_sice+k-1)) - bSin(k+1) = max(min_salin,trcrn(nt_bgc_S+k-1)) - if (trcrn(nt_bgc_S+k-1) < min_salin-puny) & - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - enddo ! k - - bSin(1) = bSin(2) - bSin(nblyr+2) = sss - - elseif (hbr_old > maxhbr*hice_old) then - - call remap_zbgc(nblyr, & - nt_bgc_S, & - trcrn, trtmp, & - 0, nblyr, & - hbr_old, & - maxhbr*hinc_old, & - bgrid(2:nblyr+1), & - bgrid(2:nblyr+1), surface_S ) - if (icepack_warnings_aborted(subname)) return - - do k = 1, nblyr - bSin(k+1) = max(min_salin,trtmp(nt_bgc_S+k-1)) - sloss = sloss + rhosi*(hbr_old*trcrn(nt_bgc_S+k-1) & - - maxhbr*hice_old*bSin(k+1))*(igrid(k+1)-igrid(k)) - trcrn(nt_bgc_S+k-1) = bSin(k+1) - if (trcrn(nt_bgc_S+k-1) < min_salin-puny) & - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - enddo ! k - - bSin(1) = bSin(2) - bSin(nblyr+2) = sss - hbr_old = maxhbr*hinc_old - - else ! old, thin ice - - do k = 1, nblyr - trcrn(nt_bgc_S+k-1) = max(min_salin,trcrn(nt_bgc_S+k-1)) - bSin (k+1) = trcrn(nt_bgc_S+k-1) - enddo ! k - - bSin (1) = bSin(2) - bSin (nblyr+2) = sss - - endif ! ice type - - !----------------------------------------------------------------- - ! sea ice temperature for bio grid - !----------------------------------------------------------------- - - do k = 1, nilyr - Tmlts = -trcrn(nt_sice+k-1)*depressT - trtmp0(nt_qice+k-1) = calculate_Tin_from_qin(trcrn(nt_qice+k-1),Tmlts) - enddo ! k - - trtmp(:) = c0 - - ! CICE to Bio: remap temperatures - call remap_zbgc (nilyr, nt_qice, & - trtmp0(1:ntrcr), trtmp, & - 0, nblyr, & - hinc_old, hbr_old, & - cgrid(2:nilyr+1), & - bgrid(2:nblyr+1), surface_S ) - if (icepack_warnings_aborted(subname)) return - - do k = 1, nblyr - Tmlts = -bSin(k+1) * depressT - bTin (k+1) = min(Tmlts,trtmp(nt_qice+k-1)) - enddo !k - - Tmlts = -min_salin* depressT - bTin (1) = min(Tmlts,(bTin(2) + trcrn(nt_Tsfc))*p5) - Tmlts = -bSin(nblyr+2)* depressT - bTin (nblyr+2) = sst - - !----------------------------------------------------------------- - ! Define ice multiphase structure - !----------------------------------------------------------------- - - call prepare_hbrine (nblyr, & - bSin, bTin, iTin, & - brine_sal, brine_rho, & - ibrine_sal, ibrine_rho, & - sice_rho, & - bphin, iphin, & - kperm, bphi_min, & - igrid, sss) - if (icepack_warnings_aborted(subname)) return - - end subroutine compute_microS - !========================================================================================== ! ! Find density difference about interface grid points @@ -1000,7 +765,8 @@ end subroutine icepack_init_hbrine !======================================================================= !autodocument_start icepack_init_zsalinity -! Initialize zSalinity +! **DEPRECATED**, all code removed +! Interface provided for backwards compatibility subroutine icepack_init_zsalinity(nblyr,ntrcr_o, Rayleigh_criteria, & Rayleigh_real, trcrn_bgc, nt_bgc_S, ncat, sss) @@ -1027,27 +793,12 @@ subroutine icepack_init_zsalinity(nblyr,ntrcr_o, Rayleigh_criteria, & ! local variables - integer (kind=int_kind) :: & - k, n - character(len=*),parameter :: subname='(icepack_init_zsalinity)' - if (nblyr .LE. 7) then - dts_b = 300.0_dbl_kind - else - dts_b = 50.0_dbl_kind - endif - - Rayleigh_criteria = .false. ! no ice initial condition - Rayleigh_real = c0 - do n = 1,ncat - do k = 1,nblyr - trcrn_bgc(nt_bgc_S+k-1-ntrcr_o,n) = sss*salt_loss - enddo ! k - enddo ! n + call icepack_warnings_add(subname//' DEPRECATED, do not use') +! call icepack_warnings_setabort(.true.,__FILE__,__LINE__) end subroutine icepack_init_zsalinity - !======================================================================= end module icepack_brine diff --git a/columnphysics/icepack_intfc.F90 b/columnphysics/icepack_intfc.F90 index 9578e83c3..4a7d466dc 100644 --- a/columnphysics/icepack_intfc.F90 +++ b/columnphysics/icepack_intfc.F90 @@ -89,7 +89,7 @@ module icepack_intfc use icepack_shortwave, only: icepack_step_radiation use icepack_brine, only: icepack_init_hbrine - use icepack_brine, only: icepack_init_zsalinity + use icepack_brine, only: icepack_init_zsalinity ! deprecated use icepack_zbgc , only: icepack_init_bgc use icepack_zbgc , only: icepack_init_zbgc diff --git a/columnphysics/icepack_itd.F90 b/columnphysics/icepack_itd.F90 index 00f7cce5b..d8b0875aa 100644 --- a/columnphysics/icepack_itd.F90 +++ b/columnphysics/icepack_itd.F90 @@ -30,10 +30,10 @@ module icepack_itd use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, Tocnfrz, rhoi use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall, rhosnew use icepack_tracers, only: nt_Tsfc, nt_qice, nt_qsno, nt_aero, nt_isosno, nt_isoice - use icepack_tracers, only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, nt_bgc_S, bio_index + use icepack_tracers, only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, bio_index use icepack_tracers, only: n_iso, tr_iso, nt_smice, nt_rsnw, nt_rhos, nt_sice use icepack_tracers, only: icepack_compute_tracers - use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers + use icepack_parameters, only: skl_bgc, z_tracers use icepack_parameters, only: kcatbound, kitd, saltflux_option, snwgrain, snwredist use icepack_therm_shared, only: Tmin, hi_min use icepack_warnings, only: warnstr, icepack_warnings_add @@ -769,7 +769,6 @@ subroutine cleanup_itd (dt, ntrcr, & fpond, fresh, & fsalt, fhocn, & faero_ocn, fiso_ocn, & - fzsal, & flux_bio, limit_aice_in) integer (kind=int_kind), intent(in) :: & @@ -823,8 +822,7 @@ subroutine cleanup_itd (dt, ntrcr, & fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) - fhocn , & ! net heat flux to ocean (W/m^2) - fzsal ! net salt flux to ocean from zsalinity (kg/m^2/s) + fhocn ! net heat flux to ocean (W/m^2) real (kind=dbl_kind), dimension (:), intent(inout), optional :: & flux_bio ! net tracer flux to ocean from biology (mmol/m^2/s) @@ -849,8 +847,7 @@ subroutine cleanup_itd (dt, ntrcr, & dfpond , & ! zapped pond water flux (kg/m^2/s) dfresh , & ! zapped fresh water flux (kg/m^2/s) dfsalt , & ! zapped salt flux (kg/m^2/s) - dfhocn , & ! zapped energy flux ( W/m^2) - dfzsal ! zapped salt flux for zsalinity (kg/m^2/s) + dfhocn ! zapped energy flux ( W/m^2) real (kind=dbl_kind), dimension (n_aero) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) @@ -883,7 +880,6 @@ subroutine cleanup_itd (dt, ntrcr, & dfaero_ocn(:) = c0 dfiso_ocn (:) = c0 dflux_bio (:) = c0 - dfzsal = c0 !----------------------------------------------------------------- ! Compute total ice area. @@ -950,7 +946,7 @@ subroutine cleanup_itd (dt, ntrcr, & tr_aero, & tr_pond_topo, & first_ice, nbtrcr, & - dfzsal, dflux_bio ) + dflux_bio ) if (icepack_warnings_aborted(subname)) then write(warnstr,*) subname, 'aice:', aice @@ -1008,8 +1004,6 @@ subroutine cleanup_itd (dt, ntrcr, & flux_bio (it) = flux_bio(it) + dflux_bio(it) enddo endif - if (present(fzsal)) & - fzsal = fzsal + dfzsal end subroutine cleanup_itd @@ -1033,9 +1027,9 @@ subroutine zap_small_areas (dt, ntrcr, & dfhocn, & dfaero_ocn, dfiso_ocn, & tr_aero, & - tr_pond_topo, & + tr_pond_topo, & first_ice, nbtrcr, & - dfzsal, dflux_bio ) + dflux_bio ) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -1065,8 +1059,7 @@ subroutine zap_small_areas (dt, ntrcr, & dfpond , & ! zapped pond water flux (kg/m^2/s) dfresh , & ! zapped fresh water flux (kg/m^2/s) dfsalt , & ! zapped salt flux (kg/m^2/s) - dfhocn , & ! zapped energy flux ( W/m^2) - dfzsal ! zapped salt flux from zsalinity(kg/m^2/s) + dfhocn ! zapped energy flux ( W/m^2) real (kind=dbl_kind), dimension (:), intent(inout) :: & dfaero_ocn ! zapped aerosol flux (kg/m^2/s) @@ -1099,7 +1092,6 @@ subroutine zap_small_areas (dt, ntrcr, & !----------------------------------------------------------------- ! I. Zap categories with very small areas. !----------------------------------------------------------------- - dfzsal = c0 do n = 1, ncat @@ -1139,15 +1131,6 @@ subroutine zap_small_areas (dt, ntrcr, & enddo endif - if (solve_zsal) then - do it = 1, nblyr - xtmp = rhosi*trcrn(nt_fbri,n)*vicen(n)*p001& - *trcrn(nt_bgc_S+it-1,n)/ & - real(nblyr,kind=dbl_kind)/dt - dfzsal = dfzsal + xtmp - enddo ! n - endif - if (skl_bgc .and. nbtrcr > 0) then blevels = 1 bvol(1) = aicen(n)*sk_l @@ -1330,21 +1313,6 @@ subroutine zap_small_areas (dt, ntrcr, & endif dfsalt = dfsalt + xtmp - if (solve_zsal) then - do k = 1,nblyr - xtmp = rhosi*trcrn(nt_fbri,n)*vicen(n)*p001& - /real(nblyr,kind=dbl_kind)*trcrn(nt_bgc_S+k-1,n) & - * (aice-c1)/aice /dt - dfzsal = dfzsal + xtmp - enddo - - if (vicen(n) > vicen(n)*trcrn(nt_fbri,n)) then - xtmp = (vicen(n)-vicen(n)*trcrn(nt_fbri,n))*(aice-c1)/& - aice*p001*rhosi*min_salin/dt - dfzsal = dfzsal + xtmp - endif - endif ! solve_zsal - aicen(n) = aicen(n) * (c1/aice) vicen(n) = vicen(n) * (c1/aice) vsnon(n) = vsnon(n) * (c1/aice) diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90 index 28b0330c9..929ec648b 100644 --- a/columnphysics/icepack_mechred.F90 +++ b/columnphysics/icepack_mechred.F90 @@ -1777,8 +1777,10 @@ subroutine icepack_step_ridge (dt, ndtd, & fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) - fhocn , & ! net heat flux to ocean (W/m^2) - fzsal ! zsalinity flux to ocean(kg/m^2/s) + fhocn ! net heat flux to ocean (W/m^2) + + real (kind=dbl_kind), intent(inout), optional :: & + fzsal ! zsalinity flux to ocean(kg/m^2/s) (deprecated) real (kind=dbl_kind), intent(inout), optional :: & closing ! rate of closing due to divergence/shear (1/s) @@ -1897,7 +1899,6 @@ subroutine icepack_step_ridge (dt, ndtd, & fpond, fresh, & fsalt, fhocn, & faero_ocn, fiso_ocn, & - fzsal, & flux_bio) if (icepack_warnings_aborted(subname)) return diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index e59067129..c0bb453c7 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -1085,16 +1085,12 @@ subroutine icepack_init_parameters( & if (present(modal_aero_in) ) modal_aero = modal_aero_in if (present(conserv_check_in) ) conserv_check = conserv_check_in if (present(skl_bgc_in) ) skl_bgc = skl_bgc_in -#ifdef UNDEPRECATE_ZSAL - if (present(solve_zsal_in) ) solve_zsal = solve_zsal_in -#else if (present(solve_zsal_in)) then + call icepack_warnings_add(subname//' WARNING: zsalinity is deprecated') if (solve_zsal_in) then - call icepack_warnings_add(subname//' zsalinity is being deprecated') call icepack_warnings_setabort(.true.,__FILE__,__LINE__) endif endif -#endif if (present(grid_o_in) ) grid_o = grid_o_in if (present(l_sk_in) ) l_sk = l_sk_in if (present(initbio_frac_in) ) initbio_frac = initbio_frac_in diff --git a/columnphysics/icepack_therm_bl99.F90 b/columnphysics/icepack_therm_bl99.F90 index 7a090ef91..03d46ee6e 100644 --- a/columnphysics/icepack_therm_bl99.F90 +++ b/columnphysics/icepack_therm_bl99.F90 @@ -17,7 +17,7 @@ module icepack_therm_bl99 use icepack_parameters, only: p01 #endif use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice - use icepack_parameters, only: conduct, calc_Tsfc, solve_zsal + use icepack_parameters, only: conduct, calc_Tsfc use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted @@ -266,8 +266,6 @@ subroutine temperature_changes (dt, & if (sw_redist) then - if (solve_zsal) sw_dtemp = p1 ! lower tolerance with dynamic salinity - do k = 1, nilyr Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 5cfb0b117..405427af3 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -26,7 +26,7 @@ module icepack_therm_itd use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss use icepack_parameters, only: rhosi, conserv_check, rhosmin, snwredist use icepack_parameters, only: kitd, ktherm - use icepack_parameters, only: z_tracers, solve_zsal, hfrazilmin + use icepack_parameters, only: z_tracers, hfrazilmin use icepack_parameters, only: saltflux_option use icepack_parameters, only: icepack_chkoptargflag @@ -887,7 +887,7 @@ subroutine lateral_melt (dt, ncat, & fside, wlat, & aicen, vicen, & vsnon, trcrn, & - fzsal, flux_bio, & + flux_bio, & nbtrcr, nblyr, & nfsd, d_afsd_latm,& floe_rad_c, floe_binwidth) @@ -928,8 +928,7 @@ subroutine lateral_melt (dt, ncat, & fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) fhocn , & ! net heat flux to ocean (W/m^2) - meltl , & ! lateral ice melt (m/step-->cm/day) - fzsal ! salt flux from zsalinity (kg/m2/s) + meltl ! lateral ice melt (m/step-->cm/day) real (kind=dbl_kind), dimension(nbtrcr), intent(inout) :: & flux_bio ! biology tracer flux from layer bgc (mmol/m^2/s) @@ -1239,11 +1238,11 @@ subroutine lateral_melt (dt, ncat, & enddo ! n - if (solve_zsal .or. z_tracers) & + if (z_tracers) & call lateral_melt_bgc(dt, & ncat, nblyr, & rside, vicen_init, & !echmod: use rsiden - trcrn, fzsal, & + trcrn, & flux_bio, nbtrcr) if (icepack_warnings_aborted(subname)) return @@ -1314,7 +1313,7 @@ subroutine add_new_ice (ncat, nilyr, & dSin0_frazil, & bgrid, cgrid, igrid, & nbtrcr, flux_bio, & - ocean_bio, fzsal, & + ocean_bio, & frazil_diag, & fiso_ocn, & HDO_ocn, H2_16O_ocn, & @@ -1402,10 +1401,6 @@ subroutine add_new_ice (ncat, nilyr, & real (kind=dbl_kind), dimension (:), intent(in) :: & ocean_bio ! ocean concentration of biological tracer - ! zsalinity - real (kind=dbl_kind), intent(inout) :: & - fzsal ! salt flux to ocean from zsalinity (kg/m^2/s) - ! water isotopes real (kind=dbl_kind), dimension(:), intent(inout), optional :: & @@ -1610,7 +1605,6 @@ subroutine add_new_ice (ncat, nilyr, & ! history diagnostics frazil = vi0new - if (solve_zsal) fzsal = fzsal - rhosi*vi0new/dt*p001*sss*salt_loss if (present(frz_onset) .and. present(yday)) then if (frazil > puny .and. frz_onset < puny) frz_onset = yday @@ -1796,7 +1790,6 @@ subroutine add_new_ice (ncat, nilyr, & trcrn(nt_qice+k-1,n) = & (trcrn(nt_qice+k-1,n)*vtmp + qi0new*vsurp) / vicen(n) ! salinity - if (.not. solve_zsal) & trcrn(nt_sice+k-1,n) = & (trcrn(nt_sice+k-1,n)*vtmp + Sprofile(k)*vsurp) / vicen(n) endif @@ -1916,7 +1909,6 @@ subroutine add_new_ice (ncat, nilyr, & trcrn(nt_qice+k-1,n) = & (trcrn(nt_qice+k-1,n)*vice1 + qi0new*vin0new(n))/vicen(n) ! salinity - if (.NOT. solve_zsal)& trcrn(nt_sice+k-1,n) = & (trcrn(nt_sice+k-1,n)*vice1 + Sprofile(k)*vin0new(n))/vicen(n) endif @@ -2038,7 +2030,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & Tf , & ! freezing temperature (C) sss , & ! sea surface salinity (ppt) rside , & ! fraction of ice that melts laterally - frzmlt ! freezing/melting potential (W/m^2) + frzmlt ! freezing/melting potential (W/m^2) integer (kind=int_kind), dimension (:), intent(in) :: & trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon @@ -2073,11 +2065,13 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) fhocn , & ! net heat flux to ocean (W/m^2) - fzsal , & ! salt flux to ocean from zsalinity (kg/m^2/s) meltl , & ! lateral ice melt (m/step-->cm/day) frazil , & ! frazil ice growth (m/step-->cm/day) frazil_diag ! frazil ice growth diagnostic (m/step-->cm/day) + real (kind=dbl_kind), intent(inout), optional :: & + fzsal ! salt flux to ocean from zsalinity (kg/m^2/s) (deprecated) + real (kind=dbl_kind), intent(in), optional :: & wlat ! lateral melt rate (m/s) @@ -2250,7 +2244,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & dSin0_frazil, bgrid, & cgrid, igrid, & nbtrcr, flux_bio, & - ocean_bio, fzsal, & + ocean_bio, & frazil_diag, fiso_ocn, & HDO_ocn, H2_16O_ocn, & H2_18O_ocn, & @@ -2276,7 +2270,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & fside, wlat, & aicen, vicen, & vsnon, trcrn, & - fzsal, flux_bio, & + flux_bio, & nbtrcr, nblyr, & nfsd, d_afsd_latm, & floe_rad_c,floe_binwidth) @@ -2325,7 +2319,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, & fpond, fresh, & fsalt, fhocn, & faero_ocn, fiso_ocn, & - fzsal, flux_bio) + flux_bio ) if (icepack_warnings_aborted(subname)) return first_call = .false. diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90 index 11b48b8de..2012cb264 100644 --- a/columnphysics/icepack_tracers.F90 +++ b/columnphysics/icepack_tracers.F90 @@ -97,7 +97,7 @@ module icepack_tracers nt_bgc_PON = 0, & ! zooplankton and detritus nt_bgc_hum = 0, & ! humic material nt_zbgc_frac = 0, & ! fraction of tracer in the mobile phase - nt_bgc_S = 0 ! Bulk salinity in fraction ice with dynamic salinity (Bio grid) + nt_bgc_S = 0 ! Bulk salinity in fraction ice with dynamic salinity (Bio grid) (deprecated) logical (kind=log_kind), public :: & tr_iage = .false., & ! if .true., use age tracer @@ -434,7 +434,7 @@ subroutine icepack_init_tracer_indices(& nlt_bgc_hum_in,& ! nlt_bgc_PON_in,& ! zooplankton and detritus nt_zbgc_frac_in,&! fraction of tracer in the mobile phase - nt_bgc_S_in, & ! Bulk salinity in fraction ice with dynamic salinity (Bio grid)) + nt_bgc_S_in, & ! (deprecated, was related to zsalinity) nlt_chl_sw_in ! points to total chla in trcrn_sw integer (kind=int_kind), dimension(:), intent(in), optional :: & @@ -796,7 +796,7 @@ subroutine icepack_query_tracer_indices(& nlt_bgc_hum_out,& ! nlt_bgc_PON_out,& ! zooplankton and detritus nt_zbgc_frac_out,&! fraction of tracer in the mobile phase - nt_bgc_S_out, & ! Bulk salinity in fraction ice with dynamic salinity (Bio grid)) + nt_bgc_S_out, & ! (deprecated, was related to zsalinity) nlt_chl_sw_out ! points to total chla in trcrn_sw integer (kind=int_kind), dimension(:), intent(out), optional :: & @@ -955,7 +955,7 @@ subroutine icepack_write_tracer_indices(iounit) write(iounit,*) " nlt_bgc_PON = ",nlt_bgc_PON write(iounit,*) " nlt_chl_sw = ",nlt_chl_sw write(iounit,*) " nt_zbgc_frac = ",nt_zbgc_frac - write(iounit,*) " nt_bgc_S = ",nt_bgc_S + write(iounit,*) " nt_bgc_S = ",nt_bgc_S," (deprecated)" write(iounit,*) " max_nbtrcr = ",max_nbtrcr do k = 1, max_nbtrcr diff --git a/columnphysics/icepack_zbgc.F90 b/columnphysics/icepack_zbgc.F90 index 715a2c514..928b4ca65 100644 --- a/columnphysics/icepack_zbgc.F90 +++ b/columnphysics/icepack_zbgc.F90 @@ -15,10 +15,10 @@ module icepack_zbgc use icepack_parameters, only: op_dep_min, fr_graze_s, fr_graze_e, fr_mort2min, fr_dFe use icepack_parameters, only: k_nitrif, t_iron_conv, max_loss, max_dfe_doc1 use icepack_parameters, only: fr_resp_s, y_sk_DMS, t_sk_conv, t_sk_ox - use icepack_parameters, only: scale_bgc, ktherm, skl_bgc, solve_zsal + use icepack_parameters, only: scale_bgc, ktherm, skl_bgc use icepack_parameters, only: z_tracers, fsal, conserv_check - use icepack_tracers, only: nt_sice, nt_bgc_S, bio_index + use icepack_tracers, only: nt_sice, bio_index use icepack_tracers, only: tr_brine, nt_fbri, nt_qice, nt_Tsfc use icepack_tracers, only: nt_zbgc_frac use icepack_tracers, only: bio_index_o, bio_index @@ -41,13 +41,10 @@ module icepack_zbgc use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted use icepack_brine, only: preflushing_changes, compute_microS_mushy - use icepack_brine, only: update_hbrine, compute_microS + use icepack_brine, only: update_hbrine use icepack_algae, only: zbio, sklbio use icepack_therm_shared, only: calculate_Tin_from_qin use icepack_itd, only: column_sum, column_conservation_check -#ifdef UNDEPRECATE_ZSAL - use icepack_zsalinity, only: zsalinity -#endif implicit none @@ -209,7 +206,7 @@ subroutine add_new_ice_bgc (dt, nblyr, & vtmp, & vsurp, sss, & nilyr, nblyr, & - solve_zsal, bgrid, & + bgrid, & cgrid, & ocean_bio, igrid, & location) @@ -246,16 +243,12 @@ subroutine add_new_ice_bgc (dt, nblyr, & vbri1, & vi0new, sss, & nilyr, nblyr, & - solve_zsal, bgrid, & + bgrid, & cgrid, & ocean_bio, igrid, & location) if (icepack_warnings_aborted(subname)) return - if (solve_zsal .and. vsnon1 .le. c0) then - Tmlts = -trcrn(nt_sice,1)*depressT - trcrn(nt_Tsfc,1) = calculate_Tin_from_qin(trcrn(nt_qice,1),Tmlts) - endif ! solve_zsal endif ! nltrcr > 0 endif ! vi0new > 0 @@ -280,7 +273,7 @@ end subroutine add_new_ice_bgc subroutine lateral_melt_bgc (dt, & ncat, nblyr, & rside, vicen, & - trcrn, fzsal, & + trcrn, & flux_bio, nbltrcr) integer (kind=int_kind), intent(in) :: & @@ -300,9 +293,6 @@ subroutine lateral_melt_bgc (dt, & real (kind=dbl_kind), intent(in) :: & rside ! fraction of ice that melts laterally - real (kind=dbl_kind), intent(inout) :: & - fzsal ! salt flux from layer Salinity (kg/m^2/s) - real (kind=dbl_kind), dimension(:), intent(inout) :: & flux_bio ! biology tracer flux from layer bgc (mmol/m^2/s) @@ -320,16 +310,6 @@ subroutine lateral_melt_bgc (dt, & zspace = c1/(real(nblyr,kind=dbl_kind)) - if (solve_zsal) then - do n = 1, ncat - do k = 1,nblyr - fzsal = fzsal + rhosi*trcrn(nt_fbri,n) & - * vicen(n)*p001*zspace*trcrn(nt_bgc_S+k-1,n) & - * rside/dt - enddo - enddo - endif - do m = 1, nbltrcr do n = 1, ncat do k = 1, nblyr+1 @@ -354,7 +334,7 @@ subroutine adjust_tracer_profile (nbtrcr, dt, ntrcr, & vtmp, & vsurp, sss, & nilyr, nblyr, & - solve_zsal, bgrid, & + bgrid, & cgrid, ocean_bio, & igrid, location) @@ -382,9 +362,6 @@ subroutine adjust_tracer_profile (nbtrcr, dt, ntrcr, & real (kind=dbl_kind), intent(in) :: & vbrin ! fbri*volume per unit area of ice (m) - logical (kind=log_kind), intent(in) :: & - solve_zsal - real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & igrid ! zbio grid @@ -433,22 +410,6 @@ subroutine adjust_tracer_profile (nbtrcr, dt, ntrcr, & hbri = vbrin hbri_old = vtmp - if (solve_zsal) then - top_conc = sss * salt_loss - do k = 1, nblyr - S_stationary(k) = trcrn(nt_bgc_S+k-1)* hbri_old - enddo - call regrid_stationary (S_stationary, hbri_old, & - hbri, dt, & - ntrcr, & - nblyr-1, top_conc, & - bgrid(2:nblyr+1), fluxb ) - if (icepack_warnings_aborted(subname)) return - do k = 1, nblyr - trcrn(nt_bgc_S+k-1) = S_stationary(k)/hbri - trtmp0(nt_sice+k-1) = trcrn(nt_bgc_S+k-1) - enddo - endif ! solve_zsal do m = 1, nbtrcr top_conc = ocean_bio(m)*zbgc_init_frac(m) @@ -466,64 +427,15 @@ subroutine adjust_tracer_profile (nbtrcr, dt, ntrcr, & enddo !k enddo !m - if (solve_zsal) then - if (aicen > c0) then - hinS_new = vbrin/aicen - hin = vicen/aicen - else - hinS_new = c0 - hin = c0 - endif ! aicen - temp_S = min_salin ! bio to cice - call remap_zbgc(nilyr, & - nt_sice, & - trtmp0(1:ntrcr), trtmp, & - 1, nblyr, & - hin, hinS_new, & - cgrid(2:nilyr+1), & - bgrid(2:nblyr+1), temp_S ) - if (icepack_warnings_aborted(subname)) return - do k = 1, nilyr - trcrn(nt_sice+k-1) = trtmp(nt_sice+k-1) - enddo ! k - endif ! solve_zsal - elseif (vbrin > c0) then ! add frazil throughout location == 0 .and. do k = 1, nblyr+1 - if (solve_zsal .and. k < nblyr + 1) then - trcrn(nt_bgc_S+k-1) = (trcrn(nt_bgc_S+k-1) * vtmp & - + sss*salt_loss * vsurp) / vbrin - trtmp0(nt_sice+k-1) = trcrn(nt_bgc_S+k-1) - endif ! solve_zsal do m = 1, nbtrcr trcrn(bio_index(m) + k-1) = (trcrn(bio_index(m) + k-1) * vtmp & + ocean_bio(m)*zbgc_init_frac(m) * vsurp) / vbrin enddo enddo - if (solve_zsal) then - if (aicen > c0) then - hinS_new = vbrin/aicen - hin = vicen/aicen - else - hinS_new = c0 - hin = c0 - endif !aicen - temp_S = min_salin ! bio to cice - call remap_zbgc(nilyr, & - nt_sice, & - trtmp0(1:ntrcr), trtmp, & - 1, nblyr, & - hin, hinS_new, & - cgrid(2:nilyr+1), & - bgrid(2:nblyr+1),temp_S ) - if (icepack_warnings_aborted(subname)) return - do k = 1, nilyr - trcrn(nt_sice+k-1) = trtmp(nt_sice+k-1) - enddo !k - endif ! solve_zsal - endif ! location end subroutine adjust_tracer_profile @@ -598,24 +510,7 @@ subroutine icepack_init_bgc(ncat, nblyr, nilyr, ntrcr_o, & else ! not skl_bgc - if (scale_bgc .and. solve_zsal) then ! bulk concentration (mmol or mg per m^3) - do n = 1,ncat - do mm = 1,nbtrcr - do k = 2, nblyr - trcrn(bio_index(mm)+k-1-ntrcr_o,n) = & - (p5*(trcrn(nt_bgc_S+k-1-ntrcr_o,n)+ trcrn(nt_bgc_S+k-2-ntrcr_o,n)) & - / sss*ocean_bio_all(bio_index_o(mm))) - enddo !k - trcrn(nt_zbgc_frac-1+mm-ntrcr_o,n) = zbgc_frac_init(mm) - trcrn(bio_index(mm)-ntrcr_o,n) = (trcrn(nt_bgc_S-ntrcr_o,n) & - / sss*ocean_bio_all(bio_index_o(mm))) - trcrn(bio_index(mm)+nblyr-ntrcr_o,n) = (trcrn(nt_bgc_S+nblyr-1-ntrcr_o,n) & - / sss*ocean_bio_all(bio_index_o(mm))) - trcrn(bio_index(mm)+nblyr+1-ntrcr_o:bio_index(mm)+nblyr+2-ntrcr_o,n) = c0 ! snow - enddo ! mm - enddo ! n - - elseif (scale_bgc .and. ktherm == 2) then + if (scale_bgc .and. ktherm == 2) then trtmp(:) = c0 do n = 1,ncat call remap_zbgc(nilyr, & @@ -868,14 +763,18 @@ subroutine icepack_biogeochemistry(dt, & grow_net , & ! Specific growth rate (/s) per grid cell PP_net , & ! Total production (mg C/m^2/s) per grid cell hbri , & ! brine height, area-averaged for comparison with hi (m) - zsal_tot , & ! Total ice salinity in per grid cell (g/m^2) - fzsal , & ! Total flux of salt to ocean at time step for conservation - fzsal_g , & ! Total gravity drainage flux upNO , & ! nitrate uptake rate (mmol/m^2/d) times aice - upNH ! ammonium uptake rate (mmol/m^2/d) times aice + upNH ! ammonium uptake rate (mmol/m^2/d) times aice + + real (kind=dbl_kind), intent(inout), optional :: & + zsal_tot ! Total ice salinity in per grid cell (g/m^2) (deprecated) - logical (kind=log_kind), intent(inout) :: & - Rayleigh_criteria ! .true. means Ra_c was reached + real (kind=dbl_kind), intent(inout), optional :: & + fzsal , & ! Total flux of salt to ocean at time step for conservation (deprecated) + fzsal_g ! Total gravity drainage flux (deprecated) + + logical (kind=log_kind), intent(inout), optional :: & + Rayleigh_criteria ! .true. means Ra_c was reached (deprecated) real (kind=dbl_kind), dimension (:,:), intent(in) :: & fswpenln ! visible SW entering ice layers (W m-2) @@ -963,8 +862,6 @@ subroutine icepack_biogeochemistry(dt, & trcrn(nt_zbgc_frac-1+mm,n) = zbgc_frac_init(mm) enddo endif - if (n == 1) Rayleigh_criteria = .false. - if (solve_zsal) trcrn(nt_bgc_S:nt_bgc_S+nblyr-1,n) = c0 endif if (aicen(n) > puny) then @@ -1002,38 +899,22 @@ subroutine icepack_biogeochemistry(dt, & hsn, first_ice(n) ) if (icepack_warnings_aborted(subname)) return - if (solve_zsal) then - - call compute_microS (n, nilyr, nblyr, & - bgrid, cgrid, igrid, & - trcrn(1:ntrcr,n), hin_old(n), hbr_old, & - sss, sst, bTiz(:,n), & - iTin, bphi(:,n), kavg, & - bphi_o, Rayleigh_criteria, & - first_ice(n), bSin, brine_sal, & - brine_rho, iphin, ibrine_rho, & - ibrine_sal, sice_rho(n), sloss) - if (icepack_warnings_aborted(subname)) return - else - - ! Requires the average ice permeability = kavg(:) - ! and the surface ice porosity = zphi_o(:) - ! computed in "compute_microS" or from "thermosaline_vertical" - - iDi(:,n) = c0 - - call compute_microS_mushy (nilyr, nblyr, & - bgrid, cgrid, igrid, & - trcrn(:,n), hin_old(n), hbr_old, & - sss, sst, bTiz(:,n), & - iTin(:), bphi(:,n), kavg, & - bphi_o, bSin(:), & - brine_sal(:), brine_rho(:), iphin(:), & - ibrine_rho(:), ibrine_sal(:), sice_rho(n), & - iDi(:,n) ) - if (icepack_warnings_aborted(subname)) return - - endif ! solve_zsal + ! Requires the average ice permeability = kavg(:) + ! and the surface ice porosity = zphi_o(:) + ! computed in "compute_microS" or from "thermosaline_vertical" + + iDi(:,n) = c0 + + call compute_microS_mushy (nilyr, nblyr, & + bgrid, cgrid, igrid, & + trcrn(:,n), hin_old(n), hbr_old, & + sss, sst, bTiz(:,n), & + iTin(:), bphi(:,n), kavg, & + bphi_o, bSin(:), & + brine_sal(:), brine_rho(:), iphin(:), & + ibrine_rho(:), ibrine_sal(:), sice_rho(n), & + iDi(:,n) ) + if (icepack_warnings_aborted(subname)) return call update_hbrine (melttn(n), & meltsn (n), dt, & @@ -1052,36 +933,6 @@ subroutine icepack_biogeochemistry(dt, & hbri = hbri + hbrin * aicen(n) -#ifdef UNDEPRECATE_ZSAL - if (solve_zsal) then - - call zsalinity (n, dt, & - nilyr, bgrid, & - cgrid, igrid, & - trcrn(nt_bgc_S:nt_bgc_S+nblyr-1,n), & - trcrn(nt_qice:nt_qice+nilyr-1,n), & - trcrn(nt_sice:nt_sice+nilyr-1,n), & - ntrcr, trcrn(nt_fbri,n), & - bSin, bTiz(:,n), & - bphi(:,n), iphin, & - iki(:,n), hbr_old, & - hbrin, hin, & - hin_old(n), iDi(:,n), & - darcy_V(n), brine_sal, & - brine_rho, ibrine_sal, & - ibrine_rho, dh_direct, & - Rayleigh_criteria, & - first_ice(n), sss, & - sst, dhbr_top(n), & - dhbr_bot(n), & - fzsal, fzsal_g, & - bphi_o, nblyr, & - vicen(n), aicen_init(n), & - zsal_tot) - if (icepack_warnings_aborted(subname)) return - - endif ! solve_zsal -#endif endif ! tr_brine !----------------------------------------------------------------- diff --git a/columnphysics/icepack_zbgc_shared.F90 b/columnphysics/icepack_zbgc_shared.F90 index 51a8a8bda..757044a89 100644 --- a/columnphysics/icepack_zbgc_shared.F90 +++ b/columnphysics/icepack_zbgc_shared.F90 @@ -158,13 +158,6 @@ module icepack_zbgc_shared real (kind=dbl_kind), parameter, public :: & bphimin = 0.03_dbl_kind ! minimum porosity for zbgc only -#ifdef UNDEPRECATE_ZSAL -! these parameters are used more generally for zbgc -!----------------------------------------------------------------------- -! Parameters for zsalinity -!----------------------------------------------------------------------- -#endif - real (kind=dbl_kind), parameter, public :: & viscos_dynamic = 2.2_dbl_kind , & ! 1.8e-3_dbl_kind (pure water at 0^oC) (kg/m/s) Dm = 1.0e-9_dbl_kind, & ! molecular diffusion (m^2/s) diff --git a/columnphysics/icepack_zsalinity.F90 b/columnphysics/icepack_zsalinity.F90 deleted file mode 100644 index 8d52b022f..000000000 --- a/columnphysics/icepack_zsalinity.F90 +++ /dev/null @@ -1,1158 +0,0 @@ -!======================================================================= -! -! Vertical salinity (trcrn(nt_bgc_S)) is solved on the bio grid (bgrid and igrid) -! with domain defined by the dynamic brine height (trcrn(nt_fbri) * vicen/aicen). -! The CICE Bitz and Lipscomb thermodynamics is solved on the cgrid with height -! vicen/aicen. -! Gravity drainage is parameterized as nonlinear advection -! Flushing is incorporated in the boundary changes and a darcy flow. -! (see Jeffery et al., JGR, 2011). -! -! authors: Nicole Jeffery, LANL -! Elizabeth C. Hunke, LANL -! - module icepack_zsalinity - -! zsalinity is being deprecated -! Additional code cleanup will follow including removal of -! solve_zsal, tr_bgc_S, grid_oS, l_skS, iki, TRZS, -! Rayleigh_real, Rayleigh_criteria, and all *zsal* -! icepack_brine.F90 subroutines icepack_init_zsalinity, compute_microS -#ifdef UNDEPRECATE_ZSAL - use icepack_kinds - use icepack_parameters, only: c0, c1, c2, p001, p5, puny, rhow, depressT, gravit - use icepack_parameters, only: rhosi, min_salin, salt_loss - use icepack_parameters, only: l_skS, grid_oS, l_sk - use icepack_parameters, only: dts_b - use icepack_tracers, only: nt_sice - use icepack_zbgc_shared, only: remap_zbgc - use icepack_zbgc_shared, only: Ra_c, k_o, viscos_dynamic, thinS, Dm, exp_h - use icepack_warnings, only: warnstr, icepack_warnings_add - use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted - - use icepack_brine, only: calculate_drho - use icepack_therm_shared, only: calculate_Tin_from_qin - - implicit none - - private - public :: zsalinity - - real (kind=dbl_kind), parameter :: & - max_salin = 200.0_dbl_kind, & !(ppt) maximum bulk salinity - lapidus_g = 0.3_dbl_kind ! constant for artificial - ! viscosity/diffusion during growth -! lapidus_m = 0.007_dbl_kind ! constant for artificial diffusion during melt - -!======================================================================= - - contains - -!======================================================================= - - subroutine zsalinity (n_cat, dt, & - nilyr, bgrid, & - cgrid, igrid, & - trcrn_S, trcrn_q, & - trcrn_Si, ntrcr, & - fbri, & - bSin, bTin, & - bphin, iphin, & - ikin, hbr_old, & - hbrin, hin, & - hin_old, iDin, & - darcy_V, brine_sal, & - brine_rho, ibrine_sal, & - ibrine_rho, dh_direct, & - Rayleigh_criteria, & - first_ice, sss, & - sst, dh_top, & - dh_bot, & - fzsal, & - fzsal_g, bphi_min, & - nblyr, vicen, & - aicen, zsal_tot) - - integer (kind=int_kind), intent(in) :: & - nilyr , & ! number of ice layers - nblyr , & ! number of bio layers - ntrcr , & ! number of tracers - n_cat ! category number - - real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & - bgrid ! biology nondimensional vertical grid points - - real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & - igrid ! biology vertical interface points - - real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: & - cgrid ! CICE vertical coordinate - - real (kind=dbl_kind), intent(in) :: & - sss , & ! ocean salinity (ppt) - sst , & ! ocean temperature (oC) - hin_old , & ! old ice thickness (m) - dh_top , & ! brine change in top and bottom for diagnostics (m) - dh_bot , & ! minimum porosity - darcy_V , & ! darcy velocity (m/s) - dt , & ! time step - fbri , & ! ratio of brine height to ice thickness - hbr_old , & ! old brine height (m) - hin , & ! new ice thickness (m) - hbrin , & ! new brine height (m) - vicen , & ! ice volume (m) - aicen , & ! ice area (m) - bphi_min , & ! - dh_direct ! flooded or runoff amount (m) - - real (kind=dbl_kind), intent(inout) :: & - zsal_tot , & ! tot salinity (psu*rhosi*total vol ice) - fzsal , & ! total flux of salt out of ice over timestep(kg/m^2/s) - fzsal_g ! gravity drainage flux of salt over timestep(kg/m^2/s) - - real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: & - bTin , & ! Ice Temperature ^oC (on bio grid) - bphin ! Ice porosity (on bio grid) - - real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: & - bSin , & ! Ice salinity ppt (on bio grid) - brine_sal , & ! brine salinity (ppt) - brine_rho ! brine density (kg/m^3) - - real (kind=dbl_kind), dimension (nblyr), intent(inout) :: & - trcrn_S ! salinity tracer ppt (on bio grid) - - real (kind=dbl_kind), dimension (nilyr), intent(inout) :: & - trcrn_q , & ! enthalpy tracer - trcrn_Si ! salinity on CICE grid - - logical (kind=log_kind), intent(inout) :: & - Rayleigh_criteria ! .true. if minimun ice thickness (Ra_c) was reached - - logical (kind=log_kind), intent(in) :: & - first_ice ! for first category ice only .true. - !initialized values should be used - - real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: & - iDin , & ! Diffusivity on the igrid (1/s) - ikin ! permeability on the igrid - - real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: & - iphin , & ! porosity on the igrid - ibrine_rho , & ! brine rho on interface - ibrine_sal ! brine sal on interface - - ! local variables - - real (kind=dbl_kind) :: & - fzsaln , & ! category flux of salt out of ice over timestep(kg/m^2/s) - fzsaln_g , & ! category gravity drainage flux of salt over timestep(kg/m^2/s) - zsal_totn ! total salt content - - character(len=*),parameter :: subname='(zsalinity)' - - call solve_zsalinity (nilyr, nblyr, n_cat, dt, & - bgrid, cgrid, igrid, & - trcrn_S, trcrn_q, & - trcrn_Si, ntrcr, & - bSin, bTin, & - bphin, iphin, & - ikin, hbr_old, & - hbrin, hin, & - hin_old, iDin, & - darcy_V, brine_sal, & - brine_rho, ibrine_sal, & - ibrine_rho, dh_direct, & - Rayleigh_criteria, & - first_ice, sss, & - sst, dh_top, & - dh_bot, & - fzsaln, & - fzsaln_g, bphi_min) - if (icepack_warnings_aborted(subname)) return - - zsal_totn = c0 - - call column_sum_zsal (zsal_totn, nblyr, & - vicen, trcrn_S, & - fbri) - if (icepack_warnings_aborted(subname)) return - - call merge_zsal_fluxes (aicen, & - zsal_totn, zsal_tot, & - fzsal, fzsaln, & - fzsal_g, fzsaln_g) - if (icepack_warnings_aborted(subname)) return - - end subroutine zsalinity - -!======================================================================= -! -! update vertical salinity -! - subroutine solve_zsalinity (nilyr, nblyr, & - n_cat, dt, & - bgrid, cgrid, igrid, & - trcrn_S, trcrn_q, & - trcrn_Si, ntrcr, & - bSin, bTin, & - bphin, iphin, & - ikin, hbr_old, & - hbrin, hin, & - hin_old, iDin, & - darcy_V, brine_sal, & - brine_rho, ibrine_sal, & - ibrine_rho, dh_direct, & - Rayleigh_criteria, & - first_ice, sss, & - sst, dh_top, & - dh_bot, & - fzsaln, & - fzsaln_g, bphi_min) - - integer (kind=int_kind), intent(in) :: & - nilyr, & ! number of ice layers - nblyr, & ! number of bio layers - ntrcr, & ! number of tracers - n_cat ! category number - - real (kind=dbl_kind), intent(in) :: & - dt ! time step - - real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & - bgrid ! biology nondimensional vertical grid points - - real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & - igrid ! biology vertical interface points - - real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: & - cgrid ! CICE vertical coordinate - - real (kind=dbl_kind), intent(in) :: & - sss , & ! ocean salinity (ppt) - sst , & ! ocean temperature (oC) - hin_old , & ! old ice thickness (m) - dh_top , & ! brine change in top and bottom for diagnostics (m) - dh_bot , & - darcy_V - - real (kind=dbl_kind), intent(in) :: & - hbr_old , & ! old brine height (m) - hin , & ! new ice thickness (m) - hbrin , & ! new brine height (m) - bphi_min , & ! - dh_direct ! flooded or runoff amount (m) - - real (kind=dbl_kind), intent(out) :: & - fzsaln , & ! total flux of salt out of ice over timestep(kg/m^2/s) - fzsaln_g ! gravity drainage flux of salt over timestep(kg/m^2/s) - - real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: & - bTin , & ! Ice Temperature ^oC (on bio grid) - bphin ! Ice porosity (on bio grid) - - real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: & - bSin , & ! Ice salinity ppt (on bio grid) - brine_sal , & ! brine salinity (ppt) - brine_rho ! brine density (kg/m^3) - - real (kind=dbl_kind), dimension (nblyr), intent(inout) :: & - trcrn_S ! salinity tracer ppt (on bio grid) - - real (kind=dbl_kind), dimension (nilyr), intent(inout) :: & - trcrn_q , & ! enthalpy tracer - trcrn_Si ! salinity on CICE grid - - logical (kind=log_kind), intent(inout) :: & - Rayleigh_criteria ! .true. if minimun ice thickness (Ra_c) was reached - - logical (kind=log_kind), intent(in) :: & - first_ice ! for first category ice only .true. - !initialized values should be used - - real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: & - iDin , & ! Diffusivity on the igrid (1/s) - ikin ! permeability on the igrid - - real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: & - iphin , & ! porosity on the igrid - ibrine_rho , & ! brine rho on interface - ibrine_sal ! brine sal on interface - - ! local variables - - integer (kind=int_kind) :: & - k, nint ! vertical biology layer index - - real (kind=dbl_kind) :: & - surface_S ! salinity of ice above hin > hbrin - - real (kind=dbl_kind), dimension(2) :: & - S_bot - - real (kind=dbl_kind) :: & - Tmlts , & ! melting temperature - dts ! local timestep (s) - - logical (kind=log_kind) :: & - Rayleigh - - real (kind=dbl_kind):: & - Ttemp ! initial temp profile on the CICE grid - - real (kind=dbl_kind), dimension (ntrcr+2) :: & - trtmp0 , & ! temporary, remapped tracers !need extra - trtmp ! temporary, remapped tracers ! - - logical (kind=log_kind) :: & - cflag - - character(len=*),parameter :: subname='(solve_zsalinity)' - - !----------------------------------------------------------------- - ! Initialize - !----------------------------------------------------------------- - - dts = dts_b - nint = max(1,INT(dt/dts)) - dts = dt/nint - - !---------------------------------------------------------------- - ! Update boundary conditions - !---------------------------------------------------------------- - - surface_S = min_salin - - Rayleigh = .true. - if (n_cat == 1 .AND. hbr_old < Ra_c) then - Rayleigh = Rayleigh_criteria ! only category 1 ice can be false - endif - - if (dh_bot + darcy_V*dt > c0) then - - bSin (nblyr+2) = sss - bTin (nblyr+2) = sst - brine_sal(nblyr+2) = sss - brine_rho(nblyr+2) = rhow - bphin (nblyr+2) = c1 - S_bot (1) = c0 - S_bot (2) = c1 - - ! bottom melt - else - bSin (nblyr+2) = bSin(nblyr+1) - Tmlts = -bSin(nblyr+2)* depressT - bTin (nblyr+2) = bTin(nblyr+1) - bphin(nblyr+2) = iphin(nblyr+1) - S_bot(1) = c1 - S_bot(2) = c0 - endif - - if (abs(dh_top) > puny .AND. abs(darcy_V) > puny) then - bSin(1) = max(min_salin,-(brine_rho(2)*brine_sal(2)/rhosi & - * darcy_V*dt - (dh_top + darcy_V*dt/bphi_min - dh_direct)*min_salin & - + max(c0,-dh_direct) * sss )/dh_top) - brine_sal(1) = brine_sal(2) - brine_rho(1) = brine_rho(2) - bphin(1) = bphi_min - else - bSin(1) = min_salin - endif - - !----------------------------------------------------------------- - ! Solve for S using CICE T. If solve_zsal = .true., then couple back - ! to the thermodynamics - !----------------------------------------------------------------- - - call solve_S_dt (cflag, nblyr, & - nint , dts , & - bSin , bTin , & - bphin , iphin , & - igrid , bgrid , & - ikin , & - hbr_old , hbrin , & - hin_old , & - iDin , darcy_V , & - brine_sal , Rayleigh , & - first_ice , sss , & - dt , dh_top , & - dh_bot , brine_rho , & - ibrine_sal , ibrine_rho , & - fzsaln , fzsaln_g , & - S_bot ) - if (icepack_warnings_aborted(subname)) return - - if (n_cat == 1) Rayleigh_criteria = Rayleigh - - trtmp0(:) = c0 - trtmp (:) = c0 - - do k = 1,nblyr ! back to bulk quantity - trcrn_S(k) = bSin(k+1) - trtmp0(nt_sice+k-1) = trcrn_S(k) - enddo ! k - - call remap_zbgc (nilyr, & - nt_sice, & - trtmp0(1:ntrcr), & - trtmp, & - 1, nblyr, & - hin, hbrin, & - cgrid(2:nilyr+1), & - bgrid(2:nblyr+1), & - surface_S ) - if (icepack_warnings_aborted(subname)) return - - do k = 1, nilyr - Tmlts = -trcrn_Si(k)*depressT - Ttemp = min(-(min_salin+puny)*depressT, & - calculate_Tin_from_qin(trcrn_q(k),Tmlts)) - trcrn_Si(k) = min(-Ttemp/depressT, max(min_salin, & - trtmp(nt_sice+k-1))) - Tmlts = - trcrn_Si(k)*depressT - ! if (cflag) trcrn_q(k) = calculate_qin_from_Sin(Ttemp,Tmlts) - enddo ! k - - end subroutine solve_zsalinity - -!======================================================================= -! -! solves salt continuity explicitly using -! Lax-Wendroff-type scheme (MacCormack) -! (Mendez-Nunez and Carroll, Monthly Weather Review, 1993) -! -! authors Nicole Jeffery, LANL -! - subroutine solve_S_dt (cflag, nblyr, nint, & - dts, bSin, bTin, & - bphin, iphin, igrid, & - bgrid, ikin, hbri_old, & - hbrin, hice_old, & - iDin, darcy_V, & - brine_sal, Rayleigh, & - first_ice, sss, & - dt, dht, & - dhb, brine_rho, & - ibrine_sal, ibrine_rho, & - fzsaln, fzsaln_g, & - S_bot ) - - integer (kind=int_kind), intent(in) :: & - nblyr , & ! number of bio layers - nint ! number of interations - - logical (kind=log_kind), intent(out) :: & - cflag ! thin or not - - real (kind=dbl_kind), intent(in) :: & - dt , & ! timestep (s) - dts , & ! local timestep (s) - sss , & ! sea surface salinity - dht , & ! change in the ice top (positive for melting) - dhb , & ! change in the ice bottom (positive for freezing) - hice_old , & ! old ice thickness (m) - hbri_old , & ! brine thickness (m) - hbrin , & ! new brine thickness (m) - darcy_V ! Darcy velocity due to a pressure head (m/s) or melt - - real (kind=dbl_kind), intent(out) :: & - fzsaln , & ! salt flux +ive to ocean (kg/m^2/s) - fzsaln_g ! gravity drainage salt flux +ive to ocean (kg/m^2/s) - - logical (kind=log_kind), intent(inout) :: & - Rayleigh ! if .true. convection is allowed; if .false. not yet - - logical (kind=log_kind), intent(in) :: & - first_ice - - real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & - brine_sal , & ! Internal brine salinity (ppt) - brine_rho , & ! Internal brine density (kg/m^3) - bgrid , & ! biology nondimensional grid layer points - bTin ! Temperature of ice layers on bio grid for history (C) - - real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: & - bphin , & ! Porosity of layers - bSin ! Bulk Salinity (ppt) contains previous timestep - ! and ocean ss - - real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & - ibrine_rho , & ! brine rho on interface - ibrine_sal , & ! brine sal on interface - igrid ! biology grid interface points - - real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: & - iphin ! Porosity of layers on interface - - real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: & - iDin , & ! Diffusivity on the igrid (1/s) with minimum bphi condition - ikin ! permeability on interface - - real (kind=dbl_kind), dimension (2), intent(in) :: & - S_bot - - ! local variables - - integer (kind=int_kind) :: & - k, m ! vertical biology layer index - - real (kind=dbl_kind), dimension (nblyr+1) :: & - iDin_p , & ! Diffusivity on the igrid (1/s)/bphi^3 - dSbdx , & ! gradient of brine rho on grid - drho , & ! brine difference rho_a-rho_b (kg/m^3) -! Ci_s , & ! - Ui_s , & ! interface function -! Vi_s , & ! for conservation check - ivel - - real (kind=dbl_kind), dimension (nblyr+2) :: & - Din_p , & ! Diffusivity on the igrid (1/s)/bphi^3 - Sintemp , & ! initial salinity - pre_sin , & ! estimate of salinity of layers - pre_sinb , & ! estimate of salinity of layers - bgrid_temp , & ! biology nondimensional grid layer points - ! with boundary values - Q_s, C_s , & ! Functions in continuity equation - V_s, U_s, F_s - - real (kind=dbl_kind) :: & - dh , & ! (m) change in hbrine over dts - dbgrid , & ! ratio of grid space to spacing across boundary - ! i.e. 1/nilyr/(dbgrid(2)-dbgrid(1)) - lapidus , & ! artificial viscosity: use lapidus_g for growth - Ssum_old,Ssum_new, & ! depth integrated salt before and after timestep - fluxcorr, & ! flux correction to prevent S < min_salin - Ssum_corr, & ! numerical boundary flux correction - fluxb , & ! bottom, top and total boundary flux (g/kg/m^2) - fluxg , & ! bottom, top and total gravity drainage flux (g/kg/m^2) - fluxm , & ! bottom, top and total molecular diffusion flux (g/kg/m^2) - sum_old,sum_new , & ! integrated salinity at t and t+dt - dh_dt, dS_dt , & - Ssum_tmp - - real (kind=dbl_kind), dimension (nblyr) :: & - vel , & ! advective velocity times dt (m) - lapidus_diff , & ! lapidus term and - flux_corr , & - lapA , & - lapB - - logical (kind=log_kind) :: & - test_conservation ! test that salt change is balanced by fluxes - - character(len=*),parameter :: subname='(solve_S_dt)' - - !----------------------------------------------------------------- - ! Initialize - !----------------------------------------------------------------- - - cflag = .false. - test_conservation = .false. - iDin_p(:) = c0 - Din_p(:) = c0 - lapA(:) = c1 - lapB(:) = c1 - lapA(nblyr) = c0 - lapB(1) = c0 - V_s(:) = c0 - U_s(:) = c0 - Q_s(:) = c0 - C_s(:) = c0 -! Ci_s(:) = c0 - F_s(:) = c0 - Ui_s(:) = c0 -! Vi_s(:) = c0 - ivel(:) = c0 - vel(:) = c0 - dh = c0 - dbgrid = c2 - fzsaln = c0 - fzsaln_g = c0 - - !----------------------------------------------------------------- - ! Find brine density gradient for gravity drainage parameterization - !----------------------------------------------------------------- - - call calculate_drho(nblyr, igrid, bgrid,& - brine_rho, ibrine_rho, drho) - if (icepack_warnings_aborted(subname)) return - - !----------------------------------------------------------------- - ! Calculate bphi diffusivity on the grid points - ! rhosi = 919-974 kg/m^2 set in bio_in - ! rhow = 1026.0 density of sea water: uses kinematic viscosity (m^2/s) in Q18 - ! dynamic viscosity divided by density = kinematic. - !----------------------------------------------------------------- - - do k = 2, nblyr+1 - iDin_p(k) = k_o*gravit*l_skS/viscos_dynamic*drho(k)/(hbri_old**2) - Din_p(k) = (iDin_p(k)*(igrid(k)-bgrid(k)) & - + iDin_p(k-1)*(bgrid(k)-igrid(k-1)))/(igrid(k)-igrid(k-1)) - enddo !k - - !----------------------------------------------------------------- - ! Critical Ra_c value is only for the onset of convection in thinS - ! ice and not throughout, therefore I need a flag to indicate the - ! Ra_c was reached for a particular ice block. - ! Using a thickness minimum (Ra_c) for simplicity. - !----------------------------------------------------------------- - - bgrid_temp(:) = bgrid(:) - Din_p(nblyr+2) = iDin_p(nblyr+1) - if (.NOT. Rayleigh .AND. hbrin < Ra_c) then - Din_p(:) = c0 - iDin_p(:) = c0 - else - Rayleigh = .true. - endif - - if (hbri_old > thinS .AND. hbrin > thinS .and. & - hice_old > thinS .AND. .NOT. first_ice) then - - cflag = .true. - - bgrid_temp(1) = c0 - bgrid_temp(nblyr+2) = c1 - dbgrid = igrid(2)/(bgrid_temp(2)-bgrid_temp(1)) - - !----------------------------------- - ! surface boundary terms - !----------------------------------- - - lapidus = lapidus_g/real(nblyr,kind=dbl_kind)**2 - ivel(1) = dht/hbri_old - U_s (1) = ivel(1)/dt*dts - Ui_s(1) = U_s(1) -! Ci_s(1) = c0 - F_s (1) = brine_rho(2)*brine_sal(2)/rhosi*darcy_V*dts/hbri_old/bSin(1) - - !----------------------------------- - ! bottom boundary terms - !----------------------------------- - - ivel(nblyr+1) = dhb/hbri_old - Ui_s(nblyr+1) = ivel(nblyr+1)/dt*dts - dSbdx(nblyr) = (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) & - - ibrine_sal(nblyr)*ibrine_rho(nblyr)) & - / (igrid(nblyr+1)-igrid(nblyr)) - C_s(nblyr+1) = Dm/brine_sal(nblyr+1)/brine_rho(nblyr+1)*dts/hbri_old**2 & - * (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) & - - ibrine_sal(nblyr)*ibrine_rho(nblyr)) & - / (igrid(nblyr+1)-igrid(nblyr)) - F_s(nblyr+1) = darcy_V*dts/hbri_old/bphin(nblyr+1) - F_s(nblyr+2) = darcy_V*dts/hbri_old/bphin(nblyr+2) - vel(nblyr) =(bgrid(nblyr+1)*(dhb) -(bgrid(nblyr+1) - c1)*dht)/hbri_old - U_s(nblyr+1) = vel(nblyr)/dt*dts - V_s(nblyr+1) = Din_p(nblyr+1)/rhosi & - * (rhosi/brine_sal(nblyr+1)/brine_rho(nblyr+1))**exp_h & - * dts*dSbdx(nblyr) - dSbdx(nblyr+1) = (brine_sal(nblyr+2)*brine_rho(nblyr+2) & - - brine_sal(nblyr+1)*brine_rho(nblyr+1)) & - / (bgrid(nblyr+2)-bgrid(nblyr+1)+ grid_oS/hbri_old ) - C_s( nblyr+2) = Dm/brine_sal(nblyr+2)/brine_rho(nblyr+2)*dts/hbri_old**2 & - * (brine_sal(nblyr+2)*brine_rho(nblyr+2) & - - brine_sal(nblyr+1)*brine_rho(nblyr+1)) & - / (bgrid(nblyr+2)-bgrid(nblyr+1) + grid_oS/hbri_old ) - U_s(nblyr+2) = ivel(nblyr+1)/dt*dts - V_s(nblyr+2) = Din_p(nblyr+2)/rhosi & - * (bphin(nblyr+1)/bSin(nblyr+2))**exp_h & - * dts*dSbdx(nblyr+1) -! Ci_s(nblyr+1) = C_s(nblyr+2) -! Vi_s(nblyr+1) = V_s(nblyr+2) - dh = (dhb-dht)/dt*dts - - do k = 2, nblyr - ivel(k) = (igrid(k)*dhb - (igrid(k)-c1)*dht)/hbri_old - Ui_s(k) = ivel(k)/dt*dts -! Vi_s(k) = iDin_p(k)/rhosi & -! *(rhosi/ibrine_rho(k)/ibrine_sal(k))**exp_h*dts & -! * (brine_sal(k+1)*brine_rho(k+1) & -! - brine_sal(k)*brine_rho(k)) & -! / (bgrid(k+1)-bgrid(k)) - dSbdx(k-1) = (ibrine_sal(k)*ibrine_rho(k) & - - ibrine_sal(k-1)*ibrine_rho(k-1))/(igrid(k)-igrid(k-1)) - F_s(k) = darcy_V*dts/hbri_old/bphin(k) - C_s(k) = Dm/brine_sal(k)/brine_rho(k)*dts/hbri_old**2 & - * (ibrine_sal(k)*ibrine_rho(k) & - - ibrine_sal(k-1)*ibrine_rho(k-1))/(igrid(k)-igrid(k-1)) -! Ci_s(k) = Dm/ibrine_sal(k)/ibrine_rho(k)*dts/hbri_old**2 & -! * (brine_sal(k+1)*brine_rho(k+1) & -! - brine_sal(k)*brine_rho(k))/(bgrid(k+1)-bgrid(k)) - vel(k-1) = (bgrid(k)*(dhb) - (bgrid(k) - c1)* dht)/hbri_old - U_s(k) = vel(k-1)/dt*dts - V_s(k) = Din_p(k)/rhosi & - * (rhosi/brine_sal(k)/brine_rho(k))**exp_h*dts*dSbdx(k-1) - C_s(2) = c0 - V_s(2) = c0 - enddo !k - - !----------------------------------------------------------------- - ! Solve - !----------------------------------------------------------------- - - do m = 1, nint - - Sintemp(:) = bSin(:) - pre_sin(:) = bSin(:) - pre_sinb(:) = bSin(:) - Ssum_old = bSin(nblyr+1)*(igrid(nblyr+1)-igrid(nblyr)) - - ! forward-difference - - do k = 2, nblyr - Ssum_old = Ssum_old + bSin(k)*(igrid(k)-igrid(k-1)) - - pre_sin(k) =bSin(k) + (Ui_s(k)*(bSin(k+1) - bSin(k)) + & - V_s(k+1)*bSin(k+1)**3 - V_s(k)*bSin(k)**3 + & - (C_s(k+1)+F_s(k+1))*bSin(k+1)-& - (C_s(k)+F_s(k))*bSin(k))/(bgrid_temp(k+1)-bgrid_temp(k)) - enddo !k - - pre_sin(nblyr+1) = bSin(nblyr+1) + (Ui_s(nblyr+1)*(bSin(nblyr+2) - & - bSin(nblyr+1)) + V_s(nblyr+2)*bSin(nblyr+2)**3 - & - V_s(nblyr+1)*bSin(nblyr+1)**3+ (C_s(nblyr+2)+F_s(nblyr+2))*& - bSin(nblyr+2)-(C_s(nblyr+1)+F_s(nblyr+1))*bSin(nblyr+1) )/& - (bgrid_temp(nblyr+2)- bgrid_temp(nblyr+1)) - - ! backward-difference - - pre_sinb(2) = p5*(bSin(2) + pre_sin(2) + (Ui_s(1)& - *(pre_sin(2) -pre_sin(1)) + & - V_s(2)*pre_sin(2)**3 - & - V_s(1)*pre_sin(1)**3 + (C_s(2)+F_s(2))*pre_sin(2)-& - (C_s(1)+F_s(1))*pre_sin(1) )/(bgrid_temp(2)-bgrid_temp(1)) ) - - do k = nblyr+1, 3, -1 !nblyr+1 - pre_sinb(k) =p5*(bSin(k) + pre_sin(k) + (Ui_s(k-1)& - *(pre_sin(k) - pre_sin(k-1)) + & - V_s(k)*pre_sin(k)**3 - & - V_s(k-1)*pre_sin(k-1)**3 + (C_s(k)+F_s(k))*pre_sin(k)-& - (C_s(k-1)+F_s(k-1))*pre_sin(k-1))/(bgrid_temp(k)-bgrid_temp(k-1)) ) - - Q_s(k) = V_s(k)*pre_sin(k)**2 + U_s(k) + C_s(k) + F_s(k) - enddo !k - - Q_s(2) = V_s(2)*pre_sin(2)**2 + U_s(2) + C_s(2) + F_s(2) - Q_s(1) = V_s(1)*pre_sin(2)**2 + Ui_s(1) + C_s(1)+ F_s(1) - Q_s(nblyr+2) = V_s(nblyr+2)*pre_sin(nblyr+1)**2 + & - Ui_s(nblyr+1) + C_s(nblyr+2) + F_s(nblyr+2) - - !----------------------------------------------------------------- - ! Add artificial viscosity [Lapidus,1967] [Lohner et al, 1985] - ! * more for melting ice - !----------------------------------------------------------------- - - lapidus_diff(:) = c0 - flux_corr(:) = c0 - Ssum_new = c0 - Ssum_corr = c0 - fluxcorr = c0 - fluxg = c0 - fluxb = c0 - fluxm = c0 - - do k = 2, nblyr+1 - - lapidus_diff(k-1) = lapidus/& ! lapidus/real(nblyr,kind=dbl_kind)**2/& - (igrid(k)-igrid(k-1))* & - ( lapA(k-1)*ABS(Q_s(k+1)-Q_s(k))*(abs(pre_sinb(k+1))-abs(pre_sinb(k)))/& - (bgrid_temp(k+1)-bgrid_temp(k) )**2 - & - lapB(k-1)*ABS(Q_s(k)-Q_s(k-1))*(abs(pre_sinb(k))-abs(pre_sinb(k-1)))/& - (bgrid_temp(k)-bgrid_temp(k-1))**2) - - bSin(k) = pre_sinb(k) + lapidus_diff(k-1) - - if (bSin(k) < min_salin) then - flux_corr(k-1) = min_salin - bSin(k) ! flux into the ice - bSin(k) = min_salin - elseif (bSin(k) > -bTin(k)/depressT) then - flux_corr(k-1) = bSin(k)+bTin(k)/depressT ! flux into the ice - bSin(k) = -bTin(k)/depressT - elseif (bSin(k) > max_salin) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - call icepack_warnings_add(subname//' bSin(k) > max_salin') - endif - - if (k == nblyr+1) bSin(nblyr+2) = S_bot(1)*bSin(nblyr+1) & - + S_bot(2)*bSin(nblyr+2) - - Ssum_new = Ssum_new + bSin(k)*(igrid(k)-igrid(k-1)) - fluxcorr = fluxcorr + (flux_corr(k-1) & - + lapidus_diff(k-1))*(igrid(k)-igrid(k-1)) - - enddo !k - - Ssum_tmp = Ssum_old - - call calc_salt_fluxes (nint, nblyr, & - Ui_s, dh,dbgrid,hbri_old,Sintemp, & - pre_sin, fluxb,fluxg,fluxm,V_s, & - C_s, F_s, Ssum_corr,fzsaln_g,fzsaln, & - Ssum_tmp,dts, Ssum_new) - if (icepack_warnings_aborted(subname)) return - - if (test_conservation) then - call check_conserve_salt(nint, m, dt, & - Ssum_tmp, Ssum_new, Ssum_corr,& - fluxcorr, fluxb, fluxg, fluxm, & - hbrin, hbri_old) - call icepack_warnings_add(subname//' check_conserve_salt fails') - if (icepack_warnings_aborted(subname)) return - endif ! test_conservation - - enddo !m - - else ! add/melt ice only - - sum_old = c0 - sum_new = c0 - dh_dt = hbrin-hbri_old - dS_dt = c0 - if (dh_dt > c0) then - dS_dt = sss*dh_dt*salt_loss - do k = 2, nblyr+1 - bSin(k) = max(min_salin,(bSin(k)*hbri_old + dS_dt)/hbrin) - enddo !k - else - do k = 2, nblyr+1 - sum_old = sum_old + bSin(k)*hbri_old*(igrid(k)-igrid(k-1)) - bSin(k) = max(min_salin,bSin(k) * (c1-abs(dh_dt)/hbri_old)) - sum_new = sum_new + bSin(k)*hbrin*(igrid(k)-igrid(k-1)) - enddo !k - endif - fzsaln = rhosi*(sum_old-sum_new + dS_dt)*p001/dt !kg/m^2/s - fzsaln_g = c0 - - endif ! (hbri_old > thinS .AND. hbrin > thinS & - ! .and. hice_old > thinS .AND. .NOT. first_ice) - - !----------------------------------------------------------------- - ! Move this to bgc calculation if using tr_salinity - ! Calculate bphin, iphin, ikin, iDin and iDin_N - !----------------------------------------------------------------- - - iDin(:) = c0 - iphin(:) = c1 - ikin(:) = c0 - - do k = 1, nblyr+1 - if (k < nblyr+1) bphin(k+1) = min(c1,max(puny, & - bSin(k+1)*rhosi/(brine_sal(k+1)*brine_rho(k+1)))) - if (k == 1) then - bphin(k) = min(c1,max(puny, bSin(k)*rhosi/(brine_sal(k)*brine_rho(k)))) - iphin(k) = bphin(2) - elseif (k == nblyr+1) then - iphin(nblyr+1) = bphin(nblyr+1) - else - iphin(k) = min(c1, max(c0,(bphin(k+1) - bphin(k))/(bgrid(k+1) & - - bgrid(k))*(igrid(k)-bgrid(k)) + bphin(k))) - endif - ikin(k) = k_o*iphin(k)**exp_h - enddo !k - - if (cflag) then - - do k = 2, nblyr+1 - iDin(k) = iphin(k)*Dm/hbri_old**2 - if (Rayleigh .AND. hbrin .GE. Ra_c) & - iDin(k) = iDin(k) + l_sk*ikin(k)*gravit/viscos_dynamic & - * drho(k)/hbri_old**2 - enddo !k - else ! .not. cflag - do k = 2, nblyr+1 - iDin(k) = iphin(k)*Dm/hbri_old**2 - enddo !k - endif - - end subroutine solve_S_dt - -!======================================================================= -! -! Calculate salt fluxes -! - subroutine calc_salt_fluxes (mint, nblyr, & - Ui_s,dh,dbgrid,hbri_old,Sintemp,pre_sin,& - fluxb,fluxg,fluxm,V_s,& - C_s,F_s,Ssum_corr,fzsaln_g,fzsaln,Ssum_old, & -! fluxcorr,dts, Ssum_new) - dts, Ssum_new) - - integer(kind=int_kind), intent(in) :: & - nblyr, & ! number of bio layers - mint ! current iteration - - real (kind=dbl_kind), intent(in) :: & - dts , & ! halodynamic timesteps (s) - ! hbrin , & ! new brine height after all iterations (m) - dh , & ! (m) change in hbrine over dts - dbgrid , & ! ratio of grid space to spacing across boundary - ! ie. 1/nilyr/(dbgrid(2)-dbgrid(1)) -! fluxcorr , & ! flux correction to ensure S >= min_salin - hbri_old ! initial brine height (m) - - real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & - Ui_s ! interface function - - real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: & - Sintemp , & ! initial salinity - pre_sin , & ! estimate of salinity of layers - C_s , & ! Functions in continuity equation - F_s , & - V_s - - real (kind=dbl_kind), intent(in) :: & - Ssum_old , & ! initial integrated salt content (ppt)/h - Ssum_new ! next integrated salt content(ppt)/h - - real (kind=dbl_kind), intent(inout) :: & - fzsaln , & ! total salt flux out of ice over timestep(kg/m^2/s) - fzsaln_g , & ! gravity drainage flux of salt over timestep(kg/m^2/s) - Ssum_corr, & ! boundary flux correction due to numerics - fluxb , & ! total boundary salt flux into the ice (+ into ice) - fluxg , & ! total gravity drainage salt flux into the ice (+ into ice) - fluxm ! total molecular diffusive salt flux into the ice (+ into ice) - - ! local variables - - real (kind=dbl_kind) :: & - Ssum_corr_flux , & ! numerical boundary flux correction - fluxb_b, fluxb_t, & ! bottom, top and total boundary flux (g/kg/m^2) - fluxg_b, fluxg_t, & ! bottom, top and total gravity drainage flux (g/kg/m^2) - fluxm_b, fluxm_t ! bottom, top and total molecular diffusion flux (g/kg/m^2) - - real (kind=dbl_kind) :: hin_old, hin_next, dhtmp !, dh - - character(len=*),parameter :: subname='(calc_salt_fluxes)' - - dhtmp = c1-dh/hbri_old - hin_next = hbri_old + real(mint,kind=dbl_kind)*dh - hin_old = hbri_old + (real(mint,kind=dbl_kind)-c1)*dh - - !----------------------------------------------------------------- - ! boundary fluxes (positive into the ice) - !--------------------------------------------- - ! without higher order numerics corrections - ! fluxb = (Ui_s(nblyr+1) + F_s(nblyr+2))*Sintemp(nblyr+2) - (Ui_s(1) + F_s(1))*Sintemp(1) - !----------------------------------------------------------------- - - fluxb_b = p5* Ui_s(nblyr+1) * (dhtmp*Sintemp(nblyr+2)*dbgrid & - + pre_sin(nblyr+1) & - + dhtmp*Sintemp(nblyr+1)*(c1-dbgrid)) & - + p5*(F_s(nblyr+2) * dhtmp*Sintemp(nblyr+2)*dbgrid & - + F_s(nblyr+1) * (pre_sin(nblyr+1) & - + dhtmp*Sintemp(nblyr+1)*(c1-dbgrid))) - - fluxb_t = -p5*Ui_s(1)*(pre_sin(1)*dbgrid + & - dhtmp*Sintemp(2) - & - (dbgrid-c1)*pre_sin(2)) - & - p5*(dbgrid*F_s(1)*pre_sin(1) + & - F_s(2)*(dhtmp*Sintemp(2) & - +(c1-dbgrid)*pre_sin(2))) - - fluxb = fluxb_b + fluxb_t - - !----------------------------------------------------------------- - ! gravity drainage fluxes (positive into the ice) - ! without higher order numerics corrections - ! fluxg = V_s(nblyr+2)*Sintemp(nblyr+1)**3 - !----------------------------------------------------------------- - - fluxg_b = p5*(dhtmp* dbgrid* & - V_s(nblyr+2)*Sintemp(nblyr+1)**3 + & - V_s(nblyr+1)*(pre_sin(nblyr+1)**3 - & - dhtmp*(dbgrid - c1)* & - Sintemp(nblyr+1)**3)) - - fluxg_t = -p5*(dbgrid*V_s(1)*pre_sin(1)**3 + & - V_s(2)*(dhtmp*Sintemp(2)**3- & - (dbgrid-c1)*pre_sin(2)**3)) - - fluxg = fluxg_b + fluxg_t - - !----------------------------------------------------------------- - ! diffusion fluxes (positive into the ice) - ! without higher order numerics corrections - ! fluxm = C_s(nblyr+2)*Sintemp(nblyr+2) - !----------------------------------------------------------------- - - fluxm_b = p5*(dhtmp*C_s(nblyr+2)* Sintemp(nblyr+2)*dbgrid & - + C_s(nblyr+1)*(pre_sin(nblyr+1) & - + dhtmp * Sintemp(nblyr+1)*(c1-dbgrid))) - - fluxm_t = -p5 * (C_s(1) * pre_sin(1)*dbgrid & - + C_s(2) * (pre_sin(2)*(c1-dbgrid) + dhtmp*Sintemp(2))) - - fluxm = fluxm_b + fluxm_t - - Ssum_corr = (-dh/hbri_old + p5*(dh/hbri_old)**2)*Ssum_old - Ssum_corr_flux = dh*Ssum_old/hin_next + Ssum_corr - Ssum_corr = Ssum_corr_flux - - fzsaln_g = fzsaln_g - hin_next * fluxg_b & - * rhosi*p001/dts - - !approximate fluxes - !fzsaln = fzsaln - hin_next * (fluxg & - ! + fluxb + fluxm + fluxcorr + Ssum_corr_flux) & - ! * rhosi*p001/dts - - fzsaln = fzsaln + (Ssum_old*hin_old - Ssum_new*hin_next) & - * rhosi*p001/dts ! positive into the ocean - - end subroutine calc_salt_fluxes - -!======================================================================= -! -! Test salt conservation: flux conservative form d(hSin)/dt = -dF(x,Sin)/dx -! - subroutine check_conserve_salt (mmax, mint, dt, & - Ssum_old, Ssum_new, Ssum_corr, & - fluxcorr, fluxb, fluxg, fluxm, & - hbrin, hbri_old) - - integer(kind=int_kind), intent(in) :: & - mint , & ! current iteration - mmax ! maximum number of iterations - - real (kind=dbl_kind), intent(in) :: & - dt , & ! thermodynamic and halodynamic timesteps (s) - hbrin , & ! (m) final brine height - hbri_old , & ! (m) initial brine height - Ssum_old , & ! initial integrated salt content - Ssum_new , & ! final integrated salt content - fluxcorr , & ! flux correction to ensure S >= min_salin - Ssum_corr , & ! boundary flux correction due to numerics - fluxb , & ! total boundary salt flux into the ice (+ into ice) - fluxg , & ! total gravity drainage salt flux into the ice (+ into ice) - fluxm ! - - ! local variables - - real (kind=dbl_kind):: & - diff2 , & ! - dsum_flux , & ! salt change in kg/m^2/s - flux_tot , & ! fluxg + fluxb - order , & ! - dh - - real (kind=dbl_kind), parameter :: & - accuracy = 1.0e-7_dbl_kind ! g/kg/m^2/s difference between boundary fluxes - - character(len=*),parameter :: subname='(check_conserve_salt)' - - dh = (hbrin-hbri_old)/real(mmax,kind=dbl_kind) - - flux_tot = (fluxb + fluxg + fluxm + fluxcorr + Ssum_corr)*& - (hbri_old + (real(mint,kind=dbl_kind))*dh)/dt - dsum_flux =(Ssum_new*(hbri_old + (real(mint,kind=dbl_kind))*dh) - & - Ssum_old*(hbri_old + (real(mint,kind=dbl_kind)-c1)* & - dh) )/dt - order = abs(dh/min(hbri_old,hbrin)) - if (abs(dsum_flux) > accuracy) then - diff2 = abs(dsum_flux - flux_tot) - if (diff2 > puny .AND. diff2 > order ) then - call icepack_warnings_setabort(.true.,__FILE__,__LINE__) - write(warnstr,*) subname, 'Poor salt conservation: check_conserve_salt' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'mint:', mint - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'Ssum_corr',Ssum_corr - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'fluxb,fluxg,fluxm,flux_tot,fluxcorr:' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, fluxb,fluxg,fluxm,flux_tot,fluxcorr - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'fluxg,',fluxg - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'dsum_flux,',dsum_flux - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'Ssum_new,Ssum_old,hbri_old,dh:' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, Ssum_new,Ssum_old,hbri_old,dh - call icepack_warnings_add(warnstr) - write(warnstr,*) subname, 'diff2,order,puny',diff2,order,puny - call icepack_warnings_add(warnstr) - endif - endif - - end subroutine check_conserve_salt - -!======================================================================= -! -! Aggregate flux information from all ice thickness categories -! - subroutine merge_zsal_fluxes(aicenS, & - zsal_totn, zsal_tot, & - fzsal, fzsaln, & - fzsal_g, fzsaln_g) - - ! single category fluxes - real (kind=dbl_kind), intent(in):: & - aicenS , & ! concentration of ice - fzsaln , & ! salt flux (kg/m**2/s) - fzsaln_g ! gravity drainage salt flux (kg/m**2/s) - - real (kind=dbl_kind), intent(in):: & - zsal_totn ! tot salinity in category (psu*volume*rhosi) - - real (kind=dbl_kind), intent(inout):: & - zsal_tot, & ! tot salinity (psu*rhosi*total vol ice) - fzsal , & ! salt flux (kg/m**2/s) - fzsal_g ! gravity drainage salt flux (kg/m**2/s) - - character(len=*),parameter :: subname='(merge_zsal_fluxes)' - - !----------------------------------------------------------------- - ! Merge fluxes - !----------------------------------------------------------------- - - zsal_tot = zsal_tot + zsal_totn ! already *aicenS - - ! ocean tot and gravity drainage salt fluxes - fzsal = fzsal + fzsaln * aicenS - fzsal_g = fzsal_g + fzsaln_g * aicenS - - end subroutine merge_zsal_fluxes - -!========================================================================== -! -! For each grid cell, sum field over all ice layers. "Net" refers to the column -! integration while "avg" is normalized by the ice thickness - - subroutine column_sum_zsal (zsal_totn, nblyr, & - vicenS, trcrn_S, fbri) - - integer (kind=int_kind), intent(in) :: & - nblyr ! number of layers - - real (kind=dbl_kind), intent(in):: & - vicenS , & ! volume of ice (m) - fbri ! brine height to ice thickness ratio - - real (kind=dbl_kind), dimension (nblyr), intent(in) :: & - trcrn_S ! input field - - real (kind=dbl_kind), intent(inout) :: & - zsal_totn ! avg salinity (psu*rhosi*vol of ice) - - ! local variables - - integer (kind=int_kind) :: & - k ! layer index - - character(len=*),parameter :: subname='(column_sum_zsal)' - - do k = 1, nblyr - zsal_totn = zsal_totn & - + rhosi * trcrn_S(k) & - * fbri & - * vicenS/real(nblyr,kind=dbl_kind) - enddo ! k - - end subroutine column_sum_zsal -#endif -!======================================================================= - - end module icepack_zsalinity - -!======================================================================= diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index 07c1d465d..7bba1d567 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -183,7 +183,6 @@ subroutine init_restart logical (kind=log_kind) :: & skl_bgc, & ! from icepack z_tracers, & ! from icepack - solve_zsal, & ! from icepack tr_brine, & ! from icepack tr_fsd ! from icepack @@ -195,7 +194,6 @@ subroutine init_restart call icepack_query_parameters(skl_bgc_out=skl_bgc) call icepack_query_parameters(z_tracers_out=z_tracers) - call icepack_query_parameters(solve_zsal_out=solve_zsal) call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -212,7 +210,7 @@ subroutine init_restart call calendar (time) endif - if (solve_zsal .or. skl_bgc .or. z_tracers) then + if (skl_bgc .or. z_tracers) then if (tr_fsd) then write (nu_diag,*) 'FSD implementation incomplete for use with BGC' call icedrv_system_abort(string=subname,file=__FILE__,line=__LINE__) diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90 index db97c79c6..e20ebf87b 100644 --- a/configuration/driver/icedrv_RunMod.F90 +++ b/configuration/driver/icedrv_RunMod.F90 @@ -98,7 +98,7 @@ subroutine ice_step use icedrv_calendar, only: dt, dt_dyn, ndtd, diagfreq, write_restart, istep use icedrv_diagnostics, only: runtime_diags, init_mass_diags ! use icedrv_diagnostics, only: icedrv_diagnostics_debug - use icedrv_diagnostics_bgc, only: hbrine_diags, zsal_diags, bgc_diags + use icedrv_diagnostics_bgc, only: hbrine_diags, bgc_diags use icedrv_flux, only: init_history_therm, init_history_bgc, & daidtt, daidtd, dvidtt, dvidtd, dagedtt, dagedtd, init_history_dyn use icedrv_history, only: history_format, history_write @@ -112,7 +112,7 @@ subroutine ice_step k ! dynamics supercycling index logical (kind=log_kind) :: & - calc_Tsfc, skl_bgc, solve_zsal, z_tracers, tr_brine, & ! from icepack + calc_Tsfc, skl_bgc, z_tracers, tr_brine, & ! from icepack tr_fsd, wave_spec, tr_snow real (kind=dbl_kind) :: & @@ -127,8 +127,7 @@ subroutine ice_step !----------------------------------------------------------------- call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) - call icepack_query_parameters(solve_zsal_out=solve_zsal, & - calc_Tsfc_out=calc_Tsfc, & + call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc, & wave_spec_out=wave_spec) call icepack_query_tracer_flags(tr_brine_out=tr_brine,tr_fsd_out=tr_fsd, & tr_snow_out=tr_snow) @@ -223,7 +222,6 @@ subroutine ice_step if (mod(istep,diagfreq) == 0) then call runtime_diags(dt) ! log file - if (solve_zsal) call zsal_diags if (skl_bgc .or. z_tracers) call bgc_diags if (tr_brine) call hbrine_diags endif @@ -234,7 +232,7 @@ subroutine ice_step if (write_restart == 1) then call dumpfile ! core variables for restarting - if (solve_zsal .or. skl_bgc .or. z_tracers) & + if (skl_bgc .or. z_tracers) & call write_restart_bgc ! biogeochemistry call final_restart endif @@ -250,7 +248,7 @@ end subroutine ice_step subroutine coupling_prep use icedrv_arrays_column, only: alvdfn, alidfn, alvdrn, alidrn, & - albicen, albsnon, albpndn, apeffn, fzsal_g, fzsal, snowfracn + albicen, albsnon, albpndn, apeffn, snowfracn use icedrv_calendar, only: dt use icedrv_domain_size, only: ncat, nx use icedrv_flux, only: alvdf, alidf, alvdr, alidr, albice, albsno, & @@ -260,7 +258,7 @@ subroutine coupling_prep fswthru_ai, fhocn, fswthru, scale_factor, snowfrac, & swvdr, swidr, swvdf, swidf, & frzmlt_init, frzmlt, & - fzsal_ai, fzsal_g_ai, flux_bio, flux_bio_ai + flux_bio, flux_bio_ai use icedrv_forcing, only: oceanmixed_ice use icedrv_state, only: aicen use icedrv_step, only: ocean_mixed_layer @@ -360,8 +358,6 @@ subroutine coupling_prep fsalt_ai (i) = fsalt (i) fhocn_ai (i) = fhocn (i) fswthru_ai(i) = fswthru(i) - fzsal_ai (i) = fzsal (i) - fzsal_g_ai(i) = fzsal_g(i) if (nbtrcr > 0) then do k = 1, nbtrcr diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90 index 66027b2d0..2e2a4a181 100644 --- a/configuration/driver/icedrv_arrays_column.F90 +++ b/configuration/driver/icedrv_arrays_column.F90 @@ -223,29 +223,13 @@ module icedrv_arrays_column darcy_V ! darcy velocity positive up (m/s) real (kind=dbl_kind), dimension (nx), public :: & - zsal_tot , & ! Total ice salinity in per grid cell (g/m^2) chl_net , & ! Total chla (mg chla/m^2) per grid cell NO_net ! Total nitrate per grid cell - logical (kind=log_kind), dimension (nx), public :: & - Rayleigh_criteria ! .true. means Ra_c was reached - - real (kind=dbl_kind), dimension (nx), public :: & - Rayleigh_real ! .true. = c1, .false. = c0 - real (kind=dbl_kind), & dimension (nx,ncat), public :: & sice_rho ! avg sea ice density (kg/m^3) ! ech: diagnostic only? - real (kind=dbl_kind), & - dimension (nx,ncat), public :: & - fzsaln , & ! category fzsal(kg/m^2/s) - fzsaln_g ! salt flux from gravity drainage only - - real (kind=dbl_kind), dimension (nx), public :: & - fzsal , & ! Total flux of salt to ocean at time step for conservation - fzsal_g ! Total gravity drainage flux - real (kind=dbl_kind), dimension (nx,nblyr+1,ncat), public :: & zfswin ! Shortwave flux into layers interpolated on bio grid (W/m^2) diff --git a/configuration/driver/icedrv_diagnostics_bgc.F90 b/configuration/driver/icedrv_diagnostics_bgc.F90 index 1eab55021..ab13a07d6 100644 --- a/configuration/driver/icedrv_diagnostics_bgc.F90 +++ b/configuration/driver/icedrv_diagnostics_bgc.F90 @@ -20,7 +20,7 @@ module icedrv_diagnostics_bgc implicit none private - public :: hbrine_diags, bgc_diags, zsal_diags + public :: hbrine_diags, bgc_diags !======================================================================= @@ -737,182 +737,6 @@ subroutine bgc_diags () end subroutine bgc_diags -!======================================================================= -! -! Writes diagnostic info (max, min, global sums, etc) to standard out -! -! authors: Nicole Jeffery, LANL - - subroutine zsal_diags () - - use icedrv_arrays_column, only: fzsal, fzsal_g, sice_rho, bTiz - use icedrv_arrays_column, only: iDi, bphi, dhbr_top, dhbr_bot, darcy_V - use icedrv_domain_size, only: nblyr, ncat, nilyr - use icedrv_state, only: aicen, aice, vice, trcr, trcrn, vicen, vsno - - ! local variables - - integer (kind=int_kind) :: & - k, n, nn - - ! fields at diagnostic points - real (kind=dbl_kind), dimension(nx) :: & - phinS, phinS1, & - phbrn,pdh_top1,pdh_bot1, psice_rho, pfzsal, & - pfzsal_g, pdarcy_V1 - - ! vertical fields of category 1 at diagnostic points for bgc layer model - real (kind=dbl_kind), dimension(nblyr+2) :: & - pphin, pphin1 - real (kind=dbl_kind), dimension(nilyr) :: & - pSice - real (kind=dbl_kind), dimension(nblyr) :: & - pSin, pSin1 - real (kind=dbl_kind), dimension(nblyr+1) :: & - pbTiz, piDin - - real (kind=dbl_kind) :: & - rhosi, rhow, rhos - - logical (kind=log_kind) :: tr_brine - integer (kind=int_kind) :: nt_fbri, nt_bgc_S, nt_sice - - character(len=*), parameter :: subname='(zsal_diags)' - - !----------------------------------------------------------------- - ! query Icepack values - !----------------------------------------------------------------- - - call icepack_query_parameters(rhosi_out=rhosi, rhow_out=rhow, rhos_out=rhos) - call icepack_query_tracer_flags(tr_brine_out=tr_brine) - call icepack_query_tracer_indices(nt_fbri_out=nt_fbri, & - nt_bgc_S_out=nt_bgc_S, nt_sice_out=nt_sice) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & - file=__FILE__,line= __LINE__) - - !----------------------------------------------------------------- - ! salinity and microstructure of the ice - !----------------------------------------------------------------- - ! NOTE these are computed for the last timestep only (not avg) - - do n = 1, nx - pfzsal = fzsal(n) - pfzsal_g = fzsal_g(n) - phinS = c0 - phinS1 = c0 - phbrn = c0 - psice_rho = c0 - pdh_top1 = c0 - pdh_bot1 = c0 - pdarcy_V1 = c0 - - do nn = 1, ncat - psice_rho = psice_rho(n) + sice_rho(n,nn)*aicen(n,nn) - enddo - if (aice(n) > c0) psice_rho = psice_rho/aice(n) - - if (tr_brine .and. aice(n) > c0) then - phinS = trcr(n,nt_fbri)*vice(n)/aice(n) - phbrn = (c1 - rhosi/rhow)*vice(n)/aice(n) & - - rhos /rhow *vsno(n)/aice(n) - endif - if (tr_brine .and. aicen(n,1) > c0) then - phinS1 = trcrn(n,nt_fbri,1)*vicen(n,1)/aicen(n,1) - pdh_top1 = dhbr_top(n,1) - pdh_bot1 = dhbr_bot(n,1) - pdarcy_V1 = darcy_V(n,1) - endif - - do k = 1, nblyr+1 - pbTiz(k) = c0 - piDin(k) = c0 - do nn = 1,ncat - pbTiz(k) = pbTiz(k) + bTiz(n,k,nn)*vicen(n,nn) - piDin(k) = piDin(k) + iDi(n,k,nn)*vicen(n,nn) - enddo - if (vice(n) > c0) then - pbTiz(k) = pbTiz(k)/vice(n) - piDin(k) = piDin(k)/vice(n) - endif - enddo ! k - do k = 1, nblyr+2 - pphin (k) = c0 - pphin1(k) = c0 - if (aicen(n,1) > c0) pphin1(k) = bphi(n,k,1) - do nn = 1,ncat - pphin(k) = pphin(k) + bphi(n,k,nn)*vicen(n,nn) - enddo - if (vice(n) > c0) pphin(k) = pphin(k)/vice(n) - enddo - do k = 1,nblyr - pSin (k) = c0 - pSin1(k) = c0 - pSin (k) = trcr(n,nt_bgc_S+k-1) - if (aicen(n,1) > c0) pSin1(k) = trcrn(n,nt_bgc_S+k-1,1) - enddo - - do k = 1,nilyr - pSice(k) = trcr(n,nt_sice+k-1) - enddo - - !----------------------------------------------------------------- - ! start spewing - !----------------------------------------------------------------- - - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,*) ' Brine height ' - write(nu_diag_out+n-1,900) 'hbrin = ',phinS(n) - write(nu_diag_out+n-1,900) 'hbrin cat 1 = ',phinS1(n) - write(nu_diag_out+n-1,900) 'Freeboard = ',phbrn(n) - write(nu_diag_out+n-1,900) 'dhbrin cat 1 top = ',pdh_top1(n) - write(nu_diag_out+n-1,900) 'dhbrin cat 1 bottom = ',pdh_bot1(n) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,*) ' zSalinity ' - write(nu_diag_out+n-1,900) 'Avg density (kg/m^3) = ',psice_rho(n) - write(nu_diag_out+n-1,900) 'Salt flux (kg/m^2/s) = ',pfzsal(n) - write(nu_diag_out+n-1,900) 'Grav. Drain. Salt flux = ',pfzsal_g(n) - write(nu_diag_out+n-1,900) 'Darcy V cat 1 (m/s) = ',pdarcy_V1(n) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,*) ' Top down bgc Layer Model' - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,803) 'bTiz ice temp' - write(nu_diag_out+n-1,*) '---------------------------------------------------' - write(nu_diag_out+n-1,802) (pbTiz(k), k = 1,nblyr+1) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,803) 'iDi diffusivity' - write(nu_diag_out+n-1,*) '---------------------------------------------------' - write(nu_diag_out+n-1,802) (piDin(k), k = 1,nblyr+1) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,803) 'bphi porosity' - write(nu_diag_out+n-1,*) '---------------------------------------------------' - write(nu_diag_out+n-1,802) (pphin(k), k = 1,nblyr+1) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,803) 'phi1 porosity' - write(nu_diag_out+n-1,*) '---------------------------------------------------' - write(nu_diag_out+n-1,802) (pphin1(k), k = 1,nblyr+1) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,803) 'zsal cat 1' - write(nu_diag_out+n-1,*) '---------------------------------------------------' - write(nu_diag_out+n-1,802) (pSin1(k), k = 1,nblyr) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,803) 'zsal Avg S' - write(nu_diag_out+n-1,*) '---------------------------------------------------' - write(nu_diag_out+n-1,802) (pSin(k), k = 1,nblyr) - write(nu_diag_out+n-1,*) ' ' - write(nu_diag_out+n-1,803) 'Sice Ice S' - write(nu_diag_out+n-1,*) '---------------------------------------------------' - write(nu_diag_out+n-1,802) (pSice(k), k = 1,nilyr) - write(nu_diag_out+n-1,*) ' ' - - enddo ! nx - -802 format (f24.17,2x) -803 format (a25,2x) -900 format (a25,2x,f24.17) - - end subroutine zsal_diags - !======================================================================= end module icedrv_diagnostics_bgc diff --git a/configuration/driver/icedrv_domain_size.F90 b/configuration/driver/icedrv_domain_size.F90 index b556f3ea1..75bb4d28e 100644 --- a/configuration/driver/icedrv_domain_size.F90 +++ b/configuration/driver/icedrv_domain_size.F90 @@ -34,8 +34,7 @@ module icedrv_domain_size ! *** add to kscavz in icepack_zbgc_shared.F90 n_bgc = (n_algae*2 + n_doc + n_dic + n_don + n_fed + n_fep + n_zaero & + 8) , & ! nit, am, sil, dmspp, dmspd, dms, pon, humic - nltrcr = (n_bgc*TRBGCZ+TRZS)*TRBRI, & ! number of zbgc (includes zaero) - ! and zsalinity tracers + nltrcr = (n_bgc*TRBGCZ)*TRBRI, & ! number of zbgc (includes zaero) max_nsw = (nilyr+nslyr+2) & ! total chlorophyll plus aerosols * (1+TRZAERO),& ! number of tracers active in shortwave calculation max_ntrcr = 1 & ! 1 = surface temperature @@ -53,7 +52,6 @@ module icedrv_domain_size + n_aero*4 & ! number of aerosols * 4 aero layers + TRBRI & ! brine height + TRBGCS*n_bgc & ! skeletal layer BGC - + TRZS *TRBRI* nblyr & ! zsalinity (off if TRBRI=0) + n_bgc*TRBGCZ*TRBRI*(nblyr+3) & ! zbgc (off if TRBRI=0) + n_bgc*TRBGCZ & ! mobile/stationary phase tracer + 1 ! for unused tracer flags diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 749e52396..06b0d06d0 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -324,10 +324,6 @@ module icedrv_flux flux_bio , & ! all bio fluxes to ocean flux_bio_ai ! all bio fluxes to ocean, averaged over grid cell - real (kind=dbl_kind), dimension (nx), public :: & - fzsal_ai, & ! salt flux to ocean from zsalinity (kg/m^2/s) - fzsal_g_ai ! gravity drainage salt flux to ocean (kg/m^2/s) - ! internal logical (kind=log_kind), public :: & @@ -782,7 +778,7 @@ subroutine init_history_bgc use icedrv_arrays_column, only: PP_net, grow_net, hbri use icedrv_arrays_column, only: ice_bio_net, snow_bio_net, fbio_snoice, fbio_atmice - use icedrv_arrays_column, only: fzsal, fzsal_g, zfswin + use icedrv_arrays_column, only: zfswin character(len=*), parameter :: subname='(init_history_bgc)' PP_net (:) = c0 @@ -794,8 +790,6 @@ subroutine init_history_bgc snow_bio_net(:,:) = c0 fbio_snoice (:,:) = c0 fbio_atmice (:,:) = c0 - fzsal (:) = c0 - fzsal_g (:) = c0 zfswin (:,:,:) = c0 fnit (:) = c0 fsil (:) = c0 diff --git a/configuration/driver/icedrv_init_column.F90 b/configuration/driver/icedrv_init_column.F90 index c32c524ae..ae21c5871 100644 --- a/configuration/driver/icedrv_init_column.F90 +++ b/configuration/driver/icedrv_init_column.F90 @@ -26,7 +26,7 @@ module icedrv_init_column use icepack_intfc, only: icepack_init_zbgc use icepack_intfc, only: icepack_init_thermo use icepack_intfc, only: icepack_step_radiation, icepack_init_orbit - use icepack_intfc, only: icepack_init_bgc, icepack_init_zsalinity + use icepack_intfc, only: icepack_init_bgc use icepack_intfc, only: icepack_init_ocean_bio, icepack_load_ocean_bio_array use icepack_intfc, only: icepack_init_hbrine use icedrv_system, only: icedrv_system_abort @@ -371,7 +371,6 @@ subroutine init_bgc() use icedrv_arrays_column, only: zfswin, trcrn_sw use icedrv_arrays_column, only: ocean_bio_all, ice_bio_net, snow_bio_net use icedrv_arrays_column, only: cgrid, igrid, bphi, iDi, bTiz, iki - use icedrv_arrays_column, only: Rayleigh_criteria, Rayleigh_real use icedrv_calendar, only: istep1 use icedrv_system, only: icedrv_system_abort use icedrv_flux, only: sss, nit, amm, sil, dmsp, dms, algalN, & @@ -389,13 +388,6 @@ subroutine init_bgc() integer (kind=int_kind) :: & max_nbtrcr, max_algae, max_don, max_doc, max_dic, max_aero, max_fe - logical (kind=log_kind) :: & - RayleighC, & - solve_zsal - - real(kind=dbl_kind) :: & - RayleighR - real(kind=dbl_kind), dimension(max_ntrcr,ncat) :: & trcrn_bgc @@ -403,8 +395,7 @@ subroutine init_bgc() sicen integer (kind=int_kind) :: & - nbtrcr, ntrcr, ntrcr_o, & - nt_sice, nt_bgc_S + nbtrcr, ntrcr, ntrcr_o, nt_sice character(len=*), parameter :: subname='(init_bgc)' @@ -426,43 +417,16 @@ subroutine init_bgc() ! query Icepack values !----------------------------------------------------------------- - call icepack_query_parameters(solve_zsal_out=solve_zsal) call icepack_query_tracer_sizes(max_nbtrcr_out=max_nbtrcr, & max_algae_out=max_algae, max_don_out=max_don, max_doc_out=max_doc, & max_dic_out=max_dic, max_aero_out=max_aero, max_fe_out=max_fe) call icepack_query_tracer_sizes(nbtrcr_out=nbtrcr, ntrcr_out=ntrcr, & ntrcr_o_out=ntrcr_o) - call icepack_query_tracer_indices(nt_sice_out=nt_sice, nt_bgc_S_out=nt_bgc_S) + call icepack_query_tracer_indices(nt_sice_out=nt_sice) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - !----------------------------------------------------------------- - ! zsalinity initialization - !----------------------------------------------------------------- - - if (solve_zsal) then - do i = 1, nx - call icepack_init_zsalinity(ncat=ncat, nblyr=nblyr, & - ntrcr_o = ntrcr_o, & - Rayleigh_criteria = RayleighC, & - Rayleigh_real = RayleighR, & - trcrn_bgc = trcrn_bgc, & - nt_bgc_S = nt_bgc_S, & - sss = sss(i)) - Rayleigh_real (i) = RayleighR - Rayleigh_criteria(i) = RayleighC - do n = 1,ncat - do k = 1, nblyr - trcrn(i,nt_bgc_S+k-1,n) = trcrn_bgc(nt_bgc_S-1+k-ntrcr_o,n) - enddo - enddo - enddo ! i - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & - file=__FILE__, line=__LINE__) - endif ! solve_zsal - !----------------------------------------------------------------- ! biogeochemistry initialization !----------------------------------------------------------------- @@ -589,7 +553,7 @@ subroutine init_zbgc ntrcr, nbtrcr, nbtrcr_sw, & ntrcr_o, nt_fbri, & nt_bgc_Nit, nt_bgc_Am, nt_bgc_Sil, & - nt_bgc_DMS, nt_bgc_PON, nt_bgc_S, & + nt_bgc_DMS, nt_bgc_PON, & nt_bgc_DMSPp, nt_bgc_DMSPd, & nt_zbgc_frac, nlt_chl_sw, & nlt_bgc_Nit, nlt_bgc_Am, nlt_bgc_Sil, & @@ -1027,45 +991,27 @@ subroutine init_zbgc !----------------------------------------------------------------- ! resolve conflicts + !----------------------------------------------------------------- + !----------------------------------------------------------------- ! zsalinity and brine !----------------------------------------------------------------- -#ifdef UNDEPRECATE_ZSAL - if (solve_zsal .and. TRZS == 0) then - write(nu_diag,*) 'WARNING: solve_zsal=T but 0 zsalinity tracers' - write(nu_diag,*) 'WARNING: setting solve_zsal = F' - solve_zsal = .false. - elseif (solve_zsal .and. nblyr < 1) then - write(nu_diag,*) 'WARNING: solve_zsal=T but 0 zsalinity tracers' - write(nu_diag,*) 'WARNING: setting solve_zsal = F' - solve_zsal = .false. - endif - - if (solve_zsal .and. ((.not. tr_brine) .or. (ktherm /= 1))) then - write(nu_diag,*) 'WARNING: solve_zsal needs tr_brine=T and ktherm=1' - write(nu_diag,*) 'WARNING: setting tr_brine=T and ktherm=1' - tr_brine = .true. - ktherm = 1 - endif -#else if (solve_zsal) then - write(nu_diag,*) 'WARNING: zsalinity is being deprecated' + call icedrv_system_abort(string=subname//'ERROR: zsalinity is deprecated ', & + file=__FILE__, line=__LINE__) endif -#endif if (tr_brine .and. TRBRI == 0 ) then write(nu_diag,*) & 'WARNING: tr_brine=T but no brine height compiled' write(nu_diag,*) & - 'WARNING: setting solve_zsal and tr_brine = F' - solve_zsal = .false. + 'WARNING: setting tr_brine = F' tr_brine = .false. elseif (tr_brine .and. nblyr < 1 ) then write(nu_diag,*) & 'WARNING: tr_brine=T but no biology layers compiled' write(nu_diag,*) & - 'WARNING: setting solve_zsal and tr_brine = F' - solve_zsal = .false. + 'WARNING: setting tr_brine = F' tr_brine = .false. endif @@ -1073,11 +1019,6 @@ subroutine init_zbgc if (tr_brine) then write(nu_diag,1005) ' phi_snow = ', phi_snow endif - if (solve_zsal) then - write(nu_diag,1010) ' solve_zsal = ', solve_zsal - write(nu_diag,1000) ' grid_oS = ', grid_oS - write(nu_diag,1005) ' l_skS = ', l_skS - endif !----------------------------------------------------------------- ! biogeochemistry @@ -1243,7 +1184,7 @@ subroutine init_zbgc fr_dFe_in=fr_dFe, k_nitrif_in=k_nitrif, t_iron_conv_in=t_iron_conv, & max_loss_in=max_loss, max_dfe_doc1_in=max_dfe_doc1, fr_resp_in=fr_resp, & fr_resp_s_in=fr_resp_s, y_sk_DMS_in=y_sk_DMS, t_sk_conv_in=t_sk_conv, & - t_sk_ox_in=t_sk_ox, modal_aero_in=modal_aero, solve_zsal_in=solve_zsal) + t_sk_ox_in=t_sk_ox, modal_aero_in=modal_aero) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) @@ -1269,21 +1210,6 @@ subroutine init_zbgc ntd = 0 ! if nt_fbri /= 0 then use fbri dependency if (nt_fbri == 0) ntd = -1 ! otherwise make tracers depend on ice volume - nt_bgc_S = 0 - if (solve_zsal) then ! .true. only if tr_brine = .true. - nt_bgc_S = ntrcr + 1 - ntrcr = ntrcr + nblyr - do k = 1,nblyr - trcr_depend(nt_bgc_S + k - 1) = 2 + nt_fbri + ntd - trcr_base (nt_bgc_S,1) = c0 ! default: ice area - trcr_base (nt_bgc_S,2) = c1 - trcr_base (nt_bgc_S,3) = c0 - n_trcr_strata(nt_bgc_S) = 1 - nt_strata(nt_bgc_S,1) = nt_fbri - nt_strata(nt_bgc_S,2) = 0 - enddo - endif - !----------------------------------------------------------------- ! biogeochemistry !----------------------------------------------------------------- @@ -1808,7 +1734,7 @@ subroutine init_zbgc call icepack_init_tracer_indices( & nt_fbri_in=nt_fbri, & nt_bgc_Nit_in=nt_bgc_Nit, nt_bgc_Am_in=nt_bgc_Am, nt_bgc_Sil_in=nt_bgc_Sil, & - nt_bgc_DMS_in=nt_bgc_DMS, nt_bgc_PON_in=nt_bgc_PON, nt_bgc_S_in=nt_bgc_S, & + nt_bgc_DMS_in=nt_bgc_DMS, nt_bgc_PON_in=nt_bgc_PON, & nt_bgc_N_in=nt_bgc_N, nt_bgc_chl_in=nt_bgc_chl, nt_bgc_hum_in=nt_bgc_hum, & nt_bgc_DOC_in=nt_bgc_DOC, nt_bgc_DON_in=nt_bgc_DON, nt_bgc_DIC_in=nt_bgc_DIC, & nt_zaero_in=nt_zaero, nt_bgc_DMSPp_in=nt_bgc_DMSPp, nt_bgc_DMSPd_in=nt_bgc_DMSPd, & diff --git a/configuration/driver/icedrv_restart.F90 b/configuration/driver/icedrv_restart.F90 index 7202da455..ed1cf45aa 100644 --- a/configuration/driver/icedrv_restart.F90 +++ b/configuration/driver/icedrv_restart.F90 @@ -90,7 +90,7 @@ subroutine dumpfile logical (kind=log_kind) :: & tr_iage, tr_FY, tr_lvl, tr_iso, tr_aero, tr_brine, & tr_pond_topo, tr_pond_lvl, tr_snow, tr_fsd -! solve_zsal, skl_bgc, z_tracers +! skl_bgc, z_tracers integer (kind=int_kind) :: dims(2) @@ -118,8 +118,7 @@ subroutine dumpfile tr_brine_out=tr_brine, & tr_pond_topo_out=tr_pond_topo, & tr_pond_lvl_out=tr_pond_lvl,tr_snow_out=tr_snow,tr_fsd_out=tr_fsd) -! call icepack_query_parameters(solve_zsal_out=solve_zsal, & -! skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) +! call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -226,7 +225,7 @@ subroutine dumpfile if (tr_brine) call write_restart_hbrine(dims) ! brine height if (tr_fsd) call write_restart_fsd(dims) ! floe size distribution ! called in icedrv_RunMod.F90 to prevent circular dependencies - ! if (solve_zsal .or. skl_bgc .or. z_tracers) & + ! if (skl_bgc .or. z_tracers) & ! call write_restart_bgc ! biogeochemistry ! Close netcdf file diff --git a/configuration/driver/icedrv_restart_bgc.F90 b/configuration/driver/icedrv_restart_bgc.F90 index c5859c898..688b833b1 100644 --- a/configuration/driver/icedrv_restart_bgc.F90 +++ b/configuration/driver/icedrv_restart_bgc.F90 @@ -34,7 +34,6 @@ module icedrv_restart_bgc subroutine write_restart_bgc() - use icedrv_arrays_column, only: Rayleigh_criteria, Rayleigh_real use icedrv_domain_size, only: n_algae, n_doc, n_dic use icedrv_domain_size, only: n_don, n_zaero, n_fed, n_fep use icedrv_flux, only: sss, nit, amm, sil, dmsp, dms, algalN @@ -48,13 +47,13 @@ subroutine write_restart_bgc() i, k , & ! horizontal, vertical indices mm ! n_algae - logical (kind=log_kind) :: skl_bgc, solve_zsal, z_tracers + logical (kind=log_kind) :: skl_bgc, z_tracers integer (kind=int_kind) :: nbtrcr logical (kind=log_kind) :: tr_bgc_Nit, tr_bgc_Am, tr_bgc_Sil, tr_bgc_hum logical (kind=log_kind) :: tr_bgc_DMS, tr_bgc_PON, tr_bgc_N, tr_bgc_C logical (kind=log_kind) :: tr_bgc_DON, tr_bgc_Fe, tr_zaero , tr_bgc_chl - integer (kind=int_kind) :: nt_bgc_S, nt_bgc_Nit, nt_bgc_AM, nt_bgc_Sil + integer (kind=int_kind) :: nt_bgc_Nit, nt_bgc_AM, nt_bgc_Sil integer (kind=int_kind) :: nt_bgc_hum, nt_bgc_PON integer (kind=int_kind) :: nt_bgc_DMSPp, nt_bgc_DMSPd, nt_bgc_DMS integer (kind=int_kind) :: nt_zbgc_frac @@ -75,7 +74,6 @@ subroutine write_restart_bgc() !----------------------------------------------------------------- call icepack_query_parameters(skl_bgc_out=skl_bgc) - call icepack_query_parameters(solve_zsal_out=solve_zsal) call icepack_query_parameters(z_tracers_out=z_tracers) call icepack_query_tracer_sizes(nbtrcr_out=nbtrcr) call icepack_query_tracer_flags(tr_bgc_Nit_out=tr_bgc_Nit, & @@ -85,7 +83,7 @@ subroutine write_restart_bgc() tr_bgc_N_out=tr_bgc_N, tr_bgc_C_out=tr_bgc_C, & tr_bgc_DON_out=tr_bgc_DON, tr_bgc_Fe_out=tr_bgc_Fe, & tr_zaero_out=tr_zaero, tr_bgc_chl_out=tr_bgc_chl) - call icepack_query_tracer_indices( nt_bgc_S_out=nt_bgc_S, & + call icepack_query_tracer_indices( & nt_bgc_N_out=nt_bgc_N, nt_bgc_AM_out=nt_bgc_AM, & nt_bgc_chl_out=nt_bgc_chl, nt_bgc_C_out=nt_bgc_C, & nt_bgc_DOC_out=nt_bgc_DOC, & @@ -100,29 +98,6 @@ subroutine write_restart_bgc() if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - !----------------------------------------------------------------- - ! Salinity and extras - !----------------------------------------------------------------- - if (solve_zsal) then - - do k = 1, nblyr - call write_restart_field(nu_dump,trcrn(:,nt_bgc_S+k-1,:),ncat) - enddo - - call write_restart_field(nu_dump,sss,1) - - do i = 1, nx - if (Rayleigh_criteria(i)) then - Rayleigh_real (i) = c1 - elseif (.NOT. Rayleigh_criteria(i)) then - Rayleigh_real (i) = c0 - endif - enddo - - call write_restart_field(nu_dump,Rayleigh_real,1) - - endif ! solve_zsal - !----------------------------------------------------------------- ! Skeletal layer BGC !----------------------------------------------------------------- @@ -334,7 +309,6 @@ end subroutine write_restart_bgc subroutine read_restart_bgc() - use icedrv_arrays_column, only: Rayleigh_real, Rayleigh_criteria use icedrv_domain_size, only: n_algae, n_doc, n_dic use icedrv_domain_size, only: n_don, n_zaero, n_fed, n_fep use icedrv_flux, only: sss, nit, amm, sil, dmsp, dms, algalN @@ -348,13 +322,13 @@ subroutine read_restart_bgc() i, k, & ! indices mm ! n_algae - logical (kind=log_kind) :: skl_bgc, solve_zsal, z_tracers + logical (kind=log_kind) :: skl_bgc, z_tracers integer (kind=int_kind) :: nbtrcr logical (kind=log_kind) :: tr_bgc_Nit, tr_bgc_Am, tr_bgc_Sil, tr_bgc_hum logical (kind=log_kind) :: tr_bgc_DMS, tr_bgc_PON, tr_bgc_N, tr_bgc_C logical (kind=log_kind) :: tr_bgc_DON, tr_bgc_Fe, tr_zaero , tr_bgc_chl - integer (kind=int_kind) :: nt_bgc_S, nt_bgc_Nit, nt_bgc_AM, nt_bgc_Sil + integer (kind=int_kind) :: nt_bgc_Nit, nt_bgc_AM, nt_bgc_Sil integer (kind=int_kind) :: nt_bgc_hum, nt_bgc_PON integer (kind=int_kind) :: nt_bgc_DMSPp, nt_bgc_DMSPd, nt_bgc_DMS integer (kind=int_kind) :: nt_zbgc_frac @@ -375,7 +349,6 @@ subroutine read_restart_bgc() !----------------------------------------------------------------- call icepack_query_parameters(skl_bgc_out=skl_bgc) - call icepack_query_parameters(solve_zsal_out=solve_zsal) call icepack_query_parameters(z_tracers_out=z_tracers) call icepack_query_tracer_sizes(nbtrcr_out=nbtrcr) @@ -387,7 +360,7 @@ subroutine read_restart_bgc() tr_bgc_N_out=tr_bgc_N, tr_bgc_C_out=tr_bgc_C, & tr_bgc_DON_out=tr_bgc_DON, tr_bgc_Fe_out=tr_bgc_Fe, & tr_zaero_out=tr_zaero, tr_bgc_chl_out=tr_bgc_chl) - call icepack_query_tracer_indices( nt_bgc_S_out=nt_bgc_S, & + call icepack_query_tracer_indices( & nt_bgc_N_out=nt_bgc_N, nt_bgc_AM_out=nt_bgc_AM, & nt_bgc_chl_out=nt_bgc_chl, nt_bgc_C_out=nt_bgc_C, & nt_bgc_DOC_out=nt_bgc_DOC, & @@ -402,32 +375,6 @@ subroutine read_restart_bgc() if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) - !----------------------------------------------------------------- - ! Salinity and extras - !----------------------------------------------------------------- - - if (solve_zsal) then - - write(nu_diag,*)'zSalinity restart' - do k = 1, nblyr - call read_restart_field(nu_restart,trcrn(:,nt_bgc_S+k-1,:),ncat) - enddo - - write(nu_diag,*) 'min/max sea surface salinity' - call read_restart_field(nu_restart,sss,1) - write(nu_diag,*) 'min/max Rayleigh' - call read_restart_field(nu_restart,Rayleigh_real,1) - - do i = 1, nx - if (Rayleigh_real (i) .GE. c1) then - Rayleigh_criteria (i) = .true. - elseif (Rayleigh_real (i) < c1) then - Rayleigh_criteria (i) = .false. - endif - enddo - - endif ! solve_zsal - !----------------------------------------------------------------- ! Skeletal Layer BGC !----------------------------------------------------------------- diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index 1d7b6ba6c..7fff4ab72 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -427,7 +427,7 @@ end subroutine step_therm1 subroutine step_therm2 (dt) - use icedrv_arrays_column, only: hin_max, fzsal, ocean_bio, & + use icedrv_arrays_column, only: hin_max, ocean_bio, & wave_sig_ht, wave_spectrum, & wavefreq, dwavefreq, & floe_rad_c, floe_binwidth, & @@ -509,7 +509,6 @@ subroutine step_therm2 (dt) bgrid=bgrid, cgrid=cgrid, & igrid=igrid, faero_ocn=faero_ocn(i,:), & first_ice=first_ice(i,:), & - fzsal=fzsal(i), & flux_bio=flux_bio(i,1:nbtrcr), & ocean_bio=ocean_bio(i,1:nbtrcr), & frazil_diag=frazil_diag(i), & @@ -814,7 +813,7 @@ end subroutine step_lateral_flux_scm subroutine step_dyn_ridge (dt, ndtd) - use icedrv_arrays_column, only: hin_max, fzsal, first_ice + use icedrv_arrays_column, only: hin_max, first_ice use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, nx use icedrv_flux, only: rdg_conv, rdg_shear, dardg1dt, dardg2dt use icedrv_flux, only: dvirdgdt, opening, closing, fpond, fresh, fhocn @@ -891,7 +890,7 @@ subroutine step_dyn_ridge (dt, ndtd) dvirdgndt=dvirdgndt(i,:), & araftn=araftn(i,:), vraftn=vraftn(i,:), & aice=aice(i), fsalt=fsalt(i), & - first_ice=first_ice(i,:), fzsal=fzsal(i), & + first_ice=first_ice(i,:), & flux_bio=flux_bio(i,1:nbtrcr), & closing=closing(i) ) @@ -934,7 +933,7 @@ subroutine step_dyn_ridge (dt, ndtd) dvirdgndt=dvirdgndt(i,:), & araftn=araftn(i,:), vraftn=vraftn(i,:), & aice=aice(i), fsalt=fsalt(i), & - first_ice=first_ice(i,:), fzsal=fzsal(i), & + first_ice=first_ice(i,:), & flux_bio=flux_bio(i,1:nbtrcr)) endif ! tmask @@ -1337,12 +1336,12 @@ end subroutine ocean_mixed_layer subroutine biogeochemistry (dt) use icedrv_arrays_column, only: upNO, upNH, iDi, iki, zfswin - use icedrv_arrays_column, only: zsal_tot, darcy_V, grow_net + use icedrv_arrays_column, only: darcy_V, grow_net use icedrv_arrays_column, only: PP_net, hbri,dhbr_bot, dhbr_top, Zoo use icedrv_arrays_column, only: fbio_snoice, fbio_atmice, ocean_bio use icedrv_arrays_column, only: first_ice, fswpenln, bphi, bTiz, ice_bio_net - use icedrv_arrays_column, only: snow_bio_net, fswthrun, Rayleigh_criteria - use icedrv_arrays_column, only: ocean_bio_all, sice_rho, fzsal, fzsal_g + use icedrv_arrays_column, only: snow_bio_net, fswthrun + use icedrv_arrays_column, only: ocean_bio_all, sice_rho use icedrv_arrays_column, only: bgrid, igrid, icgrid, cgrid use icepack_intfc, only: icepack_biogeochemistry, icepack_load_ocean_bio_array use icedrv_domain_size, only: nblyr, nilyr, nslyr, n_algae, n_zaero, ncat @@ -1462,7 +1461,6 @@ subroutine biogeochemistry (dt) iDi = iDi(i,:,:), & iki = iki(i,:,:), & zfswin = zfswin(i,:,:), & - zsal_tot = zsal_tot(i), & darcy_V = darcy_V(i,:), & grow_net = grow_net(i), & PP_net = PP_net(i), & @@ -1480,10 +1478,7 @@ subroutine biogeochemistry (dt) ice_bio_net = ice_bio_net(i,1:nbtrcr), & snow_bio_net = snow_bio_net(i,1:nbtrcr), & fswthrun = fswthrun(i,:), & - Rayleigh_criteria = Rayleigh_criteria(i), & sice_rho = sice_rho(i,:), & - fzsal = fzsal(i), & - fzsal_g = fzsal_g(i), & meltbn = meltbn(i,:), & melttn = melttn(i,:), & congeln = congeln(i,:), & diff --git a/configuration/scripts/icepack.build b/configuration/scripts/icepack.build index 43a47661b..e15d41550 100755 --- a/configuration/scripts/icepack.build +++ b/configuration/scripts/icepack.build @@ -24,7 +24,7 @@ endif if !(-d ${ICE_OBJDIR}) mkdir -p ${ICE_OBJDIR} cd ${ICE_OBJDIR} -setenv ICE_CPPDEFS "${ICE_CPPDEFS} -DNXGLOB=${ICE_NXGLOB} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRSNOW=${TRSNOW} -DTRBRI=${TRBRI} -DNTRISO=${NTRISO} -DNTRAERO=${NTRAERO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} " +setenv ICE_CPPDEFS "${ICE_CPPDEFS} -DNXGLOB=${ICE_NXGLOB} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRSNOW=${TRSNOW} -DTRBRI=${TRBRI} -DNTRISO=${NTRISO} -DNTRAERO=${NTRAERO} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} " if (${ICE_IOTYPE} == 'netcdf') then setenv ICE_CPPDEFS "${ICE_CPPDEFS} -DUSE_NETCDF" diff --git a/configuration/scripts/icepack.settings b/configuration/scripts/icepack.settings index dab7276fc..38262b19f 100755 --- a/configuration/scripts/icepack.settings +++ b/configuration/scripts/icepack.settings @@ -55,8 +55,6 @@ setenv NTRAERO 1 # number of aerosol tracers # CESM uses 3 aerosol tracers setenv NTRISO 0 # number of isotopes (up to max_iso) setenv TRBRI 0 # set to 1 for brine height tracer -setenv TRZS 0 # set to 1 for zsalinity tracer - # (needs TRBRI = 1) setenv TRBGCS 0 # set to 1 for skeletal layer tracers # (needs TRBGCZ = 0) setenv TRBGCZ 0 # set to 1 for zbgc tracers diff --git a/configuration/scripts/options/set_env.bgcispol b/configuration/scripts/options/set_env.bgcispol index 243fbeaf1..c38635a70 100644 --- a/configuration/scripts/options/set_env.bgcispol +++ b/configuration/scripts/options/set_env.bgcispol @@ -12,8 +12,6 @@ setenv NTRAERO 0 # number of aerosol tracers # (up to max_aero in ice_domain_size.F90) # CESM uses 3 aerosol tracers setenv TRBRI 1 # set to 1 for brine height tracer -setenv TRZS 0 # set to 1 for zsalinity tracer - # (needs TRBRI = 1) setenv TRBGCS 0 # set to 1 for skeletal layer tracers # (needs TRBGCZ = 0) setenv TRBGCZ 1 # set to 1 for zbgc tracers diff --git a/configuration/scripts/options/set_env.bgcnice b/configuration/scripts/options/set_env.bgcnice index 9b266b256..1d6efc405 100644 --- a/configuration/scripts/options/set_env.bgcnice +++ b/configuration/scripts/options/set_env.bgcnice @@ -12,8 +12,6 @@ setenv NTRAERO 0 # number of aerosol tracers # (up to max_aero in ice_domain_size.F90) # CESM uses 3 aerosol tracers setenv TRBRI 1 # set to 1 for brine height tracer -setenv TRZS 0 # set to 1 for zsalinity tracer - # (needs TRBRI = 1) setenv TRBGCS 0 # set to 1 for skeletal layer tracers # (needs TRBGCZ = 0) setenv TRBGCZ 1 # set to 1 for zbgc tracers diff --git a/configuration/scripts/options/set_env.bgcsklnice b/configuration/scripts/options/set_env.bgcsklnice index 1b1069862..985bacbf4 100644 --- a/configuration/scripts/options/set_env.bgcsklnice +++ b/configuration/scripts/options/set_env.bgcsklnice @@ -12,8 +12,6 @@ setenv NTRAERO 0 # number of aerosol tracers # (up to max_aero in ice_domain_size.F90) # CESM uses 3 aerosol tracers setenv TRBRI 1 # set to 1 for brine height tracer -setenv TRZS 0 # set to 1 for zsalinity tracer - # (needs TRBRI = 1) setenv TRBGCS 1 # set to 1 for skeletal layer tracers # (needs TRBGCZ = 0) setenv TRBGCZ 0 # set to 1 for zbgc tracers diff --git a/configuration/scripts/options/set_env.zsal b/configuration/scripts/options/set_env.zsal deleted file mode 100644 index a9f52bb6c..000000000 --- a/configuration/scripts/options/set_env.zsal +++ /dev/null @@ -1,10 +0,0 @@ -### Layers -setenv NICELYR 7 # number of vertical layers in the ice -setenv NSNWLYR 1 # number of vertical layers in the snow -setenv NICECAT 5 # number of ice thickness categories - -### Tracers # match icepack_in tracer_nml to conserve memory -setenv TRBRI 1 # set to 1 for brine height tracer -setenv TRZS 1 # set to 1 for zsalinity tracer - # (needs TRBRI = 1) -setenv NBGCLYR 7 # number of zbgc layers diff --git a/configuration/scripts/options/set_nml.zsal b/configuration/scripts/options/set_nml.zsal deleted file mode 100644 index 2c7eda3ff..000000000 --- a/configuration/scripts/options/set_nml.zsal +++ /dev/null @@ -1,5 +0,0 @@ -ktherm = 1 -sw_redist = .true. -tfrz_option = 'linear_salt' -tr_brine = .true. -solve_zsal = .true. diff --git a/doc/source/developer_guide/dg_adding_diagnostics.rst b/doc/source/developer_guide/dg_adding_diagnostics.rst index ef5960830..2f2f694c0 100755 --- a/doc/source/developer_guide/dg_adding_diagnostics.rst +++ b/doc/source/developer_guide/dg_adding_diagnostics.rst @@ -55,7 +55,7 @@ Icepack produces separate ASCII (text) log output for four cells, each with a di #. For BGC variables, edit **icedrv\_diagnostics\_bgc.F90**: - If the variable is already defined within the code, then add it to a "use" statement in the subroutine - ``hbrine_diags`` or ``bgc_diags`` or ``zsal_diags`` and follow a similar procedure for state variables as above. + ``hbrine_diags`` or ``bgc_diags`` and follow a similar procedure for state variables as above. - Note that the BGC needs to be activated and the particular tracer turned on. diff --git a/doc/source/developer_guide/dg_col_phys.rst b/doc/source/developer_guide/dg_col_phys.rst index 7fd5389ea..df5a182f1 100755 --- a/doc/source/developer_guide/dg_col_phys.rst +++ b/doc/source/developer_guide/dg_col_phys.rst @@ -42,7 +42,6 @@ The column physics source code contains the following files | **icepack_wavefracspec.F90** wave impact on sea ice | **icepack_zbgc.F90** driver for ice biogeochemistry and brine tracer motion | **icepack_zbgc_shared.F90** parameters and shared code for biogeochemistry and brine height -| **icepack_zsalinity.F90** vertical salinity parameterization of :cite:`Jeffery11` Coding Standard diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include index 2c641afa4..cfedb8697 100644 --- a/doc/source/user_guide/interfaces.include +++ b/doc/source/user_guide/interfaces.include @@ -110,7 +110,8 @@ icepack_init_zsalinity ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: fortran - ! Initialize zSalinity + ! **DEPRECATED**, all code removed + ! Interface provided for backwards compatibility subroutine icepack_init_zsalinity(nblyr,ntrcr_o, Rayleigh_criteria, & Rayleigh_real, trcrn_bgc, nt_bgc_S, ncat, sss) @@ -302,7 +303,7 @@ icepack_intfc.F90 use icepack_shortwave, only: icepack_step_radiation use icepack_brine, only: icepack_init_hbrine - use icepack_brine, only: icepack_init_zsalinity + use icepack_brine, only: icepack_init_zsalinity ! deprecated use icepack_zbgc , only: icepack_init_bgc use icepack_zbgc , only: icepack_init_zbgc @@ -562,8 +563,10 @@ icepack_step_ridge fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) - fhocn , & ! net heat flux to ocean (W/m^2) - fzsal ! zsalinity flux to ocean(kg/m^2/s) + fhocn ! net heat flux to ocean (W/m^2) + + real (kind=dbl_kind), intent(inout), optional :: & + fzsal ! zsalinity flux to ocean(kg/m^2/s) (deprecated) real (kind=dbl_kind), intent(inout), optional :: & closing ! rate of closing due to divergence/shear (1/s) @@ -2008,7 +2011,7 @@ icepack_step_therm2 Tf , & ! freezing temperature (C) sss , & ! sea surface salinity (ppt) rside , & ! fraction of ice that melts laterally - frzmlt ! freezing/melting potential (W/m^2) + frzmlt ! freezing/melting potential (W/m^2) integer (kind=int_kind), dimension (:), intent(in) :: & trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon @@ -2043,11 +2046,13 @@ icepack_step_therm2 fresh , & ! fresh water flux to ocean (kg/m^2/s) fsalt , & ! salt flux to ocean (kg/m^2/s) fhocn , & ! net heat flux to ocean (W/m^2) - fzsal , & ! salt flux to ocean from zsalinity (kg/m^2/s) meltl , & ! lateral ice melt (m/step-->cm/day) frazil , & ! frazil ice growth (m/step-->cm/day) frazil_diag ! frazil ice growth diagnostic (m/step-->cm/day) + real (kind=dbl_kind), intent(inout), optional :: & + fzsal ! salt flux to ocean from zsalinity (kg/m^2/s) (deprecated) + real (kind=dbl_kind), intent(in), optional :: & wlat ! lateral melt rate (m/s) @@ -2697,7 +2702,7 @@ icepack_init_tracer_indices nlt_bgc_hum_in,& ! nlt_bgc_PON_in,& ! zooplankton and detritus nt_zbgc_frac_in,&! fraction of tracer in the mobile phase - nt_bgc_S_in, & ! Bulk salinity in fraction ice with dynamic salinity (Bio grid)) + nt_bgc_S_in, & ! (deprecated, was related to zsalinity) nlt_chl_sw_in ! points to total chla in trcrn_sw integer (kind=int_kind), dimension(:), intent(in), optional :: & @@ -2801,7 +2806,7 @@ icepack_query_tracer_indices nlt_bgc_hum_out,& ! nlt_bgc_PON_out,& ! zooplankton and detritus nt_zbgc_frac_out,&! fraction of tracer in the mobile phase - nt_bgc_S_out, & ! Bulk salinity in fraction ice with dynamic salinity (Bio grid)) + nt_bgc_S_out, & ! (deprecated, was related to zsalinity) nlt_chl_sw_out ! points to total chla in trcrn_sw integer (kind=int_kind), dimension(:), intent(out), optional :: & @@ -3342,14 +3347,18 @@ icepack_biogeochemistry grow_net , & ! Specific growth rate (/s) per grid cell PP_net , & ! Total production (mg C/m^2/s) per grid cell hbri , & ! brine height, area-averaged for comparison with hi (m) - zsal_tot , & ! Total ice salinity in per grid cell (g/m^2) - fzsal , & ! Total flux of salt to ocean at time step for conservation - fzsal_g , & ! Total gravity drainage flux upNO , & ! nitrate uptake rate (mmol/m^2/d) times aice - upNH ! ammonium uptake rate (mmol/m^2/d) times aice + upNH ! ammonium uptake rate (mmol/m^2/d) times aice - logical (kind=log_kind), intent(inout) :: & - Rayleigh_criteria ! .true. means Ra_c was reached + real (kind=dbl_kind), intent(inout), optional :: & + zsal_tot ! Total ice salinity in per grid cell (g/m^2) (deprecated) + + real (kind=dbl_kind), intent(inout), optional :: & + fzsal , & ! Total flux of salt to ocean at time step for conservation (deprecated) + fzsal_g ! Total gravity drainage flux (deprecated) + + logical (kind=log_kind), intent(inout), optional :: & + Rayleigh_criteria ! .true. means Ra_c was reached (deprecated) real (kind=dbl_kind), dimension (:,:), intent(in) :: & fswpenln ! visible SW entering ice layers (W m-2) diff --git a/doc/source/user_guide/lg_protocols.rst b/doc/source/user_guide/lg_protocols.rst index 5c9b06e8d..4ebe0f8a2 100755 --- a/doc/source/user_guide/lg_protocols.rst +++ b/doc/source/user_guide/lg_protocols.rst @@ -147,8 +147,6 @@ tracers understood by Icepack and lists some of their properties. See also :ref +------------+----------+---------------+---------+---------+-----------------------------------------------------------------------------------+ | fbri | optional | tr_brine | 1 | 1 | | +------------+----------+---------------+---------+---------+-----------------------------------------------------------------------------------+ - | bgc_S | optional | | 1 | nblyr | bulk salinity in fraction ice | - +------------+----------+---------------+---------+---------+-----------------------------------------------------------------------------------+ | bgc_N | optional | tr_bgc_N | n_algae | nblyr+3 | nutrients | +------------+----------+---------------+---------+---------+-----------------------------------------------------------------------------------+ | bgc_Nit | optional | | 1 | nblyr+3 | diatoms, phaeocystis, pico/small | diff --git a/doc/source/user_guide/lg_sequence.rst b/doc/source/user_guide/lg_sequence.rst index 09127b8f2..04a72269c 100755 --- a/doc/source/user_guide/lg_sequence.rst +++ b/doc/source/user_guide/lg_sequence.rst @@ -123,10 +123,10 @@ computation and a routine to update the radiation computation:: use icepack_shortwave, only: icepack_prep_radiation use icepack_shortwave, only: icepack_step_radiation -icepack_brine address brine and zsalinity computations:: +icepack_brine addresses brine computations:: use icepack_brine, only: icepack_init_hbrine - use icepack_brine, only: icepack_init_zsalinity + use icepack_brine, only: icepack_init_zsalinity ! DEPRECATED icepack_zbgc contains several public interfaces to support initialization and computation for the skeletal layer bgc and zbgc options:: @@ -214,7 +214,6 @@ of how to do this in practice. The sample below does not include bgc:: call *icepack_init_itd_hist* loop over gridcells call *icepack_step_radiation* - call *icepack_init_zsalinity* end loop over gridcells call *icepack_init_hbrine* loop over gridcells diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 5fcc188d9..fa62e1781 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -436,7 +436,7 @@ zbgc_nml "``f_exude_l``", "real", "fraction of exudation to DOC lipids", "1.0" "``f_exude_s``", "real", "fraction of exudation to DOC saccharids", "1.0" "``grid_o``", "real", "z biology for bottom flux", "5.0" - "``grid_oS``", "real", "z salinity for bottom flux", "5.0" + "``grid_oS``", "real", "DEPRECATED", "" "``grow_Tdep_diatoms``", "real", "temperature dependence growth diatoms per degC", "0.06" "``grow_Tdep_phaeo``", "real", "temperature dependence growth phaeocystis per degC", "0.06" "``grow_Tdep_sp``", "real", "temperature dependence growth small plankton per degC", "0.06" @@ -462,7 +462,7 @@ zbgc_nml "``K_Sil_sp``", "real", "silicate half saturation small plankton mmol/m^3", "0.0" "``kn_bac_protein``", "real", "bacterial degradation of DON per day", "0.03" "``l_sk``", "real", "characteristic diffusive scale in m", "7.0" - "``l_skS``", "real", "z salinity characteristic diffusive scale in m", "7.0" + "``l_skS``", "real", "DEPRECATED", "" "``max_dfe_doc1``", "real", "max ratio of dFe to saccharides in the ice in nm Fe / muM C", "0.2" "``max_loss``", "real", "restrict uptake to percent of remaining value", "0.9" "``modal_aero``", "logical", "modal aersols", "``.false.``" @@ -546,7 +546,6 @@ zbgc_nml .. "``grid_o_t``", "real", "z biology for top flux", "5.0" .. "``restart_bgc``", "logical", "restart tracer values from file", "``.false.``" .. "``restart_hbrine``", "logical", "", "``.false.``" -.. "``restart_zsal``", "logical", "", "``.false.``" .. "``solve_zsal``", "logical", "update salinity tracer profile", "``.false.``" .. "TRZS", "0,1", "zsalinity tracer, needs TRBRI=1", "0" From 23b6c1272b50d42cad7928ffe0005d6ee673dee9 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Wed, 30 Aug 2023 12:55:49 -0700 Subject: [PATCH 2/2] Corrects thin ice/snow treatment of enthalpy and other tracers (#454) * Corrects thin ice/snow treatment of enthalpy and other tracers Adopted from https://github.com/E3SM-Project/E3SM/pull/5630 "This fix redistributes enthalpy and other tracers evenly in snow and ice when their respective thicknesses are < 1e-15 . Previously, these tracers were zero'd non-conservatively. Also corrects a bug in the zeroing of snow thicknesses, and removes snow in thickness_changes if the ice vanishes." * Clean up modifications --- columnphysics/icepack_therm_shared.F90 | 76 +++++++++++++++--------- columnphysics/icepack_therm_vertical.F90 | 5 +- 2 files changed, 52 insertions(+), 29 deletions(-) diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90 index ad00db2ac..53113c08d 100644 --- a/columnphysics/icepack_therm_shared.F90 +++ b/columnphysics/icepack_therm_shared.F90 @@ -497,7 +497,8 @@ subroutine adjust_enthalpy (nlyr, & hovlp ! overlap between old and new layers (m) real (kind=dbl_kind) :: & - rhlyr ! 1./hlyr + rhlyr, & ! 1./hlyr + qtot ! total h*q in the column real (kind=dbl_kind), dimension (nlyr) :: & hq ! h * q for a layer @@ -509,36 +510,55 @@ subroutine adjust_enthalpy (nlyr, & !----------------------------------------------------------------- rhlyr = c0 - if (hn > puny) rhlyr = c1 / hlyr - - !----------------------------------------------------------------- - ! Compute h*q for new layers (k2) given overlap with old layers (k1) - !----------------------------------------------------------------- - - do k2 = 1, nlyr - hq(k2) = c0 - enddo ! k - k1 = 1 - k2 = 1 - do while (k1 <= nlyr .and. k2 <= nlyr) - hovlp = min (z1(k1+1), z2(k2+1)) & - - max (z1(k1), z2(k2)) - hovlp = max (hovlp, c0) - hq(k2) = hq(k2) + hovlp*qn(k1) - if (z1(k1+1) > z2(k2+1)) then - k2 = k2 + 1 + if (hn > puny) then + rhlyr = c1 / hlyr + + !----------------------------------------------------------------- + ! Compute h*q for new layers (k2) given overlap with old layers (k1) + !----------------------------------------------------------------- + + do k2 = 1, nlyr + hq(k2) = c0 + enddo ! k + k1 = 1 + k2 = 1 + do while (k1 <= nlyr .and. k2 <= nlyr) + hovlp = min (z1(k1+1), z2(k2+1)) & + - max (z1(k1), z2(k2)) + hovlp = max (hovlp, c0) + hq(k2) = hq(k2) + hovlp*qn(k1) + if (z1(k1+1) > z2(k2+1)) then + k2 = k2 + 1 + else + k1 = k1 + 1 + endif + enddo ! while + + !----------------------------------------------------------------- + ! Compute new enthalpies. + !----------------------------------------------------------------- + + do k = 1, nlyr + qn(k) = hq(k) * rhlyr + enddo ! k + + else + + qtot = c0 + do k = 1, nlyr + qtot = qtot + qn(k) * (z1(k+1)-z1(k)) + enddo + if (hn > c0) then + do k = 1, nlyr + qn(k) = qtot/hn + enddo else - k1 = k1 + 1 + do k = 1, nlyr + qn(k) = c0 + enddo endif - enddo ! while - !----------------------------------------------------------------- - ! Compute new enthalpies. - !----------------------------------------------------------------- - - do k = 1, nlyr - qn(k) = hq(k) * rhlyr - enddo ! k + endif end subroutine adjust_enthalpy diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 09328c5d3..78adc7f26 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -1706,7 +1706,7 @@ subroutine thickness_changes (nilyr, nslyr, & !----------------------------------------------------------------- if (ktherm == 2) then - if (hsn <= puny) then + if (hsn <= puny .or. hin <= c0) then do k = 1, nslyr fhocnn = fhocnn & + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt) @@ -1715,8 +1715,11 @@ subroutine thickness_changes (nilyr, nslyr, & meltsliq = meltsliq + massice(k) ! add to meltponds smice(k) = rhos smliq(k) = c0 + rsnw(k) = rsnw_fall endif enddo + melts = melts + hsn + hsn = c0 hslyr = c0 endif endif