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

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

continue naming

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