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 trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 @ 4990

Last change on this file since 4990 was 4990, checked in by timgraham, 9 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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