source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90 @ 8500

Last change on this file since 8500 was 8500, checked in by clem, 3 years ago

changes in style - part3 - move output into proper files and correct a bug in debug mode

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( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula
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         ENDIF
302
303         IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
304            DO ji = 1 , nidx
305               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
306                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
307               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
308            END DO
309            DO jk = 1, nlay_i-1
310               DO ji = 1 , nidx
311                  ztcond_i(ji,jk) = rcdic +                                                                       & 
312                     &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              &
313                     &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   &
314                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 
315                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
316               END DO
317            END DO
318            DO ji = 1 , nidx
319               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
320                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
321               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
322            END DO
323         ENDIF
324         
325         !
326         !------------------------------------------------------------------------------|
327         !  5) G(he) - enhancement of thermal conductivity in mono-category case        |
328         !------------------------------------------------------------------------------|
329         !
330         ! Computation of effective thermal conductivity G(h)
331         ! Used in mono-category case only to simulate an ITD implicitly
332         ! Fichefet and Morales Maqueda, JGR 1997
333
334         zghe(:) = 1._wp
335
336         SELECT CASE ( nn_monocat )
337
338         CASE (1,3) ! LIM3
339
340            zepsilon = 0.1_wp
341            zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp
342
343            DO ji = 1, nidx
344   
345               ! Mean sea ice thermal conductivity
346               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp )
347
348               ! Effective thickness he (zhe)
349               zfac     = 1._wp / ( rn_cdsn + zkimean )
350               zratio_s = rn_cdsn   * zfac
351               zratio_i = zkimean * zfac
352               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
353
354               ! G(he)
355               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if >
356               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) )
357   
358               ! Impose G(he) < 2.
359               zghe(ji) = MIN( zghe(ji), 2._wp )
360
361            END DO
362
363         END SELECT
364
365         !
366         !------------------------------------------------------------------------------|
367         !  6) kappa factors                                                            |
368         !------------------------------------------------------------------------------|
369         !
370         !--- Snow
371         DO ji = 1, nidx
372            zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
373            zkappa_s(ji,0)        = zghe(ji) * rn_cdsn * zfac
374            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cdsn * zfac
375         END DO
376
377         DO jk = 1, nlay_s-1
378            DO ji = 1 , nidx
379               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) )
380            END DO
381         END DO
382
383         !--- Ice
384         DO jk = 1, nlay_i-1
385            DO ji = 1 , nidx
386               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) )
387            END DO
388         END DO
389
390         !--- Snow-ice interface
391         DO ji = 1 , nidx
392            zfac                  = 1./ MAX( epsi10 , zh_i(ji) )
393            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
394            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
395            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / & 
396           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) )
397            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
398         END DO
399
400         !
401         !------------------------------------------------------------------------------|
402         ! 7) Sea ice specific heat, eta factors                                        |
403         !------------------------------------------------------------------------------|
404         !
405         DO jk = 1, nlay_i
406            DO ji = 1 , nidx
407               ztitemp(ji,jk)   = t_i_1d(ji,jk)
408               zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
409               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 )
410            END DO
411         END DO
412
413         DO jk = 1, nlay_s
414            DO ji = 1 , nidx
415               ztstemp(ji,jk) = t_s_1d(ji,jk)
416               zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 )
417            END DO
418         END DO
419
420         !
421         !------------------------------------------------------------------------------|
422         ! 8) surface flux computation                                                  |
423         !------------------------------------------------------------------------------|
424         !
425         IF ( ln_dqnsice ) THEN
426            DO ji = 1 , nidx
427               ! update of the non solar flux according to the update in T_su
428               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
429            END DO
430         ENDIF
431
432         ! Update incoming flux
433         DO ji = 1 , nidx
434            ! update incoming flux
435            zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation
436               &         + qns_ice_1d(ji)                ! non solar total flux (LWup, LWdw, SH, LH)
437         END DO
438
439         !
440         !------------------------------------------------------------------------------|
441         ! 9) tridiagonal system terms                                                  |
442         !------------------------------------------------------------------------------|
443         !
444         !!layer denotes the number of the layer in the snow or in the ice
445         !!numeq denotes the reference number of the equation in the tridiagonal
446         !!system, terms of tridiagonal system are indexed as following :
447         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
448
449         !!ice interior terms (top equation has the same form as the others)
450
451         DO numeq=1,nlay_i+3
452            DO ji = 1 , nidx
453               ztrid(ji,numeq,1) = 0.
454               ztrid(ji,numeq,2) = 0.
455               ztrid(ji,numeq,3) = 0.
456               zindterm(ji,numeq)= 0.
457               zindtbis(ji,numeq)= 0.
458               zdiagbis(ji,numeq)= 0.
459            ENDDO
460         ENDDO
461
462         DO numeq = nlay_s + 2, nlay_s + nlay_i 
463            DO ji = 1 , nidx
464               jk                 = numeq - nlay_s - 1
465               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
466               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
467               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
468               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
469            END DO
470         ENDDO
471
472         numeq =  nlay_s + nlay_i + 1
473         DO ji = 1 , nidx
474            !!ice bottom term
475            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
476            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
477            ztrid(ji,numeq,3)  =  0.0
478            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
479               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
480         ENDDO
481
482
483         DO ji = 1 , nidx
484            IF ( ht_s_1d(ji) > 0.0 ) THEN
485               !
486               !------------------------------------------------------------------------------|
487               !  snow-covered cells                                                          |
488               !------------------------------------------------------------------------------|
489               !
490               !!snow interior terms (bottom equation has the same form as the others)
491               DO numeq = 3, nlay_s + 1
492                  jk                  =  numeq - 1
493                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
494                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
495                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
496                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
497               END DO
498
499               !!case of only one layer in the ice (ice equation is altered)
500               IF ( nlay_i.eq.1 ) THEN
501                  ztrid(ji,nlay_s+2,3)    =  0.0
502                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
503               ENDIF
504
505               IF ( t_su_1d(ji) < rt0 ) THEN
506
507                  !------------------------------------------------------------------------------|
508                  !  case 1 : no surface melting - snow present                                  |
509                  !------------------------------------------------------------------------------|
510                  zdifcase(ji)    =  1.0
511                  numeqmin(ji)    =  1
512                  numeqmax(ji)    =  nlay_i + nlay_s + 1
513
514                  !!surface equation
515                  ztrid(ji,1,1)  = 0.0
516                  ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0)
517                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
518                  zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
519
520                  !!first layer of snow equation
521                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
522                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
523                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
524                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
525
526               ELSE 
527                  !
528                  !------------------------------------------------------------------------------|
529                  !  case 2 : surface is melting - snow present                                  |
530                  !------------------------------------------------------------------------------|
531                  !
532                  zdifcase(ji)    =  2.0
533                  numeqmin(ji)    =  2
534                  numeqmax(ji)    =  nlay_i + nlay_s + 1
535
536                  !!first layer of snow equation
537                  ztrid(ji,2,1)  =  0.0
538                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
539                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
540                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *  &
541                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
542               ENDIF
543            ELSE
544               !
545               !------------------------------------------------------------------------------|
546               !  cells without snow                                                          |
547               !------------------------------------------------------------------------------|
548               !
549               IF ( t_su_1d(ji) < rt0 ) THEN
550                  !
551                  !------------------------------------------------------------------------------|
552                  !  case 3 : no surface melting - no snow                                       |
553                  !------------------------------------------------------------------------------|
554                  !
555                  zdifcase(ji)      =  3.0
556                  numeqmin(ji)      =  nlay_s + 1
557                  numeqmax(ji)      =  nlay_i + nlay_s + 1
558
559                  !!surface equation   
560                  ztrid(ji,numeqmin(ji),1)   =  0.0
561                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1   
562                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
563                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
564
565                  !!first layer of ice equation
566                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
567                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
568                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1) 
569                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
570
571                  !!case of only one layer in the ice (surface & ice equations are altered)
572
573                  IF ( nlay_i == 1 ) THEN
574                     ztrid(ji,numeqmin(ji),1)    =  0.0
575                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0
576                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
577                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
578                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
579                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
580
581                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1) *  &
582                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
583                  ENDIF
584
585               ELSE
586
587                  !
588                  !------------------------------------------------------------------------------|
589                  ! case 4 : surface is melting - no snow                                        |
590                  !------------------------------------------------------------------------------|
591                  !
592                  zdifcase(ji)    =  4.0
593                  numeqmin(ji)    =  nlay_s + 2
594                  numeqmax(ji)    =  nlay_i + nlay_s + 1
595
596                  !!first layer of ice equation
597                  ztrid(ji,numeqmin(ji),1)      =  0.0
598                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
599                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
600                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1) *  &
601                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
602
603                  !!case of only one layer in the ice (surface & ice equations are altered)
604                  IF ( nlay_i == 1 ) THEN
605                     ztrid(ji,numeqmin(ji),1)  =  0.0
606                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
607                     ztrid(ji,numeqmin(ji),3)  =  0.0
608                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
609                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
610                  ENDIF
611
612               ENDIF
613            ENDIF
614
615         END DO
616
617         !
618         !------------------------------------------------------------------------------|
619         ! 10) tridiagonal system solving                                               |
620         !------------------------------------------------------------------------------|
621         !
622
623         ! Solve the tridiagonal system with Gauss elimination method.
624         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
625         ! McGraw-Hill 1984. 
626
627         maxnumeqmax = 0
628         minnumeqmin = nlay_i+5
629
630         DO ji = 1 , nidx
631            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
632            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
633            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
634            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
635         END DO
636
637         DO jk = minnumeqmin+1, maxnumeqmax
638            DO ji = 1 , nidx
639               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
640               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
641               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
642            END DO
643         END DO
644
645         DO ji = 1 , nidx
646            ! ice temperatures
647            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
648         END DO
649
650         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
651            DO ji = 1 , nidx
652               jk    =  numeq - nlay_s - 1
653               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
654            END DO
655         END DO
656
657         DO ji = 1 , nidx
658            ! snow temperatures     
659            IF (ht_s_1d(ji) > 0._wp) &
660               t_s_1d(ji,nlay_s)     =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
661               &                        / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) ) 
662
663            ! surface temperature
664            isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
665            ztsubit(ji) = t_su_1d(ji)
666            IF( t_su_1d(ji) < rt0 ) &
667               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
668               &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
669         END DO
670         !
671         !--------------------------------------------------------------------------
672         !  11) Has the scheme converged ?, end of the iterative procedure         |
673         !--------------------------------------------------------------------------
674         !
675         ! check that nowhere it has started to melt
676         ! zdti(ji) is a measure of error, it has to be under zdti_bnd
677         DO ji = 1 , nidx
678            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
679            zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
680         END DO
681
682         DO jk  =  1, nlay_s
683            DO ji = 1 , nidx
684               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  )
685               zdti  (ji)    = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
686            END DO
687         END DO
688
689         DO jk  =  1, nlay_i
690            DO ji = 1 , nidx
691               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
692               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
693               zdti  (ji)    =  MAX( zdti(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
694            END DO
695         END DO
696
697         ! Compute spatial maximum over all errors
698         ! note that this could be optimized substantially by iterating only the non-converging points
699         zdti_max = 0._wp
700         DO ji = 1, nidx
701            zdti_max = MAX( zdti_max, zdti(ji) )   
702         END DO
703         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
704
705      END DO  ! End of the do while iterative procedure
706
707      ! MV SIMIP 2016
708      !--- Snow-ice interfacial temperature (diagnostic SIMIP)
709      DO ji = 1, nidx
710         zfac        = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
711         IF( zh_s(ji) >= 1.e-3 ) THEN
712            t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + &
713               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac
714         ELSE
715            t_si_1d(ji) = t_su_1d(ji)
716         ENDIF
717      END DO
718      ! END MV SIMIP 2016
719
720      IF( ln_limctl .AND. lwp ) THEN
721         WRITE(numout,*) ' zdti_max : ', zdti_max
722         WRITE(numout,*) ' iconv    : ', iconv
723      ENDIF
724
725      !
726      !-------------------------------------------------------------------------!
727      !   12) Fluxes at the interfaces                                          !
728      !-------------------------------------------------------------------------!
729      DO ji = 1, nidx
730         !                                ! surface ice conduction flux
731         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
732         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
733            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
734         !                                ! bottom ice conduction flux
735         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
736      END DO
737
738      ! MV SIMIP 2016
739      !--- Conduction fluxes (positive downwards)
740      DO ji = 1, nidx
741         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
742         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
743      END DO
744      ! END MV SIMIP 2016
745
746      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
747      CALL ice_thd_enmelt
748
749      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
750      IF ( ln_dqnsice ) THEN
751         DO ji = 1, nidx
752            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
753         END DO
754      END IF
755
756      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
757      DO ji = 1, nidx
758         zdq(ji)        = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
759            &                              SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
760         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
761            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
762         ELSE                          ! case T_su = 0degC
763            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 
764         ENDIF
765         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
766
767         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
768         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
769      END DO 
770
771      !-----------------------------------------
772      ! Heat flux used to warm/cool ice in W.m-2
773      !-----------------------------------------
774      DO ji = 1, nidx
775         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
776            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
777               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
778         ELSE                          ! case T_su = 0degC
779            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
780               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
781         ENDIF
782         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
783         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
784      END DO   
785      !
786   END SUBROUTINE ice_thd_dif
787
788
789   SUBROUTINE ice_thd_enmelt
790      !!-----------------------------------------------------------------------
791      !!                   ***  ROUTINE ice_thd_enmelt ***
792      !!                 
793      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
794      !!
795      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
796      !!-------------------------------------------------------------------
797      INTEGER  ::   ji, jk   ! dummy loop indices
798      REAL(wp) ::   ztmelts  ! local scalar
799      !!-------------------------------------------------------------------
800      !
801      DO jk = 1, nlay_i             ! Sea ice energy of melting
802         DO ji = 1, nidx
803            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0
804            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
805                                                          !   (sometimes dif scheme produces abnormally high temperatures)   
806            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           &
807               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   &
808               &                    - rcp  *         ( ztmelts-rt0 )  ) 
809         END DO
810      END DO
811      DO jk = 1, nlay_s             ! Snow energy of melting
812         DO ji = 1, nidx
813            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
814         END DO
815      END DO
816      !
817   END SUBROUTINE ice_thd_enmelt
818
819#else
820   !!----------------------------------------------------------------------
821   !!                   Dummy Module                 No LIM-3 sea-ice model
822   !!----------------------------------------------------------------------
823#endif
824
825   !!======================================================================
826END MODULE icethd_dif
Note: See TracBrowser for help on using the repository browser.