Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
MODULE icethd_ent
!!======================================================================
!! *** MODULE icethd_ent ***
!! sea-ice: redistribution of enthalpy in the ice on the new vertical grid
!! after vertical growth/melt
!!======================================================================
!! History : ! 2003-05 (M. Vancoppenolle) Original code in 1D
!! ! 2005-07 (M. Vancoppenolle) 3D version
!! 3.6 ! 2014-05 (C. Rousset) New version
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_thd_ent : ice redistribution of enthalpy
!!----------------------------------------------------------------------
USE dom_oce ! domain variables
USE domain !
USE phycst ! physical constants
USE ice ! sea-ice: variables
USE ice1D ! sea-ice: thermodynamics variables
!
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_thd_ent ! called by icethd and icethd_do
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icethd_ent.F90 14778 2021-05-03 08:58:22Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_thd_ent( ph_old, pe_old, pe_new )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_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
!!
!! ------------ cum0(0) ------------- cum1(0)
!! NEW -------------
!! ------------ cum0(1) ==> -------------
!! ... -------------
!! ------------ -------------
!! ------------ cum0(nlay_i+2) ------------- cum1(nlay_i)
!!
!!
!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(jpij,0:nlay_i+1), INTENT(in) :: ph_old, pe_old ! old tickness and enthlapy
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pe_new ! new enthlapies (J.m-3, remapped)
!
INTEGER :: ji ! dummy loop indices
INTEGER :: jk0, jk1 ! old/new layer indices
!
REAL(wp), DIMENSION(0:nlay_i+2) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces
REAL(wp), DIMENSION(0:nlay_i) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces
REAL(wp) :: zhnew ! new layers thicknesses
!!-------------------------------------------------------------------
DO ji = 1, npti
!--------------------------------------------------------------------------
! 1) Cumulative integral of old enthalpy * thickness and layers interfaces
!--------------------------------------------------------------------------
zeh_cum0(0) = 0._wp
zh_cum0 (0) = 0._wp
DO jk0 = 1, nlay_i+2
zeh_cum0(jk0) = zeh_cum0(jk0-1) + pe_old(ji,jk0-1)
zh_cum0 (jk0) = zh_cum0 (jk0-1) + ph_old(ji,jk0-1)
!------------------------------------
! 2) Interpolation on the new layers
!------------------------------------
! new layer thickesses
zhnew = SUM( ph_old(ji,0:nlay_i+1) ) * r1_nlay_i
! new layers interfaces
zh_cum1(0) = 0._wp
DO jk1 = 1, nlay_i
zh_cum1(jk1) = zh_cum1(jk1-1) + zhnew
zeh_cum1(0:nlay_i) = 0._wp
! new cumulative q*h => linear interpolation
DO jk0 = 1, nlay_i+2
DO jk1 = 1, nlay_i-1
IF( zh_cum1(jk1) <= zh_cum0(jk0) .AND. zh_cum1(jk1) > zh_cum0(jk0-1) ) THEN
zeh_cum1(jk1) = ( zeh_cum0(jk0-1) * ( zh_cum0(jk0) - zh_cum1(jk1 ) ) + &
& zeh_cum0(jk0 ) * ( zh_cum1(jk1) - zh_cum0(jk0-1) ) ) &
& / ( zh_cum0(jk0) - zh_cum0(jk0-1) )
! to ensure that total heat content is strictly conserved, set:
zeh_cum1(nlay_i) = zeh_cum0(nlay_i+2)
! new enthalpies
DO jk1 = 1, nlay_i
rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew - epsi20 ) )
pe_new(ji,jk1) = rswitch * MAX( 0._wp, zeh_cum1(jk1) - zeh_cum1(jk1-1) ) / MAX( zhnew, epsi20 ) ! max for roundoff error
! --- diag error on heat remapping --- !
! comment: if input h_old and eh_old are already multiplied by a_i (as in icethd_do),
! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
! hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice * &
! & ( SUM( pe_new(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_old(ji,0:nlay_i+1) ) )
END DO
END SUBROUTINE ice_thd_ent
#else
!!----------------------------------------------------------------------
!! Default option NO SI3 sea-ice model
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE icethd_ent