source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90 @ 8554

Last change on this file since 8554 was 8534, checked in by clem, 3 years ago

changes in style - part6 - pure cosmetics

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 :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
8   !!                 ! 2005-07 (M. Vancoppenolle) 3D version
9   !!                 ! 2006-11 (X. Fettweis) Vectorized
10   !!            3.0  ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code
11   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
12   !!             -   ! 2014-05 (C. Rousset) complete rewriting
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                       ESIM sea-ice model
17   !!----------------------------------------------------------------------
18   !!   ice_thd_ent   : ice redistribution of enthalpy
19   !!----------------------------------------------------------------------
20   USE dom_oce        ! domain variables
21   USE domain         !
22   USE phycst         ! physical constants
23   USE ice            ! sea-ice: variables
24   USE ice1D          ! sea-ice: thermodynamics variables
25   !
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ice_thd_ent         ! called by icethd and icethd_do
34
35   !!----------------------------------------------------------------------
36   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
37   !! $Id: icethd_ent.F90 8420 2017-08-08 12:18:46Z clem $
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41 
42   SUBROUTINE ice_thd_ent( qnew )
43      !!-------------------------------------------------------------------
44      !!               ***   ROUTINE ice_thd_ent  ***
45      !!
46      !! ** Purpose :
47      !!           This routine computes new vertical grids in the ice,
48      !!           and consistently redistributes temperatures.
49      !!           Redistribution is made so as to ensure to energy conservation
50      !!
51      !!
52      !! ** Method  : linear conservative remapping
53      !!           
54      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
55      !!            2) linear remapping on the new layers
56      !!
57      !! ------------ cum0(0)                        ------------- cum1(0)
58      !!                                    NEW      -------------
59      !! ------------ cum0(1)               ==>      -------------
60      !!     ...                                     -------------
61      !! ------------                                -------------
62      !! ------------ cum0(nlay_i+2)                 ------------- cum1(nlay_i)
63      !!
64      !!
65      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
66      !!-------------------------------------------------------------------
67      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   qnew             ! new enthlapies (J.m-3, remapped)
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(:,0:nlay_i+2) = 0._wp 
81      zh_cum0 (:,0:nlay_i+2) = 0._wp
82      DO jk0 = 1, nlay_i+2
83         DO ji = 1, nidx
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, nidx
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(:,0:nlay_i) = 0._wp
99      DO jk1 = 1, nlay_i
100         DO ji = 1, nidx
101            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
102         END DO
103      END DO
104
105      zeh_cum1(:,0:nlay_i) = 0._wp 
106      ! new cumulative q*h => linear interpolation
107      DO jk0 = 1, nlay_i+1
108         DO jk1 = 1, nlay_i-1
109            DO ji = 1, nidx
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(:,nlay_i) = zeh_cum0(:,nlay_i+2) 
120
121      ! new enthalpies
122      DO jk1 = 1, nlay_i
123         DO ji = 1, nidx
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      DO ji = 1, nidx
133         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  &
134            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 
135      END DO
136     
137   END SUBROUTINE ice_thd_ent
138
139#else
140   !!----------------------------------------------------------------------
141   !!   Default option                                NO ESIM sea-ice model
142   !!----------------------------------------------------------------------
143#endif
144
145   !!======================================================================
146END MODULE icethd_ent
Note: See TracBrowser for help on using the repository browser.