New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icethd_ent.F90 in NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE – NEMO

source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icethd_ent.F90 @ 13140

Last change on this file since 13140 was 13140, checked in by orioltp, 4 years ago

dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation: corrections, see #2364
  • Property svn:keywords set to Id
File size: 6.7 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, compute_hfx_err )
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      LOGICAL, INTENT(in)                     ::   compute_hfx_err  ! determines whether to compute diag.
67                                                                    ! error or not
68      !
69      INTEGER  :: ji         !  dummy loop indices
70      INTEGER  :: jk0, jk1   !  old/new layer indices
71      !
72      REAL(wp), DIMENSION(jpij,0:nlay_i+2) ::   zeh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
73      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces
74      REAL(wp), DIMENSION(jpij)            ::   zhnew               ! new layers thicknesses
75      !!-------------------------------------------------------------------
76
77      !--------------------------------------------------------------------------
78      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces
79      !--------------------------------------------------------------------------
80      zeh_cum0(1:npti,0) = 0._wp 
81      zh_cum0 (1:npti,0) = 0._wp
82      DO jk0 = 1, nlay_i+2
83         DO ji = 1, npti
84            zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + eh_i_old(ji,jk0-1)
85            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1)
86         END DO
87      END DO
88
89      !------------------------------------
90      !  2) Interpolation on the new layers
91      !------------------------------------
92      ! new layer thickesses
93      DO ji = 1, npti
94         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i 
95      END DO
96
97      ! new layers interfaces
98      zh_cum1(1:npti,0) = 0._wp
99      DO jk1 = 1, nlay_i
100         DO ji = 1, npti
101            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
102         END DO
103      END DO
104
105      zeh_cum1(1:npti,0:nlay_i) = 0._wp 
106      ! new cumulative q*h => linear interpolation
107      DO jk0 = 1, nlay_i+2
108         DO jk1 = 1, nlay_i-1
109            DO ji = 1, npti
110               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
111                  zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  &
112                     &                 zeh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  &
113                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) )
114               ENDIF
115            END DO
116         END DO
117      END DO
118      ! to ensure that total heat content is strictly conserved, set:
119      zeh_cum1(1:npti,nlay_i) = zeh_cum0(1:npti,nlay_i+2) 
120
121      ! new enthalpies
122      DO jk1 = 1, nlay_i
123         DO ji = 1, npti
124            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 
125            qnew(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
126         END DO
127      END DO
128
129      ! --- diag error on heat remapping --- !
130      ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do),
131      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
132      IF( compute_hfx_err ) THEN
133         DO ji = 1, npti
134            hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice *  &
135               &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) )
136         END DO
137      END IF
138 
139   END SUBROUTINE ice_thd_ent
140
141#else
142   !!----------------------------------------------------------------------
143   !!   Default option                                NO SI3 sea-ice model
144   !!----------------------------------------------------------------------
145#endif
146
147   !!======================================================================
148END MODULE icethd_ent
Note: See TracBrowser for help on using the repository browser.