source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 7.0 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 lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   lim_thd_ent         ! called by limthd and limthd_lac
37
38   !!----------------------------------------------------------------------
39   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44 
45   SUBROUTINE lim_thd_ent( kideb, kiut, qnew )
46      !!-------------------------------------------------------------------
47      !!               ***   ROUTINE lim_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      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
71
72      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (J.m-3, remapped)
73
74      INTEGER  :: ji         !  dummy loop indices
75      INTEGER  :: jk0, jk1   !  old/new layer indices
76      !
77      REAL(wp), DIMENSION(jpij,0:nlay_i+2) :: zqh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
78      REAL(wp), DIMENSION(jpij,0:nlay_i) :: zqh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces
79      REAL(wp), DIMENSION(jpij)   :: zhnew               ! new layers thicknesses
80      !!-------------------------------------------------------------------
81
82
83      !--------------------------------------------------------------------------
84      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces
85      !--------------------------------------------------------------------------
86      zqh_cum0(:,0:nlay_i+2) = 0._wp 
87      zh_cum0 (:,0:nlay_i+2) = 0._wp
88      DO jk0 = 1, nlay_i+2
89         DO ji = kideb, kiut
90            zqh_cum0(ji,jk0) = zqh_cum0(ji,jk0-1) + qh_i_old(ji,jk0-1)
91            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1)
92         ENDDO
93      ENDDO
94
95      !------------------------------------
96      !  2) Interpolation on the new layers
97      !------------------------------------
98      ! new layer thickesses
99      DO ji = kideb, kiut
100         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i 
101      ENDDO
102
103      ! new layers interfaces
104      zh_cum1(:,0:nlay_i) = 0._wp
105      DO jk1 = 1, nlay_i
106         DO ji = kideb, kiut
107            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
108         ENDDO
109      ENDDO
110
111      zqh_cum1(:,0:nlay_i) = 0._wp 
112      ! new cumulative q*h => linear interpolation
113      DO jk0 = 1, nlay_i+1
114         DO jk1 = 1, nlay_i-1
115            DO ji = kideb, kiut
116               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
117                  zqh_cum1(ji,jk1) = ( zqh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  &
118                     &                 zqh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  &
119                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) )
120               ENDIF
121            ENDDO
122         ENDDO
123      ENDDO
124      ! to ensure that total heat content is strictly conserved, set:
125      zqh_cum1(:,nlay_i) = zqh_cum0(:,nlay_i+2) 
126
127      ! new enthalpies
128      DO jk1 = 1, nlay_i
129         DO ji = kideb, kiut
130            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 
131            qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
132         ENDDO
133      ENDDO
134
135      ! --- diag error on heat remapping --- !
136      ! comment: if input h_i_old and qh_i_old are already multiplied by a_i (as in limthd_lac),
137      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
138      DO ji = kideb, kiut
139         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  &
140            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( qh_i_old(ji,0:nlay_i+1) ) ) 
141      END DO
142     
143      !
144      !
145   END SUBROUTINE lim_thd_ent
146
147#else
148   !!----------------------------------------------------------------------
149   !!   Default option                               NO  LIM3 sea-ice model
150   !!----------------------------------------------------------------------
151CONTAINS
152   SUBROUTINE lim_thd_ent          ! Empty routine
153   END SUBROUTINE lim_thd_ent
154#endif
155
156   !!======================================================================
157END MODULE limthd_ent
Note: See TracBrowser for help on using the repository browser.