source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5047

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

LIM3 cleaning (1): namelist

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