MODULE icevar !!====================================================================== !! *** MODULE icevar *** !! Different sets of ice model variables !! how to switch from one to another !! !! There are three sets of variables !! VGLO : global variables of the model !! - v_i (jpi,jpj,jpl) !! - v_s (jpi,jpj,jpl) !! - a_i (jpi,jpj,jpl) !! - t_s (jpi,jpj,jpl) !! - e_i (jpi,jpj,nlay_i,jpl) !! - smv_i(jpi,jpj,jpl) !! - oa_i (jpi,jpj,jpl) !! VEQV : equivalent variables sometimes used in the model !! - ht_i(jpi,jpj,jpl) !! - ht_s(jpi,jpj,jpl) !! - t_i (jpi,jpj,nlay_i,jpl) !! ... !! VAGG : aggregate variables, averaged/summed over all !! thickness categories !! - vt_i(jpi,jpj) !! - vt_s(jpi,jpj) !! - at_i(jpi,jpj) !! - et_s(jpi,jpj) !total snow heat content !! - et_i(jpi,jpj) !total ice thermal content !! - smt_i(jpi,jpj) !mean ice salinity !! - tm_i (jpi,jpj) !mean ice temperature !!====================================================================== !! History : - ! 2006-01 (M. Vancoppenolle) Original code !! 3.4 ! 2011-02 (G. Madec) dynamical allocation !! 3.5 ! 2012 (M. Vancoppenolle) add ice_var_itd !! 3.6 ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM3 sea-ice model !!---------------------------------------------------------------------- !! ice_var_agg : integrate variables over layers and categories !! ice_var_glo2eqv : transform from VGLO to VEQV !! ice_var_eqv2glo : transform from VEQV to VGLO !! ice_var_salprof : salinity profile in the ice !! ice_var_bv : brine volume !! ice_var_salprof1d : salinity profile in the ice 1D !! ice_var_zapsmall : remove very small area and volume !! ice_var_itd : convert 1-cat to multiple cat !!---------------------------------------------------------------------- USE par_oce ! ocean parameters USE phycst ! physical constants (ocean directory) USE sbc_oce , ONLY : sss_m USE ice ! ice variables USE ice1D ! ice variables (thermodynamics) ! USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) IMPLICIT NONE PRIVATE PUBLIC ice_var_agg PUBLIC ice_var_glo2eqv PUBLIC ice_var_eqv2glo PUBLIC ice_var_salprof PUBLIC ice_var_bv PUBLIC ice_var_salprof1d PUBLIC ice_var_zapsmall PUBLIC ice_var_itd !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2017) !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_var_agg( kn ) !!------------------------------------------------------------------ !! *** ROUTINE ice_var_agg *** !! !! ** Purpose : aggregates ice-thickness-category variables to !! all-ice variables, i.e. it turns VGLO into VAGG !!------------------------------------------------------------------ INTEGER, INTENT( in ) :: kn ! =1 state variables only ! ! >1 state variables + others ! INTEGER :: ji, jj, jk, jl ! dummy loop indices REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i, z1_vt_i !!------------------------------------------------------------------ ! ! ! integrated values vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 ) vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 ) at_i(:,:) = SUM( a_i(:,:,:) , dim=3 ) et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) ! MV MP 2016 IF ( ln_pnd ) THEN ! Melt pond at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ENDIF ! END MP 2016 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction IF( kn > 1 ) THEN ! ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) ) WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:) ELSEWHERE ; z1_at_i(:,:) = 0._wp END WHERE WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:) ELSEWHERE ; z1_vt_i(:,:) = 0._wp END WHERE ! ! ! mean ice/snow thickness htm_i(:,:) = vt_i(:,:) * z1_at_i(:,:) htm_s(:,:) = vt_s(:,:) * z1_at_i(:,:) ! ! ! mean temperature (K), salinity and age tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) ! tm_i (:,:) = 0._wp smt_i(:,:) = 0._wp DO jl = 1, jpl DO jk = 1, nlay_i tm_i (:,:) = tm_i (:,:) + r1_nlay_i * t_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:) smt_i(:,:) = smt_i(:,:) + r1_nlay_i * s_i(:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:) END DO END DO ! !!gm QUESTION 1 : why salinity is named smt_i and not just sm_i ? since the 4D field is named s_i. (NB for temp: tm_i, t_i) ! DEALLOCATE( z1_at_i , z1_vt_i ) ENDIF ! END SUBROUTINE ice_var_agg SUBROUTINE ice_var_glo2eqv !!------------------------------------------------------------------ !! *** ROUTINE ice_var_glo2eqv *** !! !! ** Purpose : computes equivalent variables as function of !! global variables, i.e. it turns VGLO into VEQV !!------------------------------------------------------------------ INTEGER :: ji, jj, jk, jl ! dummy loop indices REAL(wp) :: ze_i, z1_cp, z1_2cp ! local scalars REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - - REAL(wp) :: zhmax, z1_zhmax, zsm_i ! - - REAL(wp) :: zlay_i, zlay_s ! - - REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i !!------------------------------------------------------------------ !!gm Question 2: It is possible to define existence of sea-ice in a common way between !! ice area and ice volume ? !! the idea is to be able to define one for all at the begining of this routine !! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 ) !------------------------------------------------------- ! Ice thickness, snow thickness, ice salinity, ice age !------------------------------------------------------- ! !--- inverse of the ice area WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) ELSEWHERE ; z1_a_i(:,:,:) = 0._wp END WHERE ! ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness zhmax = hi_max(jpl) z1_zhmax = 1._wp / hi_max(jpl) WHERE( ht_i(:,:,jpl) > zhmax ) !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area ht_i (:,:,jpl) = zhmax a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max END WHERE ht_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) !--- snow thickness o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) !--- ice age IF( nn_icesal == 2 ) THEN !--- salinity (with a minimum value imposed everywhere) WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) ) ELSEWHERE ; sm_i(:,:,:) = rn_simin END WHERE ENDIF CALL ice_var_salprof ! salinity profile !------------------- ! Ice temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere) !------------------- zlay_i = REAL( nlay_i , wp ) ! number of layers z1_2cp = 1._wp / ( 2._wp * cpic ) DO jl = 1, jpl DO jk = 1, nlay_i DO jj = 1, jpj DO ji = 1, jpi IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area ! ze_i = e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] ! Conversion q(S,T) -> T (second order equation) zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts ! ELSE !--- no ice t_i(ji,jj,jk,jl) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) !!clem: I think we should impose rt0 instead ENDIF END DO END DO END DO END DO !-------------------- ! Snow temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere) !-------------------- zlay_s = REAL( nlay_s , wp ) z1_cp = 1._wp / cpic DO jk = 1, nlay_s WHERE( v_s(:,:,:) > epsi20 ) !--- icy area t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0 ELSEWHERE !--- no ice t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) END WHERE END DO ! integrated values vt_i (:,:) = SUM( v_i, dim=3 ) vt_s (:,:) = SUM( v_s, dim=3 ) at_i (:,:) = SUM( a_i, dim=3 ) ! MV MP 2016 ! probably should resum for melt ponds ??? ! END SUBROUTINE ice_var_glo2eqv SUBROUTINE ice_var_eqv2glo !!------------------------------------------------------------------ !! *** ROUTINE ice_var_eqv2glo *** !! !! ** Purpose : computes global variables as function of !! equivalent variables, i.e. it turns VEQV into VGLO !!------------------------------------------------------------------ ! v_i (:,:,:) = ht_i(:,:,:) * a_i(:,:,:) v_s (:,:,:) = ht_s(:,:,:) * a_i(:,:,:) smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) ! END SUBROUTINE ice_var_eqv2glo SUBROUTINE ice_var_salprof !!------------------------------------------------------------------ !! *** ROUTINE ice_var_salprof *** !! !! ** Purpose : computes salinity profile in function of bulk salinity !! !! ** Method : If bulk salinity greater than zsi1, !! the profile is assumed to be constant (S_inf) !! If bulk salinity lower than zsi0, !! the profile is linear with 0 at the surface (S_zero) !! If it is between zsi0 and zsi1, it is a !! alpha-weighted linear combination of s_inf and s_zero !! !! ** References : Vancoppenolle et al., 2007 !!------------------------------------------------------------------ INTEGER :: ji, jj, jk, jl ! dummy loop index REAL(wp) :: zsal, z1_dS REAL(wp) :: zargtemp , zs0, zs REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only REAL(wp), PARAMETER :: zsi0 = 3.5_wp REAL(wp), PARAMETER :: zsi1 = 4.5_wp !!------------------------------------------------------------------ !!gm Question: Remove the option 3 ? How many years since it last use ? SELECT CASE ( nn_icesal ) ! ! !---------------------------------------! CASE( 1 ) ! constant salinity in time and space ! ! !---------------------------------------! s_i (:,:,:,:) = rn_icesal sm_i(:,:,:) = rn_icesal ! ! !---------------------------------------------! CASE( 2 ) ! time varying salinity with linear profile ! ! !---------------------------------------------! ! ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) ) ! DO jk = 1, nlay_i s_i(:,:,jk,:) = sm_i(:,:,:) END DO ! ! Slope of the linear profile WHERE( ht_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * sm_i(:,:,:) / ht_i(:,:,:) ELSEWHERE ; z_slope_s(:,:,:) = 0._wp END WHERE ! z1_dS = 1._wp / ( zsi1 - zsi0 ) DO jl = 1, jpl DO jj = 1, jpj DO ji = 1, jpi zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - sm_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) ! ! force a constant profile when SSS too low (Baltic Sea) IF( 2._wp * sm_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp END DO END DO END DO ! Computation of the profile DO jl = 1, jpl DO jk = 1, nlay_i DO jj = 1, jpj DO ji = 1, jpi ! ! linear profile with 0 surface value zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) ! weighting the profile s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) END DO END DO END DO END DO ! DEALLOCATE( z_slope_s , zalpha ) ! ! !-------------------------------------------! CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile ! !-------------------------------------------! (mean = 2.30) ! sm_i(:,:,:) = 2.30_wp !!gm Remark: if we keep the case 3, then compute an store one for all time-step !! a array S_prof(1:nlay_i) containing the calculation and just do: ! DO jk = 1, nlay_i ! s_i(:,:,jk,:) = S_prof(jk) ! END DO !!gm end ! DO jl = 1, jpl DO jk = 1, nlay_i zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i s_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) END DO END DO ! END SELECT ! END SUBROUTINE ice_var_salprof SUBROUTINE ice_var_bv !!------------------------------------------------------------------ !! *** ROUTINE ice_var_bv *** !! !! ** Purpose : computes mean brine volume (%) in sea ice !! !! ** Method : e = - 0.054 * S (ppt) / T (C) !! !! References : Vancoppenolle et al., JGR, 2007 !!------------------------------------------------------------------ INTEGER :: ji, jj, jk, jl ! dummy loop indices !!------------------------------------------------------------------ ! !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done !! instead of setting everything to zero as just below bv_i (:,:,:) = 0._wp DO jl = 1, jpl DO jk = 1, nlay_i WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) END WHERE END DO END DO WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) ELSEWHERE ; bvm_i(:,:) = 0._wp END WHERE ! END SUBROUTINE ice_var_bv SUBROUTINE ice_var_salprof1d !!------------------------------------------------------------------- !! *** ROUTINE ice_var_salprof1d *** !! !! ** Purpose : 1d computation of the sea ice salinity profile !! Works with 1d vectors and is used by thermodynamic modules !!------------------------------------------------------------------- INTEGER :: ji, jk ! dummy loop indices REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars REAL(wp) :: zalpha, zs, zs0 ! - - ! REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s ! REAL(wp), PARAMETER :: zsi0 = 3.5_wp REAL(wp), PARAMETER :: zsi1 = 4.5_wp !!--------------------------------------------------------------------- ! SELECT CASE ( nn_icesal ) ! ! !---------------------------------------! CASE( 1 ) ! constant salinity in time and space ! ! !---------------------------------------! s_i_1d(:,:) = rn_icesal ! ! !---------------------------------------------! CASE( 2 ) ! time varying salinity with linear profile ! ! !---------------------------------------------! ! ALLOCATE( z_slope_s(jpij) ) ! DO ji = 1, nidx ! Slope of the linear profile rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) END DO z1_dS = 1._wp / ( zsi1 - zsi0 ) DO jk = 1, nlay_i DO ji = 1, nidx zalpha = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) ) IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha = 0._wp ! cst profile when SSS too low (Baltic Sea) ! zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i ! linear profile with 0 surface value zs = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji) ! weighting the profile s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) ! bounding salinity END DO END DO ! DEALLOCATE( z_slope_s ) ! !-------------------------------------------! CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile ! !-------------------------------------------! (mean = 2.30) ! sm_i_1d(:) = 2.30_wp ! !!gm cf remark in ice_var_salprof routine, CASE( 3 ) DO jk = 1, nlay_i zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) ) DO ji = 1, nidx s_i_1d(ji,jk) = zsal END DO END DO ! END SELECT ! END SUBROUTINE ice_var_salprof1d SUBROUTINE ice_var_zapsmall !!------------------------------------------------------------------- !! *** ROUTINE ice_var_zapsmall *** !! !! ** Purpose : Remove too small sea ice areas and correct fluxes !!------------------------------------------------------------------- INTEGER :: ji, jj, jl, jk ! dummy loop indices REAL(wp) :: zsal, zvi, zvs, zei, zes, zvp !!------------------------------------------------------------------- ! DO jl = 1, jpl !== loop over the categories ==! ! !----------------------------------------------------------------- ! Zap ice energy and use ocean heat to melt ice !----------------------------------------------------------------- DO jk = 1, nlay_i DO jj = 1 , jpj DO ji = 1 , jpi !!gm I think we can do better/faster : to be discussed rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch & & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch zei = e_i(ji,jj,jk,jl) e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch ) ! update exchanges with ocean hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 END DO END DO END DO DO jj = 1 , jpj DO ji = 1 , jpi rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch & & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch zsal = smv_i(ji,jj, jl) zvi = v_i (ji,jj, jl) zvs = v_s (ji,jj, jl) zes = e_s (ji,jj,1,jl) IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) zvp = v_ip (ji,jj ,jl) !----------------------------------------------------------------- ! Zap snow energy !----------------------------------------------------------------- t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch ) e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch !----------------------------------------------------------------- ! zap ice and snow volume, add water and salt to ocean !----------------------------------------------------------------- ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch ! MV MP 2016 IF ( ln_pnd ) THEN a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch IF ( ln_pnd_fw ) wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_ip(ji,jj,jl) - zvp ) * rhofw * r1_rdtice ENDIF ! END MV MP 2016 ! update exchanges with ocean sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 END DO END DO END DO ! to be sure that at_i is the sum of a_i(jl) at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) ! open water = 1 if at_i=0 WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp ! END SUBROUTINE ice_var_zapsmall SUBROUTINE ice_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i ) !!------------------------------------------------------------------ !! *** ROUTINE ice_var_itd *** !! !! ** Purpose : converting 1-cat ice to multiple ice categories !! !! ice thickness distribution follows a gaussian law !! around the concentration of the most likely ice thickness !! (similar as iceistate.F90) !! !! ** Method: Iterative procedure !! !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian !! !! 2) Check whether the distribution conserves area and volume, positivity and !! category boundaries !! !! 3) If not (input ice is too thin), the last category is empty and !! the number of categories is reduced (jpl-1) !! !! 4) Iterate until ok (SUM(itest(:) = 4) !! !! ** Arguments : zhti: 1-cat ice thickness !! zhts: 1-cat snow depth !! zai : 1-cat ice concentration !! !! ** Output : jpl-cat !! !! (Example of application: BDY forcings when input are cell averaged) !!------------------------------------------------------------------- INTEGER :: ji, jk, jl ! dummy loop indices INTEGER :: ijpij, i_fill, jl0 REAL(wp) :: zarg, zV, zconv, zdh, zdv REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables INTEGER , DIMENSION(4) :: itest !!------------------------------------------------------------------- !-------------------------------------------------------------------- ! initialisation of variables !-------------------------------------------------------------------- ijpij = SIZE( zhti , 1 ) zht_i(1:ijpij,1:jpl) = 0._wp zht_s(1:ijpij,1:jpl) = 0._wp za_i (1:ijpij,1:jpl) = 0._wp ! ---------------------------------------- ! distribution over the jpl ice categories ! ---------------------------------------- DO ji = 1, ijpij IF( zhti(ji) > 0._wp ) THEN ! find which category (jl0) the input ice thickness falls into jl0 = jpl DO jl = 1, jpl IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN jl0 = jl CYCLE ENDIF END DO ! initialisation of tests itest(:) = 0 i_fill = jpl + 1 !==================================== DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories ! iteration !==================================== i_fill = i_fill - 1 ! initialisation of ice variables for each try zht_i(ji,1:jpl) = 0._wp za_i (ji,1:jpl) = 0._wp itest(:) = 0 ! *** case very thin ice: fill only category 1 IF ( i_fill == 1 ) THEN zht_i(ji,1) = zhti(ji) za_i (ji,1) = zai (ji) ! *** case ice is thicker: fill categories >1 ELSE ! Fill ice thicknesses in the (i_fill-1) cat by hmean DO jl = 1, i_fill - 1 zht_i(ji,jl) = hi_mean(jl) END DO ! Concentrations in the (i_fill-1) categories za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) DO jl = 1, i_fill - 1 IF ( jl /= jl0 ) THEN zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) ENDIF END DO ! Concentration in the last (i_fill) category za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) ! Ice thickness in the last (i_fill) category zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) ! clem: correction if concentration of upper cat is greater than lower cat ! (it should be a gaussian around jl0 but sometimes it is not) IF ( jl0 /= jpl ) THEN DO jl = jpl, jl0+1, -1 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN zdv = zht_i(ji,jl) * za_i(ji,jl) zht_i(ji,jl ) = 0._wp za_i (ji,jl ) = 0._wp za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) END IF ENDDO ENDIF ENDIF ! case ice is thick or thin !--------------------- ! Compatibility tests !--------------------- ! Test 1: area conservation zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) IF ( zconv < epsi06 ) itest(1) = 1 ! Test 2: volume conservation zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) IF ( zconv < epsi06 ) itest(2) = 1 ! Test 3: thickness of the last category is in-bounds ? IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 4: positivity of ice concentrations itest(4) = 1 DO jl = 1, i_fill IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 END DO ! !============================ END DO ! end iteration on categories ! !============================ ENDIF ! if zhti > 0 END DO ! i loop ! ------------------------------------------------ ! Adding Snow in each category where za_i is not 0 ! ------------------------------------------------ DO jl = 1, jpl DO ji = 1, ijpij IF( za_i(ji,jl) > 0._wp ) THEN zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) ! In case snow load is in excess that would lead to transformation from snow to ice ! Then, transfer the snow excess into the ice (different from icethd_dh) zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) ! recompute ht_i, ht_s avoiding out of bounds values zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn ) ENDIF END DO END DO ! END SUBROUTINE ice_var_itd #else !!---------------------------------------------------------------------- !! Default option Dummy module NO LIM3 sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE icevar