MODULE limthd_ent #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' LIM3 sea-ice model !!---------------------------------------------------------------------- !!====================================================================== !! *** MODULE limthd_ent *** !! Redistribution of Enthalpy in the ice !! on the new vertical grid !! after vertical growth/decay !!====================================================================== !! lim_thd_ent : ice redistribution of enthalpy !! * Modules used USE par_oce ! ocean parameters USE dom_oce USE domain USE in_out_manager USE phycst USE ice_oce ! ice variables USE thd_ice USE iceini USE limistate USE ice USE limvar USE par_ice IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC lim_thd_ent ! called by lim_thd !! * Module variables REAL(wp) :: & ! constant values epsi20 = 1.e-20 , & epsi13 = 1.e-13 , & zzero = 0.e0 , & zone = 1.e0 , & epsi10 = 1.0e-10 !!---------------------------------------------------------------------- !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) !!---------------------------------------------------------------------- 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 !! !! ** Arguments !! !! ** Inputs / Outputs !! !! ** External !! !! ** References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 !! !! ** History : (05-2003) Martin V. UCL-Astr !! (07-2005) Martin for 3d adapatation !! (11-2006) Vectorized by Xavier Fettweis (ASTR) !! (03-2008) Energy conservation and clean code !! * Arguments INTEGER , INTENT(IN):: & kideb , & ! start point on which the the computation is applied kiut , & ! end point on which the the computation is applied jl ! thickness category number INTEGER :: & ji,jk , & ! dummy loop indices zji, zjj , & ! 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 INTEGER, DIMENSION(jpij) :: & snswi , & ! snow switch nbot0 , & ! old layer bottom index icsuind , & ! ice surface index icsuswi , & ! ice surface switch icboind , & ! ice bottom index icboswi , & ! ice bottom switch snicind , & ! snow ice index snicswi , & ! snow ice switch snind ! snow index REAL(wp) :: & zeps, zeps6 , & ! numerical constant very small 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 REAL(wp), DIMENSION(jpij) :: & zh_i , & ! thickness of an ice layer zh_s , & ! thickness of a snow layer zqsnow , & ! enthalpy of the snow put in snow ice zdeltah ! temporary variable REAL(wp), DIMENSION(jpij,0:jkmax+3) :: & zm0 , & ! old layer-system vertical cotes qm0 , & ! old layer-system heat content z_s , & ! new snow system vertical cotes z_i , & ! new ice system vertical cotes zthick0 ! old ice thickness REAL(wp), DIMENSION(jpij,0:jkmax+3) :: & zhl0 ! old and new layer thicknesses REAL(wp), DIMENSION(0:jkmax+3,0:jkmax+3) :: & zrl01 ! Energy conservation REAL(wp), DIMENSION(jpij) :: & zqti_in, zqts_in, & zqti_fin, zqts_fin !------------------------------------------------------------------------------| zeps = 1.0d-20 zeps6 = 1.0d-06 zthick0(:,:) = 0.0 zm0(:,:) = 0.0 qm0(:,:) = 0.0 zrl01(:,:) = 0.0 zhl0(:,:) = 0.0 z_i(:,:) = 0.0 z_s(:,:) = 0.0 ! !------------------------------------------------------------------------------| ! 1) Grid | !------------------------------------------------------------------------------| ! nlays0 = nlay_s nlayi0 = nlay_i DO ji = kideb, kiut zh_i(ji) = old_ht_i_b(ji) / nlay_i zh_s(ji) = old_ht_s_b(ji) / nlay_s ENDDO ! !------------------------------------------------------------------------------| ! 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.0 ENDDO !ji DO jk = 1, nlays0 DO ji = kideb, kiut snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-zeps))) & + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-zeps)))) zdeltah(ji)= zdeltah(ji) + zh_s(ji) END DO ! ji ENDDO ! jk ! snswi(ji) : switch which value equals 1 if snow melts ! 0 if not DO ji = kideb, kiut snswi(ji) = MAX(0,INT(-dh_s_tot(ji)/MAX(zeps,ABS(dh_s_tot(ji))))) ENDDO ! 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.0 ENDDO !ji DO jk = 1, nlayi0 DO ji = kideb, kiut icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-zeps))) & + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-zeps)))) 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,INT(-dh_i_surf(ji)/MAX(zeps , 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.0 ENDDO DO jk = nlayi0, 1, -1 DO ji = kideb, kiut icboind(ji) = (nlayi0+1-jk) & * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-zeps))) & + icboind(ji) & * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-zeps)))) zdeltah(ji) = zdeltah(ji) + zh_i(ji) END DO ENDDO DO ji = kideb, kiut ! case of total ablation with remaining snow IF ( ( ht_i_b(ji) .GT. zeps ) .AND. & ( ht_i_b(ji) - dh_snowice(ji) .LT. zeps ) ) 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,INT(dh_i_bott(ji) / MAX(zeps,ABS(dh_i_bott(ji))))) ENDDO ! 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.0 ENDDO DO jk = nlays0, 1, -1 DO ji = kideb, kiut snicind(ji) = (nlays0+1-jk) & * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-zeps))) & + snicind(ji) & * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-zeps)))) zdeltah(ji) = zdeltah(ji) + zh_s(ji) END DO ENDDO ! snicswi(ji) : switch which equals ! 1 if snow-ice forms ! 0 if not DO ji = kideb, kiut snicswi(ji) = MAX(0,INT(dh_snowice(ji)/MAX(zeps,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.0 maxnbot0 = MAX ( maxnbot0 , nbot0(ji) ) ENDDO 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) * limsum END DO ENDDO DO ji = kideb, kiut zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + & zh_s(ji) * nlays0 zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + & snswi(ji) * zm0(ji,1) ENDDO DO jk = ntop0, maxnbot0 DO ji = kideb, kiut ! layer thickness zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) END DO ENDDO zqts_in(:) = 0.0 DO ji = kideb, kiut ! layer heat content qm0(ji,1) = rhosn * ( cpic * ( rtt - ( 1. - snswi(ji) ) * ( tatm_ice_1d(ji) ) & - snswi(ji) * t_s_b(ji,1) ) & + lfus ) * zthick0(ji,1) zqts_in(ji) = zqts_in(ji) + qm0(ji,1) ENDDO 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, zeps - ht_s_b(ji) ) ) zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch END DO ! jk ENDDO ! 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.0 ENDDO DO jk = nlays0, 1, -1 DO ji = kideb, kiut zhsnow = MAX(0.0,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 ENDDO 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) / nlay_s z_s(ji,0) = 0.0 ENDDO DO jk = 1, nlay_s DO ji = kideb, kiut z_s(ji,jk) = zh_s(ji) * jk END DO ENDDO !----------------- ! Layer thickness !----------------- DO layer0 = ntop0, maxnbot0 DO ji = kideb, kiut zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) END DO ENDDO DO layer1 = ntop1, nbot1 DO ji = kideb, kiut q_s_b(ji,layer1)= 0.0 END DO ENDDO !---------------- ! Weight factors !---------------- DO layer0 = ntop0, maxnbot0 DO layer1 = ntop1, nbot1 DO ji = kideb, kiut zrl01(layer1,layer0) = 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,nbot0(ji)-layer0+zeps)) END DO END DO ENDDO ! Heat conservation zqts_fin(:) = 0.0 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 ) THEN DO ji = kideb, kiut IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 WRITE(numout,*) ' violation of heat conservation : ', & ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice WRITE(numout,*) ' ji, jj : ', zji, zjj WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) WRITE(numout,*) ' zqts_in : ', zqts_in(ji) / rdt_ice WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) / rdt_ice 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_i DO ji = kideb, kiut q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , zeps ) END DO !ji ENDDO !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, zeps - ht_s_b(ji) ) ) t_s_b(ji,jk) = rtt & + ( 1.0 - zswitch ) * & ( - zfac1 * q_s_b(ji,jk) + zfac2 ) END DO ENDDO ! !------------------------------------------------------------------------------| ! 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) ) ENDDO !------------------------- ! Cotes of old ice layers !------------------------- zm0(:,0) = 0.0 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)= icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) & + limsum * zh_i(ji) END DO ENDDO DO ji = kideb, kiut zm0(ji,nbot0(ji)) = icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) + dh_i_bott(ji) & + zh_i(ji) * nlayi0 zm0(ji,1) = snicswi(ji)*dh_snowice(ji) + (1-snicswi(ji))*zm0(ji,1) ENDDO !----------------------------- ! 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 ENDDO !--------------------------- ! 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),-zeps) ) - rcp*(ztmelts-rtt) ) & * zthick0(ji,jk) END DO ENDDO !---------------------------- ! Bottom layers heat content !---------------------------- DO ji = kideb, kiut ztmelts = ( 1.0 - icboswi(ji) ) * (-tmut * s_i_b(ji,nlayi0)) & ! case of melting ice + icboswi(ji) * (-tmut * s_i_new(ji)) & ! case of forming ice + rtt ! this temperature is in Celsius ! bottom formation temperature ztform = t_i_b(ji,nlay_i) IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji) qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice + icboswi(ji) * & ! case of forming ice rhoic*( cpic*(ztmelts-ztform) & + lfus *( 1.0-(ztmelts-rtt)/ & MIN ( (ztform-rtt) , - epsi10 ) ) & - rcp*(ztmelts-rtt) ) & *zthick0(ji,nbot0(ji)) ENDDO !----------------------------- ! 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 * 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) = snicswi(ji) * zqsnic + ( 1 - snicswi(ji) ) * qm0(ji,1) ENDDO ! 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)-zeps6+zeps) ) & * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + zeps ) ) END DO ENDDO !------------- ! NEW PROFILE !------------- !--------------- ! Vectors index !--------------- ntop1 = 1 nbot1 = nlay_i !------------------ ! Layers thickness !------------------ DO ji = kideb, kiut zh_i(ji) = ht_i_b(ji) / nlay_i ENDDO !------------- ! Layer cotes !------------- z_i(:,0) = 0.0 DO jk = 1, nlay_i DO ji = kideb, kiut z_i(ji,jk) = zh_i(ji) * jk END DO ENDDO !--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 ENDDO !------------------------ ! Weights for relayering !------------------------ q_i_b(:,:) = 0.0 DO layer0 = ntop0, maxnbot0 DO layer1 = ntop1, nbot1 DO ji = kideb, kiut zrl01(layer1,layer0) = 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)-zeps6+zeps)) & * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+zeps)) END DO END DO ENDDO !------------------------- ! Heat conservation check !------------------------- zqti_fin(:) = 0.0 DO jk = 1, nlay_i DO ji = kideb, kiut zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) END DO END DO ! DO ji = kideb, kiut IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN zji = MOD( npb(ji) - 1, jpi ) + 1 zjj = ( npb(ji) - 1 ) / jpi + 1 WRITE(numout,*) ' violation of heat conservation : ', & ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice WRITE(numout,*) ' ji, jj : ', zji, zjj WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) WRITE(numout,*) ' zqti_in : ', zqti_in(ji) / rdt_ice WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) / rdt_ice 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 !---------------------- ! Recover heat content !---------------------- DO jk = 1, nlay_i DO ji = kideb, kiut q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , zeps ) END DO !ji ENDDO !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 END SUBROUTINE lim_thd_ent #else !!====================================================================== !! *** MODULE limthd_ent *** !! no sea ice model !!====================================================================== CONTAINS SUBROUTINE lim_thd_ent ! Empty routine END SUBROUTINE lim_thd_ent #endif END MODULE limthd_ent