source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 4635

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
   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
   !!             -   ! 2014-05 (C. Rousset) complete rewriting
   !!----------------------------------------------------------------------
#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 sbc_oce        ! Surface boundary condition: ocean fields
   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   ! constant values
   REAL(wp) :: epsi10 = 1.e-10   ! constant values

   !!----------------------------------------------------------------------
   !! 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 consistently redistributes temperatures. 
      !!           Redistribution is made so as to ensure to energy conservation
      !!
      !!
      !! ** Method  : linear conservative remapping
      !!           
      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
      !!            2) linear remapping on the new layers
      !!            3) Ice salinity update + recover temperature from enthalpies
      !!
      !! 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,ii,ij   !  dummy loop indices
      INTEGER  :: jk0, jk1   !  old/new layer indices
      REAL(wp) :: ztmelts    ! temperature of melting
      REAL(wp) :: zswitch, zaaa, zbbb, zccc, zdiscrim ! converting enthalpy to temperature
      !
      REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
      REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces
      !!-------------------------------------------------------------------

      CALL wrk_alloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 )
      CALL wrk_alloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 )

      !--------------------------------------------------------------------------
      !  1) Cumulative integral of old enthalpy * thicnkess and layers interfaces
      !--------------------------------------------------------------------------
      zqh_cum0(:,0:nlay_i+2) = 0._wp 
      zh_cum0 (:,0:nlay_i+2) = 0._wp
      DO jk0 = 1, nlay_i+2
         DO ji = kideb, kiut
            zqh_cum0(ji,jk0) = zqh_cum0(ji,jk0-1) + qh_i_old(ji,jk0-1)
            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1)
         ENDDO
      ENDDO

      !------------------------------------
      !  2) Interpolation on the new layers
      !------------------------------------
      ! new layers interfaces
      zh_cum1(:,0:nlay_i) = 0._wp
      DO jk1 = 1, nlay_i
         DO ji = kideb, kiut
            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + ht_i_b(ji) / REAL( nlay_i )
         ENDDO
      ENDDO

      zqh_cum1(:,0:nlay_i) = 0._wp 
      ! new cumulative q*h => linear interpolation
      DO jk0 = 1, nlay_i+1
         DO jk1 = 1, nlay_i-1
            DO ji = kideb, kiut
               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
                  zqh_cum1(ji,jk1) = ( zqh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  &
                     &                 zqh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  &
                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) )
               ENDIF
            ENDDO
         ENDDO
      ENDDO
      ! to ensure that total heat content is strictly conserved, set:
      zqh_cum1(:,nlay_i) = zqh_cum0(:,nlay_i+2) 

      ! new enthalpies
      DO jk1 = 1, nlay_i
         DO ji = kideb, kiut
            zswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + epsi20 ) ) 
            q_i_b(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) * REAL( nlay_i ) / MAX( ht_i_b(ji), epsi20 )
         ENDDO
      ENDDO

      !---------------------------------------------------------
      !  3) Update ice salinity and recover ice 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 ice temperature
      DO jk1 = 1, nlay_i
         DO ji = kideb, kiut
            ztmelts       =  -tmut * s_i_b(ji,jk1) + rtt
            ! Conversion q(S,T) -> T (second order equation)
            zaaa          =  cpic
            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk1) / rhoic - lfus
            zccc          =  lfus * ( ztmelts - rtt )
            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
            t_i_b(ji,jk1) =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
            
            ! mask temperature
            zswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 
            t_i_b(ji,jk1) =  zswitch * t_i_b(ji,jk1) + ( 1._wp - zswitch ) * rtt
         END DO 
      END DO 

      !
      CALL wrk_dealloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 )
      CALL wrk_dealloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, 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