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

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

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

File size: 6.7 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   !
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   ice_thd_ent         ! called by icethd and icethd_lac
36
37   !!----------------------------------------------------------------------
38   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
39   !! $Id: icethd_ent.F90 8420 2017-08-08 12:18:46Z clem $
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43 
44   SUBROUTINE ice_thd_ent( qnew )
45      !!-------------------------------------------------------------------
46      !!               ***   ROUTINE ice_thd_ent  ***
47      !!
48      !! ** Purpose :
49      !!           This routine computes new vertical grids in the ice,
50      !!           and consistently redistributes temperatures.
51      !!           Redistribution is made so as to ensure to energy conservation
52      !!
53      !!
54      !! ** Method  : linear conservative remapping
55      !!           
56      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
57      !!            2) linear remapping on the new layers
58      !!
59      !! ------------ cum0(0)                        ------------- cum1(0)
60      !!                                    NEW      -------------
61      !! ------------ cum0(1)               ==>      -------------
62      !!     ...                                     -------------
63      !! ------------                                -------------
64      !! ------------ cum0(nlay_i+2)                 ------------- cum1(nlay_i)
65      !!
66      !!
67      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
68      !!-------------------------------------------------------------------
69      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (J.m-3, remapped)
70
71      INTEGER  :: ji         !  dummy loop indices
72      INTEGER  :: jk0, jk1   !  old/new layer indices
73      !
74      REAL(wp), DIMENSION(jpij,0:nlay_i+2) :: zeh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
75      REAL(wp), DIMENSION(jpij,0:nlay_i)   :: zeh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces
76      REAL(wp), DIMENSION(jpij)            :: zhnew               ! new layers thicknesses
77      !!-------------------------------------------------------------------
78
79      !--------------------------------------------------------------------------
80      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces
81      !--------------------------------------------------------------------------
82      zeh_cum0(:,0:nlay_i+2) = 0._wp 
83      zh_cum0 (:,0:nlay_i+2) = 0._wp
84      DO jk0 = 1, nlay_i+2
85         DO ji = 1, nidx
86            zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + eh_i_old(ji,jk0-1)
87            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1)
88         END DO
89      END DO
90
91      !------------------------------------
92      !  2) Interpolation on the new layers
93      !------------------------------------
94      ! new layer thickesses
95      DO ji = 1, nidx
96         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i 
97      END DO
98
99      ! new layers interfaces
100      zh_cum1(:,0:nlay_i) = 0._wp
101      DO jk1 = 1, nlay_i
102         DO ji = 1, nidx
103            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
104         END DO
105      END DO
106
107      zeh_cum1(:,0:nlay_i) = 0._wp 
108      ! new cumulative q*h => linear interpolation
109      DO jk0 = 1, nlay_i+1
110         DO jk1 = 1, nlay_i-1
111            DO ji = 1, nidx
112               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
113                  zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  &
114                     &                 zeh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  &
115                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) )
116               ENDIF
117            END DO
118         END DO
119      END DO
120      ! to ensure that total heat content is strictly conserved, set:
121      zeh_cum1(:,nlay_i) = zeh_cum0(:,nlay_i+2) 
122
123      ! new enthalpies
124      DO jk1 = 1, nlay_i
125         DO ji = 1, nidx
126            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 
127            qnew(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
128         END DO
129      END DO
130
131      ! --- diag error on heat remapping --- !
132      ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_lac),
133      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
134      DO ji = 1, nidx
135         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  &
136            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 
137      END DO
138     
139   END SUBROUTINE ice_thd_ent
140
141#else
142   !!----------------------------------------------------------------------
143   !!   Default option                               NO  LIM3 sea-ice model
144   !!----------------------------------------------------------------------
145#endif
146
147   !!======================================================================
148END MODULE icethd_ent
Note: See TracBrowser for help on using the repository browser.