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_dif.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_dif.F90 @ 8518

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

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