New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limthd_ent.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 4656

Last change on this file since 4656 was 4649, checked in by clem, 10 years ago

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