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.
icethd_zdf.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90 @ 8562

Last change on this file since 8562 was 8562, checked in by clem, 7 years ago

cosmetics only

File size: 35.8 KB
Line 
1MODULE icethd_zdf
2   !!======================================================================
3   !!                       ***  MODULE icethd_zdf ***
4   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
5   !!======================================================================
6   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
7   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
8   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
9   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
10   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
11   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
12   !!----------------------------------------------------------------------
13#if defined key_lim3
14   !!----------------------------------------------------------------------
15   !!   'key_lim3'                                       ESIM sea-ice model
16   !!----------------------------------------------------------------------
17   USE dom_oce        ! ocean space and time domain
18   USE phycst         ! physical constants (ocean directory)
19   USE ice            ! sea-ice: variables
20   USE ice1D          ! sea-ice: thermodynamics variables
21   !
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   ice_thd_zdf        ! called by icethd
30   PUBLIC   ice_thd_zdf_init   ! called by icestp
31
32   !!** namelist (namthd_zdf) **
33   LOGICAL  ::   ln_zdf_Beer      ! Heat diffusion follows a Beer Law
34   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964)
35   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007)
36   REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K]
37   REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
38   LOGICAL  ::   ln_dqns_i        ! change non-solar surface flux with changing surface temperature (T) or not (F)
39
40   !!----------------------------------------------------------------------
41   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
42   !! $Id: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE ice_thd_zdf
48      !!-------------------------------------------------------------------
49      !!                ***  ROUTINE ice_thd_zdf  ***
50      !! ** Purpose :
51      !!           This routine determines the time evolution of snow and sea-ice
52      !!           temperature profiles.
53      !! ** Method  :
54      !!           This is done by solving the heat equation diffusion with
55      !!           a Neumann boundary condition at the surface and a Dirichlet one
56      !!           at the bottom. Solar radiation is partially absorbed into the ice.
57      !!           The specific heat and thermal conductivities depend on ice salinity
58      !!           and temperature to take into account brine pocket melting. The
59      !!           numerical
60      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
61      !!           in the ice and snow system.
62      !!
63      !!           The successive steps of this routine are
64      !!           1.  initialization of ice-snow layers thicknesses
65      !!           2.  Internal absorbed and transmitted radiation
66      !!           Then iterative procedure begins
67      !!           3.  Thermal conductivity
68      !!           4.  Kappa factors
69      !!           5.  specific heat in the ice
70      !!           6.  eta factors
71      !!           7.  surface flux computation
72      !!           8.  tridiagonal system terms
73      !!           9.  solving the tridiagonal system with Gauss elimination
74      !!           Iterative procedure ends according to a criterion on evolution
75      !!           of temperature
76      !!           10. Fluxes at the interfaces
77      !!
78      !! ** Inputs / Ouputs : (global commons)
79      !!           surface temperature : t_su_1d
80      !!           ice/snow temperatures   : t_i_1d, t_s_1d
81      !!           ice salinities          : s_i_1d
82      !!           number of layers in the ice/snow: nlay_i, nlay_s
83      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
84      !!-------------------------------------------------------------------
85      INTEGER ::   ji, jk         ! spatial loop index
86      INTEGER ::   jm             ! current reference number of equation
87      INTEGER ::   jm_mint, jm_maxt
88      INTEGER ::   iconv          ! number of iterations in iterative procedure
89      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
90     
91      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
92      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
93     
94      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
95      REAL(wp) ::   zg1       =  2._wp        !
96      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
97      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
98      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
99      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
100      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C
101      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature
102      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
103      REAL(wp) ::   z1_hsu
104      REAL(wp) ::   zdti_max                  ! current maximal error on temperature
105      REAL(wp) ::   zcpi                      ! Ice specific heat
106      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
107      REAL(wp) ::   zfac                      ! dummy factor
108     
109      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
110      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
111      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
112      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
113      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
114      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
115      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function
116      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
117      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
118     
119      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
120      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
121      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
122      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
123      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
124      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
125      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
126      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
127      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
128      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
129      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
130      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
131      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
132      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
133      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
134      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
135      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
136      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
137      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
138     
139      ! Mono-category
140      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
141      REAL(wp) :: zhe           ! dummy factor
142      REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity
143      !!------------------------------------------------------------------     
144
145      ! --- diag error on heat diffusion - PART 1 --- !
146      DO ji = 1, nidx
147         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
148            &           SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
149      END DO
150
151      !------------------
152      ! 1) Initialization
153      !------------------
154      DO ji = 1, nidx
155         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
156         ! layer thickness
157         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
158         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
159      END DO
160      !
161      WHERE( zh_i(1:nidx) >= epsi10 )   ;   z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx)
162      ELSEWHERE                         ;   z1_h_i(1:nidx) = 0._wp
163      END WHERE
164
165      WHERE( zh_s(1:nidx) >= epsi10 )   ;   z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx)
166      ELSEWHERE                         ;   z1_h_s(1:nidx) = 0._wp
167      END WHERE
168      !
169      ! temperatures
170      ztsub  (1:nidx)   = t_su_1d(1:nidx)  ! temperature at the previous iteration
171      ztsold (1:nidx,:) = t_s_1d(1:nidx,:)   ! Old snow temperature
172      ztiold (1:nidx,:) = t_i_1d(1:nidx,:)   ! Old ice temperature
173      t_su_1d(1:nidx)   = MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary
174      !
175      !-------------
176      ! 2) Radiation
177      !-------------
178      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
179      DO ji = 1, nidx
180         ! --- Computation of i0 --- !
181         ! i0 describes the fraction of solar radiation which does not contribute
182         ! to the surface energy budget but rather penetrates inside the ice.
183         ! We assume that no radiation is transmitted through the snow
184         ! If there is no no snow
185         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
186         ! zftrice = io.qsr_ice       is below the surface
187         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
188         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
189         zfac = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )     
190         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )
191
192         ! --- Solar radiation absorbed / transmitted at the surface --- !
193         !     Derivative of the non solar flux
194         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
195         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
196         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
197         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
198      END DO
199
200      ! --- Transmission/absorption of solar radiation in the ice --- !
201      zradtr_s(1:nidx,0) = zftrice(1:nidx)
202      DO jk = 1, nlay_s
203         DO ji = 1, nidx
204            !                             ! radiation transmitted below the layer-th snow layer
205            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
206            !                             ! radiation absorbed by the layer-th snow layer
207            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
208         END DO
209      END DO
210         
211      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )
212      DO jk = 1, nlay_i 
213         DO ji = 1, nidx
214            !                             ! radiation transmitted below the layer-th ice layer
215            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
216            !                             ! radiation absorbed by the layer-th ice layer
217            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
218         END DO
219      END DO
220
221      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice
222      !
223      iconv    =  0          ! number of iterations
224      zdti_max =  1000._wp   ! maximal value of error on all points
225      !                                                          !============================!
226      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins !
227         !                                                       !============================!
228         iconv = iconv + 1
229         !
230         ztib(1:nidx,:) = t_i_1d(1:nidx,:)
231         ztsb(1:nidx,:) = t_s_1d(1:nidx,:)
232         !
233         !--------------------------------
234         ! 3) Sea ice thermal conductivity
235         !--------------------------------
236         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
237            !
238            DO ji = 1, nidx
239               ztcond_i(ji,0)      = rcdic + zbeta * s_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
240               ztcond_i(ji,nlay_i) = rcdic + zbeta * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
241            END DO
242            DO jk = 1, nlay_i-1
243               DO ji = 1, nidx
244                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
245                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
246               END DO
247            END DO
248            !
249         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
250            !
251            DO ji = 1, nidx
252               ztcond_i(ji,0)      = rcdic + 0.09_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
253                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
254               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )  &
255                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
256            END DO
257            DO jk = 1, nlay_i-1
258               DO ji = 1, nidx
259                  ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /          &
260                     &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   &
261                     &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
262               END DO
263            END DO
264            !
265         ENDIF
266         ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )       
267         !
268         !--- G(he) : enhancement of thermal conductivity in mono-category case
269         ! Computation of effective thermal conductivity G(h)
270         ! Used in mono-category case only to simulate an ITD implicitly
271         ! Fichefet and Morales Maqueda, JGR 1997
272         zghe(1:nidx) = 1._wp
273         !
274         SELECT CASE ( nn_monocat )
275
276         CASE ( 1 , 3 )
277
278            zepsilon = 0.1_wp
279            DO ji = 1, nidx
280               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                             ! Mean sea ice thermal conductivity
281               zhe = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )   ! Effective thickness he (zhe)
282               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
283                  zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )    ! G(he)
284               ENDIF
285            END DO
286
287         END SELECT
288         !
289         !-----------------
290         ! 4) kappa factors
291         !-----------------
292         !--- Snow
293         DO jk = 0, nlay_s-1
294            DO ji = 1, nidx
295               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
296            END DO
297         END DO
298         DO ji = 1, nidx   ! Snow-ice interface
299            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
300            IF( zfac > epsi10 ) THEN
301               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
302            ELSE
303               zkappa_s(ji,nlay_s) = 0._wp
304            ENDIF
305         END DO
306
307         !--- Ice
308         DO jk = 0, nlay_i
309            DO ji = 1, nidx
310               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
311            END DO
312         END DO
313         DO ji = 1, nidx   ! Snow-ice interface
314            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
315         END DO
316         !
317         !--------------------------------------
318         ! 5) Sea ice specific heat, eta factors
319         !--------------------------------------
320         DO jk = 1, nlay_i
321            DO ji = 1, nidx
322               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
323               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
324            END DO
325         END DO
326
327         DO jk = 1, nlay_s
328            DO ji = 1, nidx
329               zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
330            END DO
331         END DO
332         !
333         !----------------------------
334         ! 6) surface flux computation
335         !----------------------------
336         IF ( ln_dqns_i ) THEN
337            DO ji = 1, nidx
338               ! update of the non solar flux according to the update in T_su
339               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
340            END DO
341         ENDIF
342
343         DO ji = 1, nidx
344            zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
345         END DO
346         !
347         !----------------------------
348         ! 7) tridiagonal system terms
349         !----------------------------
350         !!layer denotes the number of the layer in the snow or in the ice
351         !!jm denotes the reference number of the equation in the tridiagonal
352         !!system, terms of tridiagonal system are indexed as following :
353         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
354
355         !!ice interior terms (top equation has the same form as the others)
356         ztrid   (1:nidx,:,:) = 0._wp
357         zindterm(1:nidx,:)   = 0._wp
358         zindtbis(1:nidx,:)   = 0._wp
359         zdiagbis(1:nidx,:)   = 0._wp
360
361         DO jm = nlay_s + 2, nlay_s + nlay_i 
362            DO ji = 1, nidx
363               jk = jm - nlay_s - 1
364               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
365               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
366               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
367               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
368            END DO
369         ENDDO
370
371         jm =  nlay_s + nlay_i + 1
372         DO ji = 1, nidx
373            !!ice bottom term
374            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
375            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
376            ztrid(ji,jm,3)  = 0.0
377            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
378               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
379         ENDDO
380
381
382         DO ji = 1, nidx
383            !                               !---------------------!
384            IF ( ht_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
385               !                            !---------------------!
386               ! snow interior terms (bottom equation has the same form as the others)
387               DO jm = 3, nlay_s + 1
388                  jk = jm - 1
389                  ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
390                  ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
391                  ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
392                  zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
393               END DO
394
395               ! case of only one layer in the ice (ice equation is altered)
396               IF ( nlay_i == 1 ) THEN
397                  ztrid(ji,nlay_s+2,3)  = 0.0
398                  zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
399               ENDIF
400
401               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
402
403                  jm_min(ji) = 1
404                  jm_max(ji) = nlay_i + nlay_s + 1
405
406                  ! surface equation
407                  ztrid(ji,1,1)  = 0.0
408                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
409                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
410                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
411
412                  ! first layer of snow equation
413                  ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
414                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
415                  ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1)
416                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
417
418               ELSE                            !--  case 2 : surface is melting
419                  !
420                  jm_min(ji) = 2
421                  jm_max(ji) = nlay_i + nlay_s + 1
422
423                  ! first layer of snow equation
424                  ztrid(ji,2,1)  = 0.0
425                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
426                  ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
427                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
428                     &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
429               ENDIF
430            !                               !---------------------!
431            ELSE                            ! cells without snow  !
432               !                            !---------------------!
433               !
434               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
435                  !
436                  jm_min(ji) = nlay_s + 1
437                  jm_max(ji) = nlay_i + nlay_s + 1
438
439                  ! surface equation   
440                  ztrid(ji,jm_min(ji),1)  = 0.0
441                  ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1   
442                  ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1
443                  zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji)
444
445                  ! first layer of ice equation
446                  ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
447                  ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
448                  ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
449                  zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
450
451                  ! case of only one layer in the ice (surface & ice equations are altered)
452                  IF ( nlay_i == 1 ) THEN
453                     ztrid(ji,jm_min(ji),1)    = 0.0
454                     ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
455                     ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0
456                     ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
457                     ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
458                     ztrid(ji,jm_min(ji)+1,3)  = 0.0
459                     zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  &
460                        &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
461                  ENDIF
462
463               ELSE                            !--  case 2 : surface is melting
464
465                  jm_min(ji)    =  nlay_s + 2
466                  jm_max(ji)    =  nlay_i + nlay_s + 1
467
468                  ! first layer of ice equation
469                  ztrid(ji,jm_min(ji),1)  = 0.0
470                  ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
471                  ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
472                  zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
473                     &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
474
475                  ! case of only one layer in the ice (surface & ice equations are altered)
476                  IF ( nlay_i == 1 ) THEN
477                     ztrid(ji,jm_min(ji),1)  = 0.0
478                     ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
479                     ztrid(ji,jm_min(ji),3)  = 0.0
480                     zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
481                        &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
482                  ENDIF
483
484               ENDIF
485            ENDIF
486            !
487            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
488            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
489            !
490         END DO
491         !
492         !------------------------------
493         ! 8) tridiagonal system solving
494         !------------------------------
495         ! Solve the tridiagonal system with Gauss elimination method.
496         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
497         jm_maxt = 0
498         jm_mint = nlay_i+5
499         DO ji = 1, nidx
500            jm_mint = MIN(jm_min(ji),jm_mint)
501            jm_maxt = MAX(jm_max(ji),jm_maxt)
502         END DO
503
504         DO jk = jm_mint+1, jm_maxt
505            DO ji = 1, nidx
506               jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
507               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
508               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
509            END DO
510         END DO
511
512         DO ji = 1, nidx
513            ! ice temperatures
514            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
515         END DO
516
517         DO jm = nlay_i + nlay_s, nlay_s + 2, -1
518            DO ji = 1, nidx
519               jk    =  jm - nlay_s - 1
520               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
521            END DO
522         END DO
523
524         DO ji = 1, nidx
525            ! snow temperatures     
526            IF( ht_s_1d(ji) > 0._wp ) THEN
527               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
528                  &                / zdiagbis(ji,nlay_s+1)
529            ENDIF
530            ! surface temperature
531            ztsub(ji) = t_su_1d(ji)
532            IF( t_su_1d(ji) < rt0 ) THEN
533               t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
534                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
535            ENDIF
536         END DO
537         !
538         !--------------------------------------------------------------
539         ! 9) Has the scheme converged ?, end of the iterative procedure
540         !--------------------------------------------------------------
541         ! check that nowhere it has started to melt
542         ! zdti_max is a measure of error, it has to be under zdti_bnd
543         zdti_max = 0._wp
544         DO ji = 1, nidx
545            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
546            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
547         END DO
548
549         DO jk  =  1, nlay_s
550            DO ji = 1, nidx
551               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
552               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
553            END DO
554         END DO
555
556         DO jk  =  1, nlay_i
557            DO ji = 1, nidx
558               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
559               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
560               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
561            END DO
562         END DO
563
564         ! Compute spatial maximum over all errors
565         ! note that this could be optimized substantially by iterating only the non-converging points
566         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
567
568      END DO  ! End of the do while iterative procedure
569
570      IF( ln_icectl .AND. lwp ) THEN
571         WRITE(numout,*) ' zdti_max : ', zdti_max
572         WRITE(numout,*) ' iconv    : ', iconv
573      ENDIF
574      !
575      !-----------------------------
576      ! 10) Fluxes at the interfaces
577      !-----------------------------
578      DO ji = 1, nidx
579         !                                ! surface ice conduction flux
580         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
581            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
582         !                                ! bottom ice conduction flux
583         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
584      END DO
585
586      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
587      CALL ice_thd_enmelt
588
589      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
590      IF ( ln_dqns_i ) THEN
591         DO ji = 1, nidx
592            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
593         END DO
594      END IF
595
596      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
597      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
598      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
599      DO ji = 1, nidx
600         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
601            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
602
603         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
604            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
605         ELSE                          ! case T_su = 0degC
606            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
607         ENDIF
608         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
609
610         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
611         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
612
613      END DO   
614
615      ! --- Diagnostics SIMIP --- !
616      DO ji = 1, nidx
617         !--- Conduction fluxes (positive downwards)
618         diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
619         diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
620
621         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
622         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
623         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
624            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
625               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
626         ELSE
627            t_si_1d(ji) = t_su_1d(ji)
628         ENDIF
629      END DO
630      !
631   END SUBROUTINE ice_thd_zdf
632
633
634   SUBROUTINE ice_thd_enmelt
635      !!-------------------------------------------------------------------
636      !!                   ***  ROUTINE ice_thd_enmelt ***
637      !!                 
638      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
639      !!
640      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
641      !!-------------------------------------------------------------------
642      INTEGER  ::   ji, jk   ! dummy loop indices
643      REAL(wp) ::   ztmelts  ! local scalar
644      !!-------------------------------------------------------------------
645      !
646      DO jk = 1, nlay_i             ! Sea ice energy of melting
647         DO ji = 1, nidx
648            ztmelts      = - tmut  * s_i_1d(ji,jk)
649            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
650                                                                !   (sometimes dif scheme produces abnormally high temperatures)   
651            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
652               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
653               &                    - rcp  *   ztmelts )
654         END DO
655      END DO
656      DO jk = 1, nlay_s             ! Snow energy of melting
657         DO ji = 1, nidx
658            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
659         END DO
660      END DO
661      !
662   END SUBROUTINE ice_thd_enmelt
663
664
665   SUBROUTINE ice_thd_zdf_init
666      !!-----------------------------------------------------------------------
667      !!                   ***  ROUTINE ice_thd_zdf_init ***
668      !!                 
669      !! ** Purpose :   Physical constants and parameters associated with
670      !!                ice thermodynamics
671      !!
672      !! ** Method  :   Read the namthd_zdf namelist and check the parameters
673      !!                called at the first timestep (nit000)
674      !!
675      !! ** input   :   Namelist namthd_zdf
676      !!-------------------------------------------------------------------
677      INTEGER  ::   ios   ! Local integer output status for namelist read
678      !!
679      NAMELIST/namthd_zdf/ ln_zdf_Beer, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i
680      !!-------------------------------------------------------------------
681      !
682      REWIND( numnam_ice_ref )              ! Namelist namthd_zdf in reference namelist : Ice thermodynamics
683      READ  ( numnam_ice_ref, namthd_zdf, IOSTAT = ios, ERR = 901)
684901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in reference namelist', lwp )
685
686      REWIND( numnam_ice_cfg )              ! Namelist namthd_zdf in configuration namelist : Ice thermodynamics
687      READ  ( numnam_ice_cfg, namthd_zdf, IOSTAT = ios, ERR = 902 )
688902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in configuration namelist', lwp )
689      IF(lwm) WRITE ( numoni, namthd_zdf )
690      !
691      !
692      IF(lwp) THEN                          ! control print
693         WRITE(numout,*) 'ice_thd_zdf_init: Ice vertical heat diffusion'
694         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
695         WRITE(numout,*) '   Namelist namthd_zdf:'
696         WRITE(numout,*) '      Diffusion follows a Beer Law                            ln_zdf_Beer  = ', ln_zdf_Beer
697         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64
698         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07
699         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
700         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
701         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i
702     ENDIF
703      !
704      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
705         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' )
706      ENDIF
707      !
708   END SUBROUTINE ice_thd_zdf_init
709
710#else
711   !!----------------------------------------------------------------------
712   !!   Default option       Dummy Module             No ESIM sea-ice model
713   !!----------------------------------------------------------------------
714#endif
715
716   !!======================================================================
717END MODULE icethd_zdf
Note: See TracBrowser for help on using the repository browser.