SUBROUTINE ice_sal_adv(nlay_i,kideb,kiut) !!------------------------------------------------------------------ !! *** ROUTINE ice_sal_adv *** !! !! ** Purpose : !! This routine computes new salinities in the ice !! !! ** Method : Vertical salinity profile computation !! Resolves brine transport equation !! !! ** Steps !! !! ** Arguments !! !! ** Inputs / Outputs !! !! ** External !! !! ** References : Vancop. et al., 2008 !! !! ** History : !! (06-2003) Martin Vancop. LIM1D !! (06-2008) Martin Vancop. BIO-LIM !! (09-2008) Martin Vancop. Explicit gravity drainage !! !!------------------------------------------------------------------ USE lib_fortran INCLUDE 'type.com' INCLUDE 'para.com' INCLUDE 'const.com' INCLUDE 'ice.com' INCLUDE 'thermo.com' REAL(8), DIMENSION(nlay_i) :: & z_ms_i , !: mass of salt times thickness & z_sbr_i !: brine salinity REAL(8), DIMENSION(nlay_i) :: !: dummy factors for tracer equation & za , !: all & zb , !: gravity drainage & zc , !: upward advective flow & zRae , !: effective Ra & ze , !: downward advective flow & zind , !: independent term in the tridiag system & zindtbis , !: & zdiagbis !: REAL(8), DIMENSION(nlay_i,3) :: !: dummy factors for tracer equation & ztrid !: tridiagonal matrix REAL(8) :: & zdummy1 , !: dummy factors & zdummy2 , !: & zdummy3 , !: & zswitchs , !: switch for summer drainage & zeps = 1.0e-20 !: numerical limit ! Rayleigh number computation REAL(8) :: & ze_i_min , !: minimum brine volume & zcp , !: temporary scalar for sea ice specific heat & zk , !: temporary scalar for sea ice thermal conductivity & zalphara !: multiplicator for diffusivity REAL(8), DIMENSION(nlay_i) :: & zsigma , !: brine salinity at layer interfaces & zperm , !: permeability & zpermin , !: minimum permeability & zperm_eff , !: minimum permeability & zrhodiff , !: density difference & zlevel , !: height of the water column & zthdiff !: thermal diffusivity REAL(8), DIMENSION(nlay_i+1) :: & z_sbr_int !: brine salinity at layer interfaces INTEGER :: & layer2 , !: layer loop index & indtr !: index of tridiagonal system CHARACTER(len=4) :: & bc = 'conc' !: Boundary condition 'conc' or 'flux' REAL(8) :: & z_ms_i_ini , !: initial mass of salt & z_ms_i_fin , !: final mass of salt & z_fs_b , !: basal flux of salt & z_fs_su , !: surface flux of salt & z_dms_i !: mass variation LOGICAL :: & ln_write , & ln_con , & ln_sal ln_write = .TRUE. ! write outputs ln_con = .TRUE. ! conservation check ln_sal = .TRUE. ! compute salinity variations or not IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' ** ice_sal_adv : ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ ' WRITE(numout,*) ' ln_sal = ', ln_sal WRITE(numout,*) ' ln_grd = ', ln_grd WRITE(numout,*) ' ln_flu = ', ln_flu WRITE(numout,*) ' ln_flo = ', ln_flo WRITE(numout,*) ' c_gravdr = ', c_gravdr WRITE(numout,*) ' c_sbr = ', c_sbr WRITE(numout,*) ' c_perm = ', c_perm WRITE(numout,*) ' c_permeff= ', c_permeff ENDIF ji = 1 IF ( ln_sal ) THEN ! !------------------------------------------------------------------------------| ! 1) Initialization !------------------------------------------------------------------------------| ! IF ( ln_write ) THEN WRITE(numout,*) ' - Initialization ... ' ENDIF ! brine diffusivity diff_br(:) = 0.0 ! Darcy velocity w_adv_br(:) = 0.0 !-------------------- ! Conservation check !-------------------- IF ( ln_con ) THEN CALL ice_sal_column( kideb , kiut , z_ms_i_ini , & s_i_b(1,1:nlay_i), & deltaz_i_phy, nlay_i, .FALSE. ) ENDIF ! ln_con IF ( ln_write ) THEN WRITE(numout,*) ' ' WRITE(numout,*) ' nlay_i : ', nlay_i WRITE(numout,*) ' kideb : ', kideb WRITE(numout,*) ' kiut : ', kiut WRITE(numout,*) ' ' WRITE(numout,*) ' deltaz_i_phy : ', ( deltaz_i_phy(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' z_i_phy : ', ( z_i_phy(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' s_i_b : ', ( s_i_b (ji,layer), & layer = 1, nlay_i ) WRITE(numout,*) ' t_i_b : ', ( t_i_b (ji,layer), & layer = 1, nlay_i ) WRITE(numout,*) ' t_i_int : ', ( t_i_int (ji,layer), & layer = 1, nlay_i+1 ) WRITE(numout,*) ENDIF ! ln_write ! !------------------------------------------------------------------------------| ! 2) Brine salinity at layer mid points and interfaces, brine fraction !------------------------------------------------------------------------------| ! DO layer = 1, nlay_i e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer) & - tpw ) END DO IF ( c_sbr .EQ. 'LIN' ) THEN ! Linear liquidus DO layer = 1, nlay_i + 1 zTc = t_i_int(ji,layer) - tpw z_sbr_int(layer) = - zTc / tmut !--- interfacial value END DO DO layer = 1, nlay_i zTc = t_i_b(ji,layer) - tpw z_sbr_i(layer) = - zTc / tmut !--- mid-point value END DO ENDIF IF ( c_sbr .EQ. 'WEA' ) THEN ! Weast (1971) 3rd order liquidus DO layer = 1, nlay_i + 1 zTc = t_i_int(ji,layer) - tpw z_sbr_int(layer) = -17.6*zTc -0.389*zTc**2. -0.00362*zTc**3. END DO DO layer = 1, nlay_i zTc = t_i_b(ji,layer) - tpw z_sbr_i(layer) = -17.6*zTc -0.389*zTc**2. -0.00362*zTc**3. END DO ENDIF IF ( c_sbr .EQ. 'NTZ' ) THEN ! Notz (2005) 3rd order liquidus DO layer = 1, nlay_i + 1 zTc = t_i_int(ji,layer) - tpw z_sbr_int(layer) = -21.4*zTc - 0.886*zTc**2. - 0.017*zTc**3. END DO DO layer = 1, nlay_i zTc = t_i_b(ji,layer) - tpw z_sbr_i(layer) = -21.4*zTc - 0.886*zTc**2. - 0.017*zTc**3. END DO ENDIF IF ( ln_write ) THEN WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i (layer), & layer = 1, nlay_i ) WRITE(numout,*) ' z_sbr_int : ', ( z_sbr_int (layer), & layer = 1, nlay_i + 1 ) WRITE(numout,*) ' e_i_b : ', ( e_i_b (layer), & layer = 1, nlay_i ) ENDIF ! ln_write ! !------------------------------------------------------------------------------| ! 3) Effective permeability at layer mid-points !------------------------------------------------------------------------------| ! IF ( c_permeff .EQ. 'HAR' ) THEN ! Harmonic mean WRITE(numout,*) ENDIF IF ( c_permeff .EQ. 'MIN' ) THEN ! Minimum permeability DO layer = 1, nlay_i ze_i_min = 99999.0 DO layer2 = layer, nlay_i ze_i_min = MIN( ze_i_min , e_i_b(layer2) ) IF ( c_perm .EQ. 'FRE' ) ! Freitag (1999) & zperm_eff(layer) = 1.995e-8 * ze_i_min**3.1 IF ( c_perm .EQ. 'RJW' ) ! Rees-Jones and Worster (2014) & zperm_eff(layer) = 1.0e-8 * ze_i_min**3 END DO END DO ! layer ENDIF ! !------------------------------------------------------------------------------| ! 4) Rayleigh number at mid-points (see Vancoppenolle et al TCD2013) !------------------------------------------------------------------------------| ! z1 = gpes / ( thdiff_br * visc_br ); DO layer = 1, nlay_i z2 = beta_ocs * ( z_sbr_i(layer) - oce_sal ) ! Delta rho z3 = ht_i_b(ji) - z_i_phy(layer) z4 = zperm_eff(layer); rayleigh(layer) = MAX( z1 * z2 * z3 * z4, 0.0) END DO ! !------------------------------------------------------------------------------| ! 5) Brine diffusivity !------------------------------------------------------------------------------| ! !------------------- ! Brine Diffusivity !------------------- DO layer = 1, nlay_i zalphara = ( TANH( ra_smooth * ( rayleigh(layer) - ra_c ) ) & + 1 ) / 2.0 diff_br(layer) = ( 1.0 - zalphara ) * d_br_mol + & zalphara * ( d_br_tur ) IF ( .NOT. ln_grd ) diff_br(layer) = 0. END DO IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' rayleigh : ', ( rayleigh(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' diff_br : ', ( diff_br(layer), & layer = 1, nlay_i ) WRITE(numout,*) ENDIF ! !------------------------------------------------------------------------------| ! 6) Brine velocity !------------------------------------------------------------------------------| ! ! c_wadv = 'GN' or 'RW' alpha_GN = 1.56e-3 rho_br_GN = 1020. Rc_GN = 1.01 !alpha_GN = 1.0e-3 !Rc_GN = 1.0 w_adv_br(:) = 0.0 zRae(:) = 0. DO layer = 1, nlay_i zRae(layer) = MAX( rayleigh(layer) - Rc_GN, 0.) ! correction from the orig scheme END DO DO layer = 1, nlay_i w_adv_br(layer) = - alpha_GN / rho_br_GN * & SUM ( zRae(1:layer)*deltaz_i_phy(1:layer) ) END DO IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' w_adv_br : ', ( w_adv_br(layer), & layer = 1, nlay_i ) ENDIF ! !------------------------------------------------------------------------------| ! 7) New salinities !------------------------------------------------------------------------------| ! DO layer = 1, nlay_i za(layer) = w_adv_br(layer) * ddtb / deltaz_i_phy(layer) END DO ! first layer sn_i_b(1) = z_sbr_i(1) * ( e_i_b(1) + za(1) ) + & z_sbr_i(2) * ( - za(1) ) ! inner layers DO layer = 2, nlay_i - 1 sn_i_b(layer) = z_sbr_i(layer-1) * ( za(layer)/2. ) + & z_sbr_i(layer) * e_i_b(layer) + & z_sbr_i(layer+1) * ( - za(layer)/2. ) END DO ! lowermost layer sn_i_b(nlay_i) = z_sbr_i(nlay_i-1) * ( za(nlay_i)/2. ) + & z_sbr_i(nlay_i) * ( e_i_b(nlay_i) + & za(nlay_i)/2. ) - za(nlay_i) * oce_sal IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' sn_i_b : ', ( sn_i_b(layer) , & layer = 1, nlay_i ) ENDIF ! !----------------------------------------------------------------------- ! 8) Mass of salt conserved ? !----------------------------------------------------------------------- ! ! Final mass of salt CALL ice_sal_column( kideb , kiut , z_ms_i_fin , & sn_i_b(1:nlay_i), & deltaz_i_phy, nlay_i, .FALSE. ) ! Bottom flux ( positive upwards for conservation routine ) zfb = - e_i_b(nlay_i) * & ( diff_br(nlay_i) * 2.0 / deltaz_i_phy(nlay_i) * & ( z_sbr_i(nlay_i) - oce_sal ) + w_flood * ( z_flood * & oce_sal + ( 1. - z_flood ) * z_sbr_i(nlay_i) ) ) & - e_i_b(nlay_i) * w_flush * z_sbr_i(nlay_i) ! & - qsummer * z_sbr_i(nlay_i) / ddtb fsb = - zfb * rhog / 1000. ! ice-ocean salt flux is positive downwards IF ( ln_write ) THEN WRITE(numout,*) ' fsb : ', fsb WRITE(numout,*) ENDIF ! Surface flux of salt zfsu = 0.0 ! conservation check zerror = 1.0e-15 CALL ice_sal_conserv(kideb,kiut,'ice_sal_adv : ',zerror, & z_ms_i_ini,z_ms_i_fin, & zfb , zfsu , ddtb) ENDIF ! ln_sal ! !------------------------------------------------------------------------------| ! End of la sous-routine WRITE(numout,*) END SUBROUTINE