source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_ent.F90 @ 12726

Last change on this file since 12726 was 12726, checked in by clem, 8 months ago

clean useless variables

  • Property svn:keywords set to Id
File size: 6.5 KB
Line 
1MODULE icethd_ent
2   !!======================================================================
3   !!                       ***  MODULE icethd_ent   ***
4   !!   sea-ice: redistribution of enthalpy in the ice on the new vertical grid
5   !!                       after vertical growth/melt
6   !!======================================================================
7   !! History :       !  2003-05  (M. Vancoppenolle) Original code in 1D
8   !!                 !  2005-07  (M. Vancoppenolle) 3D version
9   !!            3.6  !  2014-05  (C. Rousset)       New version
10   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
11   !!----------------------------------------------------------------------
12#if defined key_si3
13   !!----------------------------------------------------------------------
14   !!   'key_si3'                                       SI3 sea-ice model
15   !!----------------------------------------------------------------------
16   !!   ice_thd_ent   : ice redistribution of enthalpy
17   !!----------------------------------------------------------------------
18   USE dom_oce        ! domain variables
19   USE domain         !
20   USE phycst         ! physical constants
21   USE ice            ! sea-ice: variables
22   USE ice1D          ! sea-ice: thermodynamics variables
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   ice_thd_ent         ! called by icethd and icethd_do
32
33   !!----------------------------------------------------------------------
34   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
35   !! $Id$
36   !! Software governed by the CeCILL license (see ./LICENSE)
37   !!----------------------------------------------------------------------
38CONTAINS
39 
40   SUBROUTINE ice_thd_ent( qnew )
41      !!-------------------------------------------------------------------
42      !!               ***   ROUTINE ice_thd_ent  ***
43      !!
44      !! ** Purpose :
45      !!           This routine computes new vertical grids in the ice,
46      !!           and consistently redistributes temperatures.
47      !!           Redistribution is made so as to ensure to energy conservation
48      !!
49      !!
50      !! ** Method  : linear conservative remapping
51      !!           
52      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
53      !!            2) linear remapping on the new layers
54      !!
55      !! ------------ cum0(0)                        ------------- cum1(0)
56      !!                                    NEW      -------------
57      !! ------------ cum0(1)               ==>      -------------
58      !!     ...                                     -------------
59      !! ------------                                -------------
60      !! ------------ cum0(nlay_i+2)                 ------------- cum1(nlay_i)
61      !!
62      !!
63      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
64      !!-------------------------------------------------------------------
65      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   qnew             ! new enthlapies (J.m-3, remapped)
66      !
67      INTEGER  :: ji         !  dummy loop indices
68      INTEGER  :: jk0, jk1   !  old/new layer indices
69      !
70      REAL(wp), DIMENSION(jpij,0:nlay_i+2) ::   zeh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
71      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces
72      REAL(wp), DIMENSION(jpij)            ::   zhnew               ! new layers thicknesses
73      !!-------------------------------------------------------------------
74
75      !--------------------------------------------------------------------------
76      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces
77      !--------------------------------------------------------------------------
78      zeh_cum0(1:npti,0) = 0._wp 
79      zh_cum0 (1:npti,0) = 0._wp
80      DO jk0 = 1, nlay_i+2
81         DO ji = 1, npti
82            zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + eh_i_old(ji,jk0-1)
83            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1)
84         END DO
85      END DO
86
87      !------------------------------------
88      !  2) Interpolation on the new layers
89      !------------------------------------
90      ! new layer thickesses
91      DO ji = 1, npti
92         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i 
93      END DO
94
95      ! new layers interfaces
96      zh_cum1(1:npti,0) = 0._wp
97      DO jk1 = 1, nlay_i
98         DO ji = 1, npti
99            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
100         END DO
101      END DO
102
103      zeh_cum1(1:npti,0:nlay_i) = 0._wp 
104      ! new cumulative q*h => linear interpolation
105      DO jk0 = 1, nlay_i+2
106         DO jk1 = 1, nlay_i-1
107            DO ji = 1, npti
108               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
109                  zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  &
110                     &                 zeh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  &
111                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) )
112               ENDIF
113            END DO
114         END DO
115      END DO
116      ! to ensure that total heat content is strictly conserved, set:
117      zeh_cum1(1:npti,nlay_i) = zeh_cum0(1:npti,nlay_i+2) 
118
119      ! new enthalpies
120      DO jk1 = 1, nlay_i
121         DO ji = 1, npti
122            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 
123            qnew(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
124         END DO
125      END DO
126
127      ! --- diag error on heat remapping --- !
128      ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do),
129      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
130      !DO ji = 1, npti
131      !   hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  &
132      !      &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) )
133      !END DO
134     
135   END SUBROUTINE ice_thd_ent
136
137#else
138   !!----------------------------------------------------------------------
139   !!   Default option                                NO SI3 sea-ice model
140   !!----------------------------------------------------------------------
141#endif
142
143   !!======================================================================
144END MODULE icethd_ent
Note: See TracBrowser for help on using the repository browser.