source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

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