MODULE limthd_ent !!====================================================================== !! *** MODULE limthd_ent *** !! Redistribution of Enthalpy in the ice !! on the new vertical grid !! after vertical growth/decay !!====================================================================== !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D !! ! 2005-07 (M. Vancoppenolle) 3D version !! ! 2006-11 (X. Fettweis) Vectorized !! 3.0 ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code !! 4.0 ! 2011-02 (G. Madec) dynamical allocation !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM3 sea-ice model !!---------------------------------------------------------------------- !! lim_thd_ent : ice redistribution of enthalpy !!---------------------------------------------------------------------- USE par_oce ! ocean parameters USE dom_oce ! domain variables USE domain ! USE phycst ! physical constants USE ice ! LIM variables USE par_ice ! LIM parameters USE thd_ice ! LIM thermodynamics USE limvar ! LIM variables USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE wrk_nemo ! work arrays USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) IMPLICIT NONE PRIVATE PUBLIC lim_thd_ent ! called by lim_thd REAL(wp) :: epsi20 = 1.e-20_wp ! constant values REAL(wp) :: epsi10 = 1.e-10_wp ! REAL(wp) :: zzero = 0._wp ! REAL(wp) :: zone = 1._wp ! !!---------------------------------------------------------------------- !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_thd_ent( kideb, kiut, jl ) !!------------------------------------------------------------------- !! *** ROUTINE lim_thd_ent *** !! !! ** Purpose : !! This routine computes new vertical grids !! in the ice and in the snow, and consistently redistributes !! temperatures in the snow / ice. !! Redistribution is made so as to ensure to energy conservation !! !! !! ** Method : linear conservative remapping !! !! ** Steps : 1) Grid !! 2) Switches !! 3) Snow redistribution !! 4) Ice enthalpy redistribution !! 5) Ice salinity, recover temperature !! !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 !!------------------------------------------------------------------- INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied INTEGER , INTENT(in) :: jl ! Thickness cateogry number INTEGER :: ji,jk ! dummy loop indices INTEGER :: ii, ij , & ! dummy indices ntop0 , & ! old layer top index nbot1 , & ! new layer bottom index ntop1 , & ! new layer top index limsum , & ! temporary loop index nlayi0,nlays0 , & ! old number of layers maxnbot0 , & ! old layer bottom index layer0, layer1 ! old/new layer indexes REAL(wp) :: & ztmelts , & ! ice melting point zqsnic , & ! enthalpy of snow ice layer zhsnow , & ! temporary snow thickness variable zswitch , & ! dummy switch argument zfac1 , & ! dummy factor zfac2 , & ! dummy factor ztform , & !: bottom formation temperature zaaa , & !: dummy factor zbbb , & !: dummy factor zccc , & !: dummy factor zdiscrim !: dummy factor INTEGER, POINTER, DIMENSION(:) :: snswi ! snow switch INTEGER, POINTER, DIMENSION(:) :: nbot0 ! old layer bottom index INTEGER, POINTER, DIMENSION(:) :: icsuind ! ice surface index INTEGER, POINTER, DIMENSION(:) :: icsuswi ! ice surface switch INTEGER, POINTER, DIMENSION(:) :: icboind ! ice bottom index INTEGER, POINTER, DIMENSION(:) :: icboswi ! ice bottom switch INTEGER, POINTER, DIMENSION(:) :: snicind ! snow ice index INTEGER, POINTER, DIMENSION(:) :: snicswi ! snow ice switch INTEGER, POINTER, DIMENSION(:) :: snind ! snow index ! REAL(wp), POINTER, DIMENSION(:) :: zh_i ! thickness of an ice layer REAL(wp), POINTER, DIMENSION(:) :: zh_s ! thickness of a snow layer REAL(wp), POINTER, DIMENSION(:) :: zqsnow ! enthalpy of the snow put in snow ice REAL(wp), POINTER, DIMENSION(:) :: zdeltah ! temporary variable REAL(wp), POINTER, DIMENSION(:) :: zqti_in, zqts_in REAL(wp), POINTER, DIMENSION(:) :: zqti_fin, zqts_fin REAL(wp), POINTER, DIMENSION(:,:) :: zm0 ! old layer-system vertical cotes REAL(wp), POINTER, DIMENSION(:,:) :: qm0 ! old layer-system heat content REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! new snow system vertical cotes REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! new ice system vertical cotes REAL(wp), POINTER, DIMENSION(:,:) :: zthick0 ! old ice thickness REAL(wp), POINTER, DIMENSION(:,:) :: zhl0 ! old and new layer thicknesses REAL(wp), POINTER, DIMENSION(:,:) :: zrl01 REAL(wp) :: zinda !!------------------------------------------------------------------- CALL wrk_alloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind ) ! integer CALL wrk_alloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin ) ! real CALL wrk_alloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 ) CALL wrk_alloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 ) zthick0(:,:) = 0._wp zm0 (:,:) = 0._wp qm0 (:,:) = 0._wp zrl01 (:,:) = 0._wp zhl0 (:,:) = 0._wp z_i (:,:) = 0._wp z_s (:,:) = 0._wp ! !------------------------------------------------------------------------------| ! 1) Grid | !------------------------------------------------------------------------------| nlays0 = nlay_s nlayi0 = nlay_i DO ji = kideb, kiut zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i ) zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s ) END DO ! !------------------------------------------------------------------------------| ! 2) Switches | !------------------------------------------------------------------------------| ! 2.1 snind(ji), snswi(ji) ! snow surface behaviour : computation of snind(ji)-snswi(ji) ! snind(ji) : index which equals ! 0 if snow is accumulating ! 1 if 1st layer is melting ! 2 if 2nd layer is melting ... DO ji = kideb, kiut snind (ji) = 0 zdeltah(ji) = 0._wp ENDDO !ji DO jk = 1, nlays0 DO ji = kideb, kiut snind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) & + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji))))) zdeltah(ji)= zdeltah(ji) + zh_s(ji) END DO ! ji END DO ! jk ! snswi(ji) : switch which value equals 1 if snow melts ! 0 if not DO ji = kideb, kiut snswi(ji) = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) END DO ! ji ! 2.2 icsuind(ji), icsuswi(ji) ! ice surface behaviour : computation of icsuind(ji)-icsuswi(ji) ! icsuind(ji) : index which equals ! 0 if nothing happens at the surface ! 1 if first layer is melting ! 2 if 2nd layer is reached by melt ... DO ji = kideb, kiut icsuind(ji) = 0 zdeltah(ji) = 0._wp END DO !ji DO jk = 1, nlayi0 DO ji = kideb, kiut icsuind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) & + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji))))) zdeltah(ji) = zdeltah(ji) + zh_i(ji) END DO ! ji ENDDO !jk ! icsuswi(ji) : switch which equals ! 1 if ice melts at the surface ! 0 if not DO ji = kideb, kiut icsuswi(ji) = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) ENDDO ! 2.3 icboind(ji), icboswi(ji) ! ice bottom behaviour : computation of icboind(ji)-icboswi(ji) ! icboind(ji) : index which equals ! 0 if accretion is on the way ! 1 if last layer has started to melt ! 2 if penultiem layer is melting ... and so on ! N+1 if all layers melt and that snow transforms into ice DO ji = kideb, kiut icboind(ji) = 0 zdeltah(ji) = 0._wp END DO DO jk = nlayi0, 1, -1 DO ji = kideb, kiut icboind(ji) = (nlayi0+1-jk) * NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) & & + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji))))) zdeltah(ji) = zdeltah(ji) + zh_i(ji) END DO END DO DO ji = kideb, kiut ! case of total ablation with remaining snow IF ( ( ht_i_b(ji) .GT. epsi20 ) .AND. & ( ht_i_b(ji) - dh_snowice(ji) .LT. epsi20 ) ) icboind(ji) = nlay_i + 1 END DO ! icboswi(ji) : switch which equals ! 1 if ice accretion is on the way ! 0 if ablation is on the way DO ji = kideb, kiut icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) END DO ! 2.4 snicind(ji), snicswi(ji) ! snow ice formation : calcul de snicind(ji)-snicswi(ji) ! snicind(ji) : index which equals ! 0 if no snow-ice forms ! 1 if last layer of snow has started to melt ! 2 if penultiem layer ... DO ji = kideb, kiut snicind(ji) = 0 zdeltah(ji) = 0._wp END DO DO jk = nlays0, 1, -1 DO ji = kideb, kiut snicind(ji) = (nlays0+1-jk) & * NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji) & * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji))))) zdeltah(ji) = zdeltah(ji) + zh_s(ji) END DO END DO ! snicswi(ji) : switch which equals ! 1 if snow-ice forms ! 0 if not DO ji = kideb, kiut snicswi(ji) = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) ENDDO ! !------------------------------------------------------------------------------| ! 3) Snow redistribution | !------------------------------------------------------------------------------| ! !------------- ! Old profile !------------- ! by 'old', it is meant that layers coming from accretion are included, ! and that interfacial layers which were partly melted are reduced ! indexes of the vectors !------------------------ ntop0 = 1 maxnbot0 = 0 DO ji = kideb, kiut nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) ! cotes of the top of the layers zm0(ji,0) = 0._wp maxnbot0 = MAX ( maxnbot0 , nbot0(ji) ) END DO IF( lk_mpp ) CALL mpp_max( maxnbot0, kcom=ncomm_ice ) DO jk = 1, maxnbot0 DO ji = kideb, kiut !change limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) limsum = MIN( limsum , nlay_s ) zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * REAL( limsum ) END DO END DO DO ji = kideb, kiut zm0(ji,nbot0(ji)) = dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 ) zm0(ji,1) = dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1) END DO DO jk = ntop0, maxnbot0 DO ji = kideb, kiut zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) ! layer thickness END DO END DO zqts_in(:) = 0._wp DO ji = kideb, kiut ! layer heat content qm0 (ji,1) = rhosn * ( cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji) & & - REAL( snswi(ji) ) * t_s_b (ji,1) ) & & + lfus ) * zthick0(ji,1) zqts_in(ji) = zqts_in(ji) + qm0(ji,1) END DO DO jk = 2, maxnbot0 DO ji = kideb, kiut limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) limsum = MIN( limsum , nlay_s ) qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) ) zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch END DO ! jk END DO ! ji !------------------------------------------------ ! Energy given by the snow in snow-ice formation !------------------------------------------------ ! zqsnow, enthalpy of the flooded snow DO ji = kideb, kiut zqsnow (ji) = rhosn * lfus zdeltah(ji) = 0._wp END DO DO jk = nlays0, 1, -1 DO ji = kideb, kiut zhsnow = MAX( 0._wp , dh_snowice(ji)-zdeltah(ji) ) zqsnow (ji) = zqsnow (ji) + rhosn*cpic*(rtt-t_s_b(ji,jk)) zdeltah(ji) = zdeltah(ji) + zh_s(ji) END DO END DO DO ji = kideb, kiut zqsnow(ji) = zqsnow(ji) * dh_snowice(ji) END DO !------------------ ! new snow profile !------------------ !-------------- ! Vector index !-------------- ntop1 = 1 nbot1 = nlay_s !------------------- ! Layer coordinates !------------------- DO ji = kideb, kiut zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) z_s(ji,0) = 0._wp ENDDO DO jk = 1, nlay_s DO ji = kideb, kiut z_s(ji,jk) = zh_s(ji) * REAL( jk ) END DO END DO !----------------- ! Layer thickness !----------------- DO layer0 = ntop0, maxnbot0 DO ji = kideb, kiut zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) END DO END DO DO layer1 = ntop1, nbot1 DO ji = kideb, kiut q_s_b(ji,layer1) = 0._wp END DO END DO !---------------- ! Weight factors !---------------- DO layer0 = ntop0, maxnbot0 DO layer1 = ntop1, nbot1 DO ji = kideb, kiut zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) & & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & & * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) END DO END DO END DO ! Heat conservation zqts_fin(:) = 0._wp DO jk = 1, nlay_s DO ji = kideb, kiut zqts_fin(ji) = zqts_fin(ji) + q_s_b(ji,jk) END DO END DO IF ( con_i .AND. jiindex_1d > 0 ) THEN DO ji = kideb, kiut IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN ii = MOD( npb(ji) - 1, jpi ) + 1 ij = ( npb(ji) - 1 ) / jpi + 1 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice WRITE(numout,*) ' ji, jj : ', ii, ij WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) WRITE(numout,*) ' zqts_in : ', zqts_in (ji) * r1_rdtice WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) * r1_rdtice WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji) WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) WRITE(numout,*) ' snswi : ', snswi(ji) ENDIF END DO ENDIF !--------------------- ! Recover heat content !--------------------- DO jk = 1, nlay_s DO ji = kideb, kiut zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) ) q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 ) END DO !ji END DO !jk !--------------------- ! Recover temperature !--------------------- zfac1 = 1. / ( rhosn * cpic ) zfac2 = lfus / cpic DO jk = 1, nlay_s DO ji = kideb, kiut zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) END DO END DO ! !------------------------------------------------------------------------------| ! 4) Ice redistribution | !------------------------------------------------------------------------------| ! !------------- ! OLD PROFILE !------------- !---------------- ! Vector indexes !---------------- ntop0 = 1 maxnbot0 = 0 DO ji = kideb, kiut ! reference number of the bottommost layer nbot0(ji) = MAX( 1 , MIN( nlayi0 + ( 1 - icboind(ji) ) + & & ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , nlay_i + 2 ) ) ! maximum reference number of the bottommost layer over all domain maxnbot0 = MAX( maxnbot0 , nbot0(ji) ) END DO !------------------------- ! Cotes of old ice layers !------------------------- zm0(:,0) = 0._wp DO jk = 1, maxnbot0 DO ji = kideb, kiut ! jk goes from 1 to nbot0 ! the ice layer number goes from 1 to nlay_i ! limsum is the real ice layer number corresponding to present jk limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) zm0(ji,jk)= REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) & + REAL(limsum) * zh_i(ji) END DO END DO DO ji = kideb, kiut zm0(ji,nbot0(ji)) = REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) & + zh_i(ji) * REAL(nlayi0) zm0(ji,1) = REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1) END DO !----------------------------- ! Thickness of old ice layers !----------------------------- DO jk = ntop0, maxnbot0 DO ji = kideb, kiut zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) END DO END DO !--------------------------- ! Inner layers heat content !--------------------------- qm0(:,:) = 0.0 zqti_in(:) = 0.0 DO jk = ntop0, maxnbot0 DO ji = kideb, kiut limsum = MAX(1,MIN(snicswi(ji)*(jk-1) + icsuswi(ji)*(jk-1+icsuind(ji)) + & (1-icsuswi(ji))*(1-snicswi(ji))*jk,nlay_i)) ztmelts = -tmut * s_i_b(ji,limsum) + rtt qm0(ji,jk) = rhoic * ( cpic * (ztmelts-t_i_b(ji,limsum)) + lfus * ( 1.0-(ztmelts-rtt)/ & MIN((t_i_b(ji,limsum)-rtt),-epsi20) ) - rcp*(ztmelts-rtt) ) & * zthick0(ji,jk) END DO END DO !---------------------------- ! Bottom layers heat content !---------------------------- DO ji = kideb, kiut ztmelts = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice & + REAL( icboswi(ji) ) * (-tmut * s_i_new(ji) ) & ! case of forming ice & + rtt ! in Kelvin ! bottom formation temperature ztform = t_i_b(ji,nlay_i) IF( num_sal == 2 ) ztform = t_bo_b(ji) qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice & + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) & - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji) ) END DO !----------------------------- ! Snow ice layer heat content !----------------------------- DO ji = kideb, kiut ! energy of the flooding seawater zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive ! Heat conservation diagnostic qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic qldif_1d(ji) = qldif_1d(ji) + zqsnic * a_i_b(ji) ! enthalpy of the newly formed snow-ice layer ! = enthalpy of snow + enthalpy of frozen water zqsnic = zqsnow(ji) + zqsnic qm0(ji,1) = REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1) END DO ! ji DO jk = ntop0, maxnbot0 DO ji = kideb, kiut ! Heat conservation zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) & & * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) END DO END DO !------------- ! NEW PROFILE !------------- !--------------- ! Vectors index !--------------- ntop1 = 1 nbot1 = nlay_i !------------------ ! Layers thickness !------------------ DO ji = kideb, kiut zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) ENDDO !------------- ! Layer cotes !------------- z_i(:,0) = 0._wp DO jk = 1, nlay_i DO ji = kideb, kiut z_i(ji,jk) = zh_i(ji) * jk END DO END DO !--thicknesses of the layers DO layer0 = ntop0, maxnbot0 DO ji = kideb, kiut zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) ! thicknesses of the layers END DO END DO !------------------------ ! Weights for relayering !------------------------ q_i_b(:,:) = 0._wp DO layer0 = ntop0, maxnbot0 DO layer1 = ntop1, nbot1 DO ji = kideb, kiut zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) q_i_b(ji,layer1) = q_i_b(ji,layer1) & + zrl01(layer1,layer0)*qm0(ji,layer0) & * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) & * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) END DO END DO END DO !------------------------- ! Heat conservation check !------------------------- zqti_fin(:) = 0._wp DO jk = 1, nlay_i DO ji = kideb, kiut zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) END DO END DO ! IF ( con_i .AND. jiindex_1d > 0 ) THEN DO ji = kideb, kiut IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN ii = MOD( npb(ji) - 1, jpi ) + 1 ij = ( npb(ji) - 1 ) / jpi + 1 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice WRITE(numout,*) ' ji, jj : ', ii, ij WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) WRITE(numout,*) ' zqti_in : ', zqti_in (ji) * r1_rdtice WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) WRITE(numout,*) ' icsuswi : ', icsuswi(ji) WRITE(numout,*) ' icboswi : ', icboswi(ji) WRITE(numout,*) ' snicswi : ', snicswi(ji) ENDIF END DO ENDIF !---------------------- ! Recover heat content !---------------------- DO jk = 1, nlay_i DO ji = kideb, kiut zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) ) q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 ) END DO !ji END DO !jk ! Heat conservation zqti_fin(:) = 0.0 DO jk = 1, nlay_i DO ji = kideb, kiut zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) * zh_i(ji) END DO END DO ! !------------------------------------------------------------------------------| ! 5) Update salinity and recover temperature | !------------------------------------------------------------------------------| ! ! Update salinity (basal entrapment, snow ice formation) DO ji = kideb, kiut sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) END DO !ji ! Recover temperature DO jk = 1, nlay_i DO ji = kideb, kiut ztmelts = -tmut*s_i_b(ji,jk) + rtt !Conversion q(S,T) -> T (second order equation) zaaa = cpic zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus zccc = lfus * ( ztmelts - rtt ) zdiscrim = SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) ) t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa ) END DO !ji END DO !jk ! CALL wrk_dealloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind ) ! integer CALL wrk_dealloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin ) ! real CALL wrk_dealloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 ) CALL wrk_dealloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 ) ! END SUBROUTINE lim_thd_ent #else !!---------------------------------------------------------------------- !! Default option NO LIM3 sea-ice model !!---------------------------------------------------------------------- CONTAINS SUBROUTINE lim_thd_ent ! Empty routine END SUBROUTINE lim_thd_ent #endif !!====================================================================== END MODULE limthd_ent