SUBROUTINE ice_sal_diff(nlay_i,kideb,kiut) !!------------------------------------------------------------------ !! *** ROUTINE ice_sal_diff *** !! !! ** 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 & 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 & zrhodiff , !: density difference & zlevel , !: height of the water column & zthdiff !: thermal diffusivity 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_diff : ' 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 ENDIF WRITE(numout,*) " nlay_i : ", nlay_i IF ( ln_sal ) THEN ! !------------------------------------------------------------------------------| ! 1) Initialization !------------------------------------------------------------------------------| ! IF ( ln_write ) THEN WRITE(numout,*) ' - Initialization ... ' ENDIF DO 10 ji = kideb, kiut ! brine diffusivity diff_br(:) = 0.0 !--------------------------- ! Brine volume and salinity !--------------------------- DO layer = 1, nlay_i e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer) & - tpw ) z_sbr_i(layer) = s_i_b(ji,layer) / e_i_b(layer) END DO !-------------------- ! 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,*) ' 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,*) ' e_i_b : ', ( e_i_b (layer), & layer = 1, nlay_i ) WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i (layer), & layer = 1, nlay_i ) WRITE(numout,*) ENDIF ! ln_write 10 CONTINUE ! !------------------------------------------------------------------------------| ! 2) Rayleigh-number-based diffusivity !------------------------------------------------------------------------------| ! ! Diffusivity is a function of the local Rayleigh number ! see Notz and Worster, JGR 2008 ! Diffusivity, layer represents the interface' ! between layer and layer-1 ' ! IF ( ln_write ) THEN WRITE(numout,*) ' - Rayleigh-number based diffusivity ... ' WRITE(numout,*) ' ' ENDIF DO 20 ji = kideb, kiut !----------------------------------------- ! Brine salinity between layer interfaces !----------------------------------------- DO layer = 1, nlay_i - 1 zdummy1 = t_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * & ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) / & ( z_i_phy(layer+1) - z_i_phy(layer) ) - tpw zsigma(layer) = - zdummy1 / tmut END DO zsigma(nlay_i) = - ( t_i_b(ji,nlay_i) - tpw ) / tmut !-------------------- ! Density difference !-------------------- DO layer = 1, nlay_i zrhodiff(layer) = - beta_ocs * ( oce_sal - zsigma(layer) ) END DO !------------------------------------------ ! Minimum permeability under current level !------------------------------------------ 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) ) zpermin(layer) = 1.0e-17 * ( ( 1000. * ze_i_min )**3.1 ) END DO END DO ! layer !------------------------------------------------ ! length of the water column under current level !------------------------------------------------ DO layer = nlay_i, 1, -1 zlevel(layer) = 0.0 DO layer2 = layer, nlay_i zlevel(layer) = zlevel(layer) + deltaz_i_phy(layer2) END DO END DO zlevel(nlay_i) = deltaz_i_phy(nlay_i) / 2.0 !--------------------- ! Thermal diffusivity !--------------------- zkimin = 0.1 DO layer = 1, nlay_i - 1 zdummy1 = t_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * & ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) / & ( z_i_phy(layer+1) - z_i_phy(layer) ) - tpw zdummy2 = s_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * & ( s_i_b(ji,layer+1) - s_i_b(ji,layer) ) / & ( z_i_phy(layer+1) - z_i_phy(layer) ) zcp = cpg + lfus * tmut * zdummy2 / & MAX( zdummy1 * zdummy1 , zeps ) zk = xkg + betak1 * zdummy2 / & MIN( -zeps , zdummy1 ) - betak2 * zdummy1 zk = MAX( zk, zkimin ) zthdiff(layer) = zk / ( rhog * zcp ) END DO zcp = cpg + lfus * tmut * s_i_b(ji,nlay_i) / & MAX( ( t_i_b(ji,nlay_i) - tpw ) * ( t_i_b(ji,nlay_i) - tpw ), & zeps ) zk = xkg + betak1 * s_i_b(ji,nlay_i) / & MIN( -zeps , t_i_b(ji,nlay_i) - tpw ) & - betak2 * ( t_i_b(ji,nlay_i) - tpw ) zk = MAX( zk, zkimin ) zthdiff(nlay_i) = zk / ( rhog * zcp ) !----------------- ! Rayleigh number !----------------- DO layer = 1, nlay_i rayleigh(layer) = gpes * MAX(zrhodiff(layer),0.0) * & zpermin(layer) * zlevel(layer) / & ( zthdiff(layer) * visc_br ) END DO !------------------- ! 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,*) ' zsigma : ', ( zsigma(layer), & layer = 1, nlay_i) WRITE(numout,*) ' zrhodiff : ', ( zrhodiff(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' zpermin : ', ( zpermin(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' zthdiff : ', ( zthdiff(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' zlevel : ', ( zlevel(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' rayleigh : ', ( rayleigh(layer), & layer = 1, nlay_i ) WRITE(numout,*) ' diff_br : ', ( diff_br(layer), & layer = 1, nlay_i ) WRITE(numout,*) ENDIF 20 CONTINUE ! !------------------------------------------------------------------------------| ! 3) Flooding and flushing velocities !------------------------------------------------------------------------------| ! IF ( ln_write ) THEN WRITE(numout,*) ' - Flooding and flushing velocities ' WRITE(numout,*) ' ' ENDIF DO 30 ji = kideb, kiut !----------------------- ! Permeability switches !----------------------- ! Permeability switch = 1 if brine volume fraction > e_thr_flu zswitch_per = 1.0 zbvmin = 1.0 DO layer = 1, nlay_i zbvmin = MIN( e_i_b(layer) , zbvmin ) ! minimum brine volume END DO IF ( zbvmin .LT. e_thr_flu ) zswitch_per = 0.0 ! Summer switch = 1 if Tsu ge tpw and min brine volume superior than e_thr_flu zswitchs = MAX( 0.0, SIGN ( 1.0d0, t_su_b(ji) - tpw ) ) ! 0 if winter, 1 if summer zswitchs = zswitchs * zswitch_per !------------------- ! Flooding velocity !------------------- zdhitot = dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) zdhstot = dh_s_melt(ji) w_flood = ( ( rho0 - rhog ) / rho0 * zdhitot - & rhon / rho0 * zdhstot ) / ddtb * zswitch_per IF ( w_flood .GT. 0 ) THEN z_flood = 0.0 ELSE z_flood = 1.0 ENDIF IF ( .NOT. ln_flo ) w_flood = 0. !------------------ ! Percolating flux !------------------ ! Percolating flow ( rho dh * beta * switch / rho0 ) qsummer = ( - rhog * MIN ( dh_i_surf(ji) , 0.0 ) & - rhon * MIN ( dh_s_tot(ji) , 0.0 ) ) qsummer = qsummer * flu_beta * zswitchs / 1000.0 ! 1000 is a ref density for brine w_flush = qsummer / ddtb / e_i_b(1) IF ( .NOT. ln_flu ) THEN w_flush = 0. qsummer = 0. ENDIF IF ( ln_write ) THEN WRITE(numout,*) ' zswitchs : ', zswitchs WRITE(numout,*) ' w_flush : ', w_flush WRITE(numout,*) ' w_flood : ', w_flood ENDIF 30 CONTINUE ! !------------------------------------------------------------------------------| ! 4) Compute dummy factors for tracer diffusion equation !------------------------------------------------------------------------------| ! IF ( ln_write ) THEN WRITE(numout,*) ' - Compute dummy factors for tracer diffusion' WRITE(numout,*) ' ' ENDIF DO 40 ji = kideb, kiut !---------------- ! za factors !---------------- DO layer = 1, nlay_i za(layer) = ddtb / ( deltaz_i_phy(layer) * e_i_b(layer) ) END DO !-------------------- ! zb, zc, ze factors !-------------------- DO layer = 1, nlay_i - 1 ! interpolate brine volume at the interface between layers zdummy1 = ( e_i_b(layer + 1 ) - e_i_b(layer) ) / & ( z_i_phy(layer + 1) - z_i_phy(layer) ) zdummy2 = deltaz_i_phy(layer) / 2.0 zdummy3 = e_i_b(layer) + zdummy1 * zdummy2 zb(layer) = zdummy3 * diff_br(layer) / & ( z_i_phy(layer + 1) - z_i_phy(layer) ) zc(layer) = w_flood * zdummy3 * z_flood ze(layer) = ( w_flood * ( 1. - z_flood ) + w_flush ) * zdummy3 ! old qsummer scheme ! ze(layer) = ( w_flood * ( 1. - z_flood ) ) * zdummy3 + ! & zswitchs * qsummer / ddtb END DO ! Fixed boundary condition (imposed cc.) zb(nlay_i) = 2. * e_i_b(nlay_i) / & deltaz_i_phy(nlay_i) * diff_br(nlay_i) zc(nlay_i) = w_flood * e_i_b(nlay_i) * z_flood ze(nlay_i) = ( w_flood * ( 1. - z_flood ) + w_flush ) * & e_i_b(nlay_i) ! ze(nlay_i) = ( w_flood * ( 1. - z_flood ) ) * e_i_b(nlay_i) + ! & zswitchs * qsummer / ddtb IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' -Dummy factors ' WRITE(numout,*) ' za : ', ( za (layer), & layer = 1, nlay_i) WRITE(numout,*) ' zb : ', ( zb (layer), & layer = 1, nlay_i) WRITE(numout,*) ' zc : ', ( zc (layer), & layer = 1, nlay_i) WRITE(numout,*) ' ze : ', ( ze (layer), & layer = 1, nlay_i) ENDIF 40 CONTINUE ! !----------------------------------------------------------------------- ! 5) Tridiagonal system terms for tracer diffusion equation !----------------------------------------------------------------------- ! DO 50 ji = kideb, kiut !---------------- ! first equation !---------------- ztrid(1,1) = 0.0 ztrid(1,2) = 1.0 + za(1) * ( zb(1) + ze(1) ) ztrid(1,3) = za(1) * ( -zb(1) + zc(1) ) zind(1) = z_sbr_i(1) !----------------- ! inner equations !----------------- DO layer = 2, nlay_i - 1 ztrid(layer,1) = - za(layer) * ( zb(layer-1) + ze(layer-1) ) ztrid(layer,2) = 1.0 + za(layer) * ( zb(layer) + & ze(layer) + zb(layer-1) - zc(layer-1) ) ztrid(layer,3) = za(layer) * ( -zb(layer) + zc(layer) ) zind(layer) = z_sbr_i(layer) END DO !---------------- ! last equation !---------------- WRITE(numout,*) " nlay_i : ", nlay_i ztrid(nlay_i,1) = - za(nlay_i) * ( zb(nlay_i-1) + ze(nlay_i-1) ) ztrid(nlay_i,2) = 1.0 + za(nlay_i) * ( zb(nlay_i) + & ze(nlay_i) + zb(nlay_i-1) - zc(nlay_i-1) ) ztrid(nlay_i,3) = 0. zind(nlay_i) = z_sbr_i(nlay_i) + za(nlay_i) * ( zb(nlay_i) & - zc(nlay_i) ) * oce_sal IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' -Tridiag terms, winter ... ' WRITE(numout,*) DO layer = 1, nlay_i WRITE(numout,*) ' layer : ', layer WRITE(numout,*) ' ztrid : ', ztrid(layer,1), & ztrid(layer,2), ztrid(layer,3) WRITE(numout,*) ' zind : ', zind(layer) END DO ENDIF 50 CONTINUE ! !----------------------------------------------------------------------- ! 6) Solving the tridiagonal system !----------------------------------------------------------------------- ! DO 60 ji = kideb, kiut ! The tridiagonal system is solved with Gauss elimination ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, ! McGraw-Hill 1984. zindtbis(1) = zind(1) zdiagbis(1) = ztrid(1,2) DO layer = 2, nlay_i zdiagbis(layer) = ztrid(layer,2) - ztrid(layer,1) * & ztrid(layer-1,3) / zdiagbis(layer-1) zindtbis(layer) = zind(layer) - ztrid(layer,1) * & zindtbis(layer-1) / zdiagbis(layer-1) END DO ! Recover brine salinity z_sbr_i(nlay_i) = zindtbis(nlay_i) / zdiagbis(nlay_i) DO layer = nlay_i - 1 , 1 , -1 z_sbr_i(layer) = ( zindtbis(layer) - ztrid(layer,3)* & z_sbr_i(layer+1)) / zdiagbis(layer) END DO ! Recover ice salinity DO layer = 1, nlay_i sn_i_b(layer) = z_sbr_i(layer) * e_i_b(layer) ! s_i_b(ji,layer) = z_sbr_i(layer) * e_i_b(layer) END DO IF ( ln_write ) THEN WRITE(numout,*) WRITE(numout,*) ' -Solving the tridiagonal system ... ' WRITE(numout,*) WRITE(numout,*) ' zdiagbis: ', ( zdiagbis(layer) , & layer = 1, nlay_i ) WRITE(numout,*) ' zindtbis: ', ( zdiagbis(layer) , & layer = 1, nlay_i ) WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i(layer) , & layer = 1, nlay_i ) WRITE(numout,*) ' sn_i_b : ', ( sn_i_b(layer) , & layer = 1, nlay_i ) ENDIF 60 CONTINUE ! !----------------------------------------------------------------------- ! 7) Mass of salt conserved ? !----------------------------------------------------------------------- ! DO 70 ji = kideb, kiut ! 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_diff : ',zerror, & z_ms_i_ini,z_ms_i_fin, & zfb , zfsu , ddtb) 70 CONTINUE ENDIF ! ln_sal ! !------------------------------------------------------------------------------| ! End of la sous-routine WRITE(numout,*) END SUBROUTINE