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 @ 8514

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

changes in style - part5 - almost done
File size: 40.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) ::   zhsu
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     
100      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
101      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! old surface temperature (before the iterative procedure )
102      REAL(wp), DIMENSION(jpij)     ::   ztsubit     ! surface temperature at previous iteration
103      REAL(wp), DIMENSION(jpij)     ::   zh_i        ! ice layer thickness
104      REAL(wp), DIMENSION(jpij)     ::   zh_s        ! snow layer thickness
105      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
106      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
107      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function
108      REAL(wp), DIMENSION(jpij)     ::   dzf         ! derivative of the surface flux function
109      REAL(wp), DIMENSION(jpij)     ::   zdti        ! current error on temperature
110      REAL(wp), DIMENSION(jpij)     ::   zdifcase    ! case of the equation resolution (1->4)
111      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
112      REAL(wp), DIMENSION(jpij)     ::   zihic
113     
114      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
115      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
116      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
117      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
118      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztib        ! Old temperature in the ice
119      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
120      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
121      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zspeche_i   ! Ice specific heat
122      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice
123      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
124      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
125      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
126      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
127      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
128      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztsb        ! Temporary temperature in the snow
129      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow
130      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
131      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
132      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
133      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
134      REAL(wp), DIMENSION(jpij)            ::   zdq, zq_ini, zhfx_err ! diag errors on heat
135      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
136     
137      ! Mono-category
138      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
139      REAL(wp) :: zratio_s      ! dummy factor
140      REAL(wp) :: zratio_i      ! dummy factor
141      REAL(wp) :: zh_thres      ! thickness thres. for G(h) computation
142      REAL(wp) :: zhe           ! dummy factor
143      REAL(wp) :: zkimean       ! mean sea ice thermal conductivity
144      REAL(wp) :: zfac          ! dummy factor
145      REAL(wp) :: zihe          ! dummy factor
146      REAL(wp) :: zheshth       ! dummy factor
147      !!------------------------------------------------------------------     
148
149      ! --- diag error on heat diffusion - PART 1 --- !
150      zdq(:) = 0._wp ; zq_ini(:) = 0._wp     
151      DO ji = 1, nidx
152         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
153            &           SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
154      END DO
155
156      !------------------------------------------------------------------------------!
157      ! 1) Initialization                                                            !
158      !------------------------------------------------------------------------------!
159      DO ji = 1 , nidx
160         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
161         ! layer thickness
162         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
163         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
164      END DO
165
166      !--------------------
167      ! Ice / snow layers
168      !--------------------
169
170      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
171      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
172
173      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
174         DO ji = 1 , nidx
175            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
176         END DO
177      END DO
178
179      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
180         DO ji = 1 , nidx
181            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
182         END DO
183      END DO
184      !
185      !------------------------------------------------------------------------------|
186      ! 2) Radiation                                                       |
187      !------------------------------------------------------------------------------|
188      !
189      !-------------------
190      ! Computation of i0
191      !-------------------
192      ! i0 describes the fraction of solar radiation which does not contribute
193      ! to the surface energy budget but rather penetrates inside the ice.
194      ! We assume that no radiation is transmitted through the snow
195      ! If there is no no snow
196      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
197      ! zftrice = io.qsr_ice       is below the surface
198      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
199      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
200      zhsu = 0.1_wp ! threshold for the computation of i0
201      DO ji = 1 , nidx
202         ! switches
203         isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
204         ! hs > 0, isnow = 1
205         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )     
206
207         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
208      END DO
209
210      !-------------------------------------------------------
211      ! Solar radiation absorbed / transmitted at the surface
212      ! Derivative of the non solar flux
213      !-------------------------------------------------------
214      DO ji = 1 , nidx
215         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
216         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
217         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
218         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value
219      END DO
220
221      !---------------------------------------------------------
222      ! Transmission - absorption of solar radiation in the ice
223      !---------------------------------------------------------
224
225      DO ji = 1, nidx           ! snow initialization
226         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
227      END DO
228
229      DO jk = 1, nlay_s          ! Radiation through snow
230         DO ji = 1, nidx
231            !                             ! radiation transmitted below the layer-th snow layer
232            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
233            !                             ! radiation absorbed by the layer-th snow layer
234            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
235         END DO
236      END DO
237
238      DO ji = 1, nidx           ! ice initialization
239         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
240      END DO
241
242      DO jk = 1, nlay_i          ! Radiation through ice
243         DO ji = 1, nidx
244            !                             ! radiation transmitted below the layer-th ice layer
245            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
246            !                             ! radiation absorbed by the layer-th ice layer
247            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
248         END DO
249      END DO
250
251      DO ji = 1, nidx           ! Radiation transmitted below the ice
252         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
253      END DO
254
255      !------------------------------------------------------------------------------|
256      !  3) Iterative procedure begins                                               |
257      !------------------------------------------------------------------------------|
258      !
259      DO ji = 1, nidx        ! Old surface temperature
260         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
261         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
262         t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary
263         zdti   (ji) =  1000._wp                                 ! initial value of error
264      END DO
265
266      DO jk = 1, nlay_s       ! Old snow temperature
267         DO ji = 1 , nidx
268            ztsb(ji,jk) =  t_s_1d(ji,jk)
269         END DO
270      END DO
271
272      DO jk = 1, nlay_i       ! Old ice temperature
273         DO ji = 1 , nidx
274            ztib(ji,jk) =  t_i_1d(ji,jk)
275         END DO
276      END DO
277
278      iconv    =  0           ! number of iterations
279      zdti_max =  1000._wp    ! maximal value of error on all points
280
281      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )
282         !
283         iconv = iconv + 1
284         !
285         !------------------------------------------------------------------------------|
286         ! 4) Sea ice thermal conductivity                                              |
287         !------------------------------------------------------------------------------|
288         !
289         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
290            DO ji = 1 , nidx
291               ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
292               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
293            END DO
294            DO jk = 1, nlay_i-1
295               DO ji = 1 , nidx
296                  ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
297                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0)
298                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
299               END DO
300            END DO
301
302         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
303            DO ji = 1 , nidx
304               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
305                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
306               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
307            END DO
308            DO jk = 1, nlay_i-1
309               DO ji = 1 , nidx
310                  ztcond_i(ji,jk) = rcdic +                                                                       & 
311                     &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              &
312                     &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   &
313                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 
314                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
315               END DO
316            END DO
317            DO ji = 1 , nidx
318               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
319                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
320               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
321            END DO
322         ENDIF
323         
324         !
325         !------------------------------------------------------------------------------|
326         !  5) G(he) - enhancement of thermal conductivity in mono-category case        |
327         !------------------------------------------------------------------------------|
328         !
329         ! Computation of effective thermal conductivity G(h)
330         ! Used in mono-category case only to simulate an ITD implicitly
331         ! Fichefet and Morales Maqueda, JGR 1997
332
333         zghe(:) = 1._wp
334
335         SELECT CASE ( nn_monocat )
336
337         CASE (1,3) ! LIM3
338
339            zepsilon = 0.1_wp
340            zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp
341
342            DO ji = 1, nidx
343   
344               ! Mean sea ice thermal conductivity
345               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp )
346
347               ! Effective thickness he (zhe)
348               zfac     = 1._wp / ( rn_cnd_s + zkimean )
349               zratio_s = rn_cnd_s   * zfac
350               zratio_i = zkimean * zfac
351               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
352
353               ! G(he)
354               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if >
355               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) )
356   
357               ! Impose G(he) < 2.
358               zghe(ji) = MIN( zghe(ji), 2._wp )
359
360            END DO
361
362         END SELECT
363
364         !
365         !------------------------------------------------------------------------------|
366         !  6) kappa factors                                                            |
367         !------------------------------------------------------------------------------|
368         !
369         !--- Snow
370         DO ji = 1, nidx
371            zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
372            zkappa_s(ji,0)        = zghe(ji) * rn_cnd_s * zfac
373            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cnd_s * zfac
374         END DO
375
376         DO jk = 1, nlay_s-1
377            DO ji = 1 , nidx
378               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cnd_s / MAX( epsi10, 2.0 * zh_s(ji) )
379            END DO
380         END DO
381
382         !--- Ice
383         DO jk = 1, nlay_i-1
384            DO ji = 1 , nidx
385               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) )
386            END DO
387         END DO
388
389         !--- Snow-ice interface
390         DO ji = 1 , nidx
391            zfac                  = 1./ MAX( epsi10 , zh_i(ji) )
392            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
393            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
394            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cnd_s * ztcond_i(ji,0) / & 
395           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cnd_s * zh_i(ji) ) )
396            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
397         END DO
398
399         !
400         !------------------------------------------------------------------------------|
401         ! 7) Sea ice specific heat, eta factors                                        |
402         !------------------------------------------------------------------------------|
403         !
404         DO jk = 1, nlay_i
405            DO ji = 1 , nidx
406               ztitemp(ji,jk)   = t_i_1d(ji,jk)
407               zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
408               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 )
409            END DO
410         END DO
411
412         DO jk = 1, nlay_s
413            DO ji = 1 , nidx
414               ztstemp(ji,jk) = t_s_1d(ji,jk)
415               zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 )
416            END DO
417         END DO
418
419         !
420         !------------------------------------------------------------------------------|
421         ! 8) surface flux computation                                                  |
422         !------------------------------------------------------------------------------|
423         !
424         IF ( ln_dqns_i ) THEN
425            DO ji = 1 , nidx
426               ! update of the non solar flux according to the update in T_su
427               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
428            END DO
429         ENDIF
430
431         ! Update incoming flux
432         DO ji = 1 , nidx
433            ! update incoming flux
434            zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation
435               &         + qns_ice_1d(ji)                ! non solar total flux (LWup, LWdw, SH, LH)
436         END DO
437
438         !
439         !------------------------------------------------------------------------------|
440         ! 9) tridiagonal system terms                                                  |
441         !------------------------------------------------------------------------------|
442         !
443         !!layer denotes the number of the layer in the snow or in the ice
444         !!numeq denotes the reference number of the equation in the tridiagonal
445         !!system, terms of tridiagonal system are indexed as following :
446         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
447
448         !!ice interior terms (top equation has the same form as the others)
449
450         DO numeq=1,nlay_i+3
451            DO ji = 1 , nidx
452               ztrid(ji,numeq,1) = 0.
453               ztrid(ji,numeq,2) = 0.
454               ztrid(ji,numeq,3) = 0.
455               zindterm(ji,numeq)= 0.
456               zindtbis(ji,numeq)= 0.
457               zdiagbis(ji,numeq)= 0.
458            ENDDO
459         ENDDO
460
461         DO numeq = nlay_s + 2, nlay_s + nlay_i 
462            DO ji = 1 , nidx
463               jk                 = numeq - nlay_s - 1
464               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
465               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
466               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
467               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
468            END DO
469         ENDDO
470
471         numeq =  nlay_s + nlay_i + 1
472         DO ji = 1 , nidx
473            !!ice bottom term
474            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
475            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
476            ztrid(ji,numeq,3)  =  0.0
477            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
478               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
479         ENDDO
480
481
482         DO ji = 1 , nidx
483            IF ( ht_s_1d(ji) > 0.0 ) THEN
484               !
485               !------------------------------------------------------------------------------|
486               !  snow-covered cells                                                          |
487               !------------------------------------------------------------------------------|
488               !
489               !!snow interior terms (bottom equation has the same form as the others)
490               DO numeq = 3, nlay_s + 1
491                  jk                  =  numeq - 1
492                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
493                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
494                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
495                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
496               END DO
497
498               !!case of only one layer in the ice (ice equation is altered)
499               IF ( nlay_i.eq.1 ) THEN
500                  ztrid(ji,nlay_s+2,3)    =  0.0
501                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
502               ENDIF
503
504               IF ( t_su_1d(ji) < rt0 ) THEN
505
506                  !------------------------------------------------------------------------------|
507                  !  case 1 : no surface melting - snow present                                  |
508                  !------------------------------------------------------------------------------|
509                  zdifcase(ji)    =  1.0
510                  numeqmin(ji)    =  1
511                  numeqmax(ji)    =  nlay_i + nlay_s + 1
512
513                  !!surface equation
514                  ztrid(ji,1,1)  = 0.0
515                  ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0)
516                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
517                  zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
518
519                  !!first layer of snow equation
520                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
521                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
522                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
523                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
524
525               ELSE 
526                  !
527                  !------------------------------------------------------------------------------|
528                  !  case 2 : surface is melting - snow present                                  |
529                  !------------------------------------------------------------------------------|
530                  !
531                  zdifcase(ji)    =  2.0
532                  numeqmin(ji)    =  2
533                  numeqmax(ji)    =  nlay_i + nlay_s + 1
534
535                  !!first layer of snow equation
536                  ztrid(ji,2,1)  =  0.0
537                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
538                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
539                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *  &
540                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
541               ENDIF
542            ELSE
543               !
544               !------------------------------------------------------------------------------|
545               !  cells without snow                                                          |
546               !------------------------------------------------------------------------------|
547               !
548               IF ( t_su_1d(ji) < rt0 ) THEN
549                  !
550                  !------------------------------------------------------------------------------|
551                  !  case 3 : no surface melting - no snow                                       |
552                  !------------------------------------------------------------------------------|
553                  !
554                  zdifcase(ji)      =  3.0
555                  numeqmin(ji)      =  nlay_s + 1
556                  numeqmax(ji)      =  nlay_i + nlay_s + 1
557
558                  !!surface equation   
559                  ztrid(ji,numeqmin(ji),1)   =  0.0
560                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1   
561                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
562                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
563
564                  !!first layer of ice equation
565                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
566                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
567                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1) 
568                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
569
570                  !!case of only one layer in the ice (surface & ice equations are altered)
571
572                  IF ( nlay_i == 1 ) THEN
573                     ztrid(ji,numeqmin(ji),1)    =  0.0
574                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0
575                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
576                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
577                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
578                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
579
580                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1) *  &
581                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
582                  ENDIF
583
584               ELSE
585
586                  !
587                  !------------------------------------------------------------------------------|
588                  ! case 4 : surface is melting - no snow                                        |
589                  !------------------------------------------------------------------------------|
590                  !
591                  zdifcase(ji)    =  4.0
592                  numeqmin(ji)    =  nlay_s + 2
593                  numeqmax(ji)    =  nlay_i + nlay_s + 1
594
595                  !!first layer of ice equation
596                  ztrid(ji,numeqmin(ji),1)      =  0.0
597                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
598                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
599                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1) *  &
600                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
601
602                  !!case of only one layer in the ice (surface & ice equations are altered)
603                  IF ( nlay_i == 1 ) THEN
604                     ztrid(ji,numeqmin(ji),1)  =  0.0
605                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
606                     ztrid(ji,numeqmin(ji),3)  =  0.0
607                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
608                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
609                  ENDIF
610
611               ENDIF
612            ENDIF
613
614         END DO
615
616         !
617         !------------------------------------------------------------------------------|
618         ! 10) tridiagonal system solving                                               |
619         !------------------------------------------------------------------------------|
620         !
621
622         ! Solve the tridiagonal system with Gauss elimination method.
623         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
624         ! McGraw-Hill 1984. 
625
626         maxnumeqmax = 0
627         minnumeqmin = nlay_i+5
628
629         DO ji = 1 , nidx
630            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
631            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
632            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
633            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
634         END DO
635
636         DO jk = minnumeqmin+1, maxnumeqmax
637            DO ji = 1 , nidx
638               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
639               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
640               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
641            END DO
642         END DO
643
644         DO ji = 1 , nidx
645            ! ice temperatures
646            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
647         END DO
648
649         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
650            DO ji = 1 , nidx
651               jk    =  numeq - nlay_s - 1
652               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
653            END DO
654         END DO
655
656         DO ji = 1 , nidx
657            ! snow temperatures     
658            IF (ht_s_1d(ji) > 0._wp) &
659               t_s_1d(ji,nlay_s)     =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
660               &                        / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) ) 
661
662            ! surface temperature
663            isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
664            ztsubit(ji) = t_su_1d(ji)
665            IF( t_su_1d(ji) < rt0 ) &
666               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
667               &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
668         END DO
669         !
670         !--------------------------------------------------------------------------
671         !  11) Has the scheme converged ?, end of the iterative procedure         |
672         !--------------------------------------------------------------------------
673         !
674         ! check that nowhere it has started to melt
675         ! zdti(ji) is a measure of error, it has to be under zdti_bnd
676         DO ji = 1 , nidx
677            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
678            zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
679         END DO
680
681         DO jk  =  1, nlay_s
682            DO ji = 1 , nidx
683               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  )
684               zdti  (ji)    = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
685            END DO
686         END DO
687
688         DO jk  =  1, nlay_i
689            DO ji = 1 , nidx
690               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
691               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
692               zdti  (ji)    =  MAX( zdti(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
693            END DO
694         END DO
695
696         ! Compute spatial maximum over all errors
697         ! note that this could be optimized substantially by iterating only the non-converging points
698         zdti_max = 0._wp
699         DO ji = 1, nidx
700            zdti_max = MAX( zdti_max, zdti(ji) )   
701         END DO
702         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
703
704      END DO  ! End of the do while iterative procedure
705
706      ! MV SIMIP 2016
707      !--- Snow-ice interfacial temperature (diagnostic SIMIP)
708      DO ji = 1, nidx
709         zfac        = 1. / MAX( epsi10 , rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
710         IF( zh_s(ji) >= 1.e-3 ) THEN
711            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) + &
712               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac
713         ELSE
714            t_si_1d(ji) = t_su_1d(ji)
715         ENDIF
716      END DO
717      ! END MV SIMIP 2016
718
719      IF( ln_icectl .AND. lwp ) THEN
720         WRITE(numout,*) ' zdti_max : ', zdti_max
721         WRITE(numout,*) ' iconv    : ', iconv
722      ENDIF
723
724      !
725      !-------------------------------------------------------------------------!
726      !   12) Fluxes at the interfaces                                          !
727      !-------------------------------------------------------------------------!
728      DO ji = 1, nidx
729         !                                ! surface ice conduction flux
730         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
731         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
732            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
733         !                                ! bottom ice conduction flux
734         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
735      END DO
736
737      ! MV SIMIP 2016
738      !--- Conduction fluxes (positive downwards)
739      DO ji = 1, nidx
740         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
741         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
742      END DO
743      ! END MV SIMIP 2016
744
745      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
746      CALL ice_thd_enmelt
747
748      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
749      IF ( ln_dqns_i ) THEN
750         DO ji = 1, nidx
751            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
752         END DO
753      END IF
754
755      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
756      DO ji = 1, nidx
757         zdq(ji)        = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
758            &                              SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
759         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
760            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
761         ELSE                          ! case T_su = 0degC
762            zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
763         ENDIF
764         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
765
766         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
767         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
768      END DO 
769
770      !-----------------------------------------
771      ! Heat flux used to warm/cool ice in W.m-2
772      !-----------------------------------------
773      DO ji = 1, nidx
774         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
775            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
776               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
777         ELSE                          ! case T_su = 0degC
778            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
779               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
780         ENDIF
781         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
782         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
783      END DO   
784      !
785   END SUBROUTINE ice_thd_dif
786
787
788   SUBROUTINE ice_thd_enmelt
789      !!-----------------------------------------------------------------------
790      !!                   ***  ROUTINE ice_thd_enmelt ***
791      !!                 
792      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
793      !!
794      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
795      !!-------------------------------------------------------------------
796      INTEGER  ::   ji, jk   ! dummy loop indices
797      REAL(wp) ::   ztmelts  ! local scalar
798      !!-------------------------------------------------------------------
799      !
800      DO jk = 1, nlay_i             ! Sea ice energy of melting
801         DO ji = 1, nidx
802            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0
803            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
804                                                          !   (sometimes dif scheme produces abnormally high temperatures)   
805            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           &
806               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   &
807               &                    - rcp  *         ( ztmelts-rt0 )  ) 
808         END DO
809      END DO
810      DO jk = 1, nlay_s             ! Snow energy of melting
811         DO ji = 1, nidx
812            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
813         END DO
814      END DO
815      !
816   END SUBROUTINE ice_thd_enmelt
817
818#else
819   !!----------------------------------------------------------------------
820   !!                   Dummy Module                 No LIM-3 sea-ice model
821   !!----------------------------------------------------------------------
822#endif
823
824   !!======================================================================
825END MODULE icethd_dif
Note: See TracBrowser for help on using the repository browser.