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

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

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