source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 4626

Last change on this file since 4626 was 4626, checked in by gm, 7 years ago

dev_r4028_CNRS_LIM3_MV2014 : minor corrections on LIM3

  • Property svn:keywords set to Id
File size: 37.1 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
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   lim_thd_dif   ! called by lim_thd
32
33   REAL(wp) ::   epsi10      =  1.e-10_wp    !
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 , jl )
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_b
77      !!           ice/snow temperatures   : t_i_b, t_s_b
78      !!           ice salinities          : s_i_b
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_b, ht_s_b
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   ! Start point on which the  the computation is applied
94      INTEGER , INTENT (in) ::   kiut    ! End point on which the  the computation is applied
95      INTEGER , INTENT (in) ::   jl      ! Category number
96
97      !! * Local variables
98      INTEGER ::   ji          ! spatial loop index
99      INTEGER ::   ii, ij      ! temporary dummy loop index
100      INTEGER ::   numeq       ! current reference number of equation
101      INTEGER ::   layer       ! vertical dummy loop index
102      INTEGER ::   nconv       ! number of iterations in iterative procedure
103      INTEGER ::   minnumeqmin, maxnumeqmax
104      INTEGER, DIMENSION(kiut) ::   numeqmin   ! reference number of top equation
105      INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation
106      INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow
107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
108      REAL(wp) ::   zg1       =  2._wp        !
109      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
110      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
111      REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow
112      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
113      REAL(wp) ::   ztmelt_i    ! ice melting temperature
114      REAL(wp) ::   zerritmax   ! current maximal error on temperature
115      REAL(wp), DIMENSION(kiut) ::   ztfs        ! ice melting point
116      REAL(wp), DIMENSION(kiut) ::   ztsuold     ! old surface temperature (before the iterative procedure )
117      REAL(wp), DIMENSION(kiut) ::   ztsuoldit   ! surface temperature at previous iteration
118      REAL(wp), DIMENSION(kiut) ::   zh_i        ! ice layer thickness
119      REAL(wp), DIMENSION(kiut) ::   zh_s        ! snow layer thickness
120      REAL(wp), DIMENSION(kiut) ::   zfsw        ! solar radiation absorbed at the surface
121      REAL(wp), DIMENSION(kiut) ::   zf          ! surface flux function
122      REAL(wp), DIMENSION(kiut) ::   dzf         ! derivative of the surface flux function
123      REAL(wp), DIMENSION(kiut) ::   zerrit      ! current error on temperature
124      REAL(wp), DIMENSION(kiut) ::   zdifcase    ! case of the equation resolution (1->4)
125      REAL(wp), DIMENSION(kiut) ::   zftrice     ! solar radiation transmitted through the ice
126      REAL(wp), DIMENSION(kiut) ::   zihic, zhsu
127      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity
128      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice
129      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice
130      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice
131      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztiold      ! Old temperature in the ice
132      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zeta_i      ! Eta factor in the ice
133      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
134      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zspeche_i   ! Ice specific heat
135      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   z_i         ! Vertical cotes of the layers in the ice
136      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow
137      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow
138      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow
139      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zeta_s       ! Eta factor in the snow
140      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztstemp      ! Temporary temperature in the snow to check the convergence
141      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztsold       ! Temporary temperature in the snow
142      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   z_s          ! Vertical cotes of the layers in the snow
143      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindterm   ! Independent term
144      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindtbis   ! temporary independent term
145      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis
146      REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms
147      REAL(wp) ::   ztemp   ! local scalar
148      !!------------------------------------------------------------------     
149      !
150      !------------------------------------------------------------------------------!
151      ! 1) Initialization                                                            !
152      !------------------------------------------------------------------------------!
153      !
154      DO ji = kideb , kiut
155         ! is there snow or not
156         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  )
157         ! surface temperature of fusion
158!!gm ???  ztfs(ji) = rtt !!!????
159         ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt
160         ! layer thickness
161         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
162         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
163      END DO
164
165      !--------------------
166      ! Ice / snow layers
167      !--------------------
168
169      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
170      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
171
172      DO layer = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
173         DO ji = kideb , kiut
174            z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )
175         END DO
176      END DO
177
178      DO layer = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
179         DO ji = kideb , kiut
180            z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )
181         END DO
182      END DO
183      !
184      !------------------------------------------------------------------------------|
185      ! 2) Radiations                                                                |
186      !------------------------------------------------------------------------------|
187      !
188      !-------------------
189      ! Computation of i0
190      !-------------------
191      ! i0 describes the fraction of solar radiation which does not contribute
192      ! to the surface energy budget but rather penetrates inside the ice.
193      ! We assume that no radiation is transmitted through the snow
194      ! If there is no no snow
195      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
196      ! zftrice = io.qsr_ice       is below the surface
197      ! fstbif  = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
198
199      DO ji = kideb , kiut
200         ! switches
201         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  ) 
202         ! hs > 0, isnow = 1
203         zhsu (ji) = hnzst  ! threshold for the computation of i0
204         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )     
205
206         i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
207         !fr1_i0_1d = i0 for a thin ice surface
208         !fr1_i0_2d = i0 for a thick ice surface
209         !            a function of the cloud cover
210         !
211         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0)
212         !formula used in Cice
213      END DO
214
215      !-------------------------------------------------------
216      ! Solar radiation absorbed / transmitted at the surface
217      ! Derivative of the non solar flux
218      !-------------------------------------------------------
219      DO ji = kideb , kiut
220         zfsw   (ji) =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
221         zftrice(ji) =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
222         dzf    (ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
223      END DO
224
225      !---------------------------------------------------------
226      ! Transmission - absorption of solar radiation in the ice
227      !---------------------------------------------------------
228
229      DO ji = kideb, kiut           ! snow initialization
230         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
231      END DO
232
233      DO layer = 1, nlay_s          ! Radiation through snow
234         DO ji = kideb, kiut
235            !                             ! radiation transmitted below the layer-th snow layer
236            zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )
237            !                             ! radiation absorbed by the layer-th snow layer
238            zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
239         END DO
240      END DO
241
242      DO ji = kideb, kiut           ! ice initialization
243         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) )
244      END DO
245
246      DO layer = 1, nlay_i          ! Radiation through ice
247         DO ji = kideb, kiut
248            !                             ! radiation transmitted below the layer-th ice layer
249            zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )
250            !                             ! radiation absorbed by the layer-th ice layer
251            zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
252         END DO
253      END DO
254
255      DO ji = kideb, kiut           ! Radiation transmitted below the ice
256         fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif
257      END DO
258
259      ! +++++
260      ! just to check energy conservation
261      DO ji = kideb, kiut
262         ii = MOD( npb(ji) - 1 , jpi ) + 1
263         ij =    ( npb(ji) - 1 ) / jpi + 1
264         fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif
265      END DO
266      ! +++++
267
268      DO layer = 1, nlay_i
269         DO ji = kideb, kiut
270            radab(ji,layer) = zradab_i(ji,layer)
271         END DO
272      END DO
273
274      !
275      !------------------------------------------------------------------------------|
276      !  3) Iterative procedure begins                                               |
277      !------------------------------------------------------------------------------|
278      !
279      DO ji = kideb, kiut        ! Old surface temperature
280         ztsuold  (ji) =  t_su_b(ji)                              ! temperature at the beg of iter pr.
281         ztsuoldit(ji) =  t_su_b(ji)                              ! temperature at the previous iter
282         t_su_b   (ji) =  MIN( t_su_b(ji), ztfs(ji)-0.00001 )     ! necessary
283         zerrit   (ji) =  1000._wp                                ! initial value of error
284      END DO
285
286      DO layer = 1, nlay_s       ! Old snow temperature
287         DO ji = kideb , kiut
288            ztsold(ji,layer) =  t_s_b(ji,layer)
289         END DO
290      END DO
291
292      DO layer = 1, nlay_i       ! Old ice temperature
293         DO ji = kideb , kiut
294            ztiold(ji,layer) =  t_i_b(ji,layer)
295         END DO
296      END DO
297
298      nconv     =  0           ! number of iterations
299      zerritmax =  1000._wp    ! maximal value of error on all points
300
301      DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd )
302         !
303         nconv = nconv + 1
304         !
305         !------------------------------------------------------------------------------|
306         ! 4) Sea ice thermal conductivity                                              |
307         !------------------------------------------------------------------------------|
308         !
309         IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula
310            DO ji = kideb , kiut
311               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)
312               ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin)
313            END DO
314            DO layer = 1, nlay_i-1
315               DO ji = kideb , kiut
316                  ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  &
317                     MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
318                  ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
319               END DO
320            END DO
321         ENDIF
322
323         IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
324            DO ji = kideb , kiut
325               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   &
326                  &                   - 0.011_wp * ( t_i_b(ji,1) - rtt ) 
327               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
328            END DO
329            DO layer = 1, nlay_i-1
330               DO ji = kideb , kiut
331                  ztemp = t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2._wp * rtt
332                  ztcond_i(ji,layer) = rcdic + 0.0900_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
333                     &                                   / MIN( -2.0_wp * epsi10, ztemp )   &
334                     &                       - 0.0055_wp * ztemp
335                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
336               END DO
337            END DO
338            DO ji = kideb , kiut
339               ztemp = t_bo_b(ji) - rtt
340               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN( -epsi10, ztemp )   &
341                  &                        - 0.011_wp * ztemp 
342               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
343            END DO
344         ENDIF
345         !
346         !------------------------------------------------------------------------------|
347         !  5) kappa factors                                                            |
348         !------------------------------------------------------------------------------|
349         !
350         DO ji = kideb, kiut
351
352            !-- Snow kappa factors
353            zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji))
354            zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji))
355         END DO
356
357         DO layer = 1, nlay_s-1
358            DO ji = kideb , kiut
359               zkappa_s(ji,layer)  = 2.0 * rcdsn / &
360                  MAX(epsi10,2.0*zh_s(ji))
361            END DO
362         END DO
363
364         DO layer = 1, nlay_i-1
365            DO ji = kideb , kiut
366               !-- Ice kappa factors
367               zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ &
368                  MAX(epsi10,2.0*zh_i(ji)) 
369            END DO
370         END DO
371
372         DO ji = kideb , kiut
373            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))
374            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))
375            !-- Interface
376            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, &
377               (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji)))
378            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) &
379               + zkappa_i(ji,0)*REAL( 1 - isnow(ji) )
380         END DO
381         !
382         !------------------------------------------------------------------------------|
383         ! 6) Sea ice specific heat, eta factors                                        |
384         !------------------------------------------------------------------------------|
385         !
386         DO layer = 1, nlay_i
387            DO ji = kideb , kiut
388               ztitemp(ji,layer)   = t_i_b(ji,layer)
389               zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
390                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)
391               zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &
392                  epsi10)
393            END DO
394         END DO
395
396         DO layer = 1, nlay_s
397            DO ji = kideb , kiut
398               ztstemp(ji,layer) = t_s_b(ji,layer)
399               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)
400            END DO
401         END DO
402         !
403         !------------------------------------------------------------------------------|
404         ! 7) surface flux computation                                                  |
405         !------------------------------------------------------------------------------|
406         !
407         DO ji = kideb , kiut
408
409            ! update of the non solar flux according to the update in T_su
410            qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * & 
411               ( t_su_b(ji) - ztsuoldit(ji) )
412
413            ! update incoming flux
414            zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation
415               + qnsr_ice_1d(ji)           ! non solar total flux
416            ! (LWup, LWdw, SH, LH)
417
418         END DO
419
420         !
421         !------------------------------------------------------------------------------|
422         ! 8) tridiagonal system terms                                                  |
423         !------------------------------------------------------------------------------|
424         !
425         !!layer denotes the number of the layer in the snow or in the ice
426         !!numeq denotes the reference number of the equation in the tridiagonal
427         !!system, terms of tridiagonal system are indexed as following :
428         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
429
430         !!ice interior terms (top equation has the same form as the others)
431
432         DO numeq=1,jkmax+2
433            DO ji = kideb , kiut
434               ztrid(ji,numeq,1) = 0.
435               ztrid(ji,numeq,2) = 0.
436               ztrid(ji,numeq,3) = 0.
437               zindterm(ji,numeq)= 0.
438               zindtbis(ji,numeq)= 0.
439               zdiagbis(ji,numeq)= 0.
440            ENDDO
441         ENDDO
442
443         DO numeq = nlay_s + 2, nlay_s + nlay_i 
444            DO ji = kideb , kiut
445               layer              = numeq - nlay_s - 1
446               ztrid(ji,numeq,1)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
447               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
448                  zkappa_i(ji,layer))
449               ztrid(ji,numeq,3)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer)
450               zindterm(ji,numeq) =  ztiold(ji,layer) + zeta_i(ji,layer)* &
451                  zradab_i(ji,layer)
452            END DO
453         ENDDO
454
455         numeq =  nlay_s + nlay_i + 1
456         DO ji = kideb , kiut
457            !!ice bottom term
458            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
459            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 &
460               +  zkappa_i(ji,nlay_i-1) )
461            ztrid(ji,numeq,3)  =  0.0
462            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
463               ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
464               *  t_bo_b(ji) ) 
465         ENDDO
466
467
468         DO ji = kideb , kiut
469            IF ( ht_s_b(ji).gt.0.0 ) THEN
470               !
471               !------------------------------------------------------------------------------|
472               !  snow-covered cells                                                          |
473               !------------------------------------------------------------------------------|
474               !
475               !!snow interior terms (bottom equation has the same form as the others)
476               DO numeq = 3, nlay_s + 1
477                  layer =  numeq - 1
478                  ztrid(ji,numeq,1)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
479                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
480                     zkappa_s(ji,layer) )
481                  ztrid(ji,numeq,3)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer)
482                  zindterm(ji,numeq)  =  ztsold(ji,layer) + zeta_s(ji,layer)* &
483                     zradab_s(ji,layer)
484               END DO
485
486               !!case of only one layer in the ice (ice equation is altered)
487               IF ( nlay_i.eq.1 ) THEN
488                  ztrid(ji,nlay_s+2,3)    =  0.0
489                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
490                     t_bo_b(ji) 
491               ENDIF
492
493               IF ( t_su_b(ji) .LT. rtt ) THEN
494
495                  !------------------------------------------------------------------------------|
496                  !  case 1 : no surface melting - snow present                                  |
497                  !------------------------------------------------------------------------------|
498                  zdifcase(ji)    =  1.0
499                  numeqmin(ji)    =  1
500                  numeqmax(ji)    =  nlay_i + nlay_s + 1
501
502                  !!surface equation
503                  ztrid(ji,1,1) = 0.0
504                  ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
505                  ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
506                  zindterm(ji,1) = dzf(ji)*t_su_b(ji)   - zf(ji)
507
508                  !!first layer of snow equation
509                  ztrid(ji,2,1)  =  - zkappa_s(ji,0)*zg1s*zeta_s(ji,1)
510                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s)
511                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
512                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
513
514               ELSE 
515                  !
516                  !------------------------------------------------------------------------------|
517                  !  case 2 : surface is melting - snow present                                  |
518                  !------------------------------------------------------------------------------|
519                  !
520                  zdifcase(ji)    =  2.0
521                  numeqmin(ji)    =  2
522                  numeqmax(ji)    =  nlay_i + nlay_s + 1
523
524                  !!first layer of snow equation
525                  ztrid(ji,2,1)  =  0.0
526                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + &
527                     zkappa_s(ji,0) * zg1s )
528                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
529                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *            &
530                     ( zradab_s(ji,1) +                         &
531                     zkappa_s(ji,0) * zg1s * t_su_b(ji) ) 
532               ENDIF
533            ELSE
534               !
535               !------------------------------------------------------------------------------|
536               !  cells without snow                                                          |
537               !------------------------------------------------------------------------------|
538               !
539               IF (t_su_b(ji) .LT. rtt) THEN
540                  !
541                  !------------------------------------------------------------------------------|
542                  !  case 3 : no surface melting - no snow                                       |
543                  !------------------------------------------------------------------------------|
544                  !
545                  zdifcase(ji)      =  3.0
546                  numeqmin(ji)      =  nlay_s + 1
547                  numeqmax(ji)      =  nlay_i + nlay_s + 1
548
549                  !!surface equation   
550                  ztrid(ji,numeqmin(ji),1)   =  0.0
551                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1   
552                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
553                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_b(ji) - zf(ji)
554
555                  !!first layer of ice equation
556                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
557                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) & 
558                     + zkappa_i(ji,0) * zg1 )
559                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1) 
560                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 
561
562                  !!case of only one layer in the ice (surface & ice equations are altered)
563
564                  IF (nlay_i.eq.1) THEN
565                     ztrid(ji,numeqmin(ji),1)    =  0.0
566                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0)*2.0
567                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0)*2.0
568                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0)*2.0*zeta_i(ji,1)
569                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + &
570                        zkappa_i(ji,1))
571                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
572
573                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1)* &
574                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
575                  ENDIF
576
577               ELSE
578
579                  !
580                  !------------------------------------------------------------------------------|
581                  ! case 4 : surface is melting - no snow                                        |
582                  !------------------------------------------------------------------------------|
583                  !
584                  zdifcase(ji)    =  4.0
585                  numeqmin(ji)    =  nlay_s + 2
586                  numeqmax(ji)    =  nlay_i + nlay_s + 1
587
588                  !!first layer of ice equation
589                  ztrid(ji,numeqmin(ji),1)      =  0.0
590                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* &
591                     zg1) 
592                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
593                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
594                     zkappa_i(ji,0) * zg1 * t_su_b(ji) ) 
595
596                  !!case of only one layer in the ice (surface & ice equations are altered)
597                  IF (nlay_i.eq.1) THEN
598                     ztrid(ji,numeqmin(ji),1)  =  0.0
599                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + &
600                        zkappa_i(ji,1))
601                     ztrid(ji,numeqmin(ji),3)  =  0.0
602                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1)* &
603                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
604                        + t_su_b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0
605                  ENDIF
606
607               ENDIF
608            ENDIF
609
610         END DO
611
612         !
613         !------------------------------------------------------------------------------|
614         ! 9) tridiagonal system solving                                                |
615         !------------------------------------------------------------------------------|
616         !
617
618         ! Solve the tridiagonal system with Gauss elimination method.
619         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
620         ! McGraw-Hill 1984. 
621
622         maxnumeqmax = 0
623         minnumeqmin = jkmax+4
624
625         DO ji = kideb , kiut
626            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
627            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
628            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
629            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
630         END DO
631
632         DO layer = minnumeqmin+1, maxnumeqmax
633            DO ji = kideb , kiut
634               numeq               =  min(max(numeqmin(ji)+1,layer),numeqmax(ji))
635               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
636                  ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
637               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1)* &
638                  zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
639            END DO
640         END DO
641
642         DO ji = kideb , kiut
643            ! ice temperatures
644            t_i_b(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
645         END DO
646
647         DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
648            DO ji = kideb , kiut
649               layer    =  numeq - nlay_s - 1
650               t_i_b(ji,layer)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
651                  t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
652            END DO
653         END DO
654
655         DO ji = kideb , kiut
656            ! snow temperatures     
657            IF (ht_s_b(ji).GT.0._wp) &
658               t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
659               *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
660               *        MAX(0.0,SIGN(1.0,ht_s_b(ji))) 
661
662            ! surface temperature
663            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  )
664            ztsuoldit(ji) = t_su_b(ji)
665            IF( t_su_b(ji) < ztfs(ji) ) &
666               t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1)   &
667               &          + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
668         END DO
669         !
670         !--------------------------------------------------------------------------
671         !  10) Has the scheme converged ?, end of the iterative procedure         |
672         !--------------------------------------------------------------------------
673         !
674         ! check that nowhere it has started to melt
675         ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
676         DO ji = kideb , kiut
677            t_su_b(ji) =  MAX(  MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp  )
678            zerrit(ji) =  ABS( t_su_b(ji) - ztsuoldit(ji) )     
679         END DO
680
681         DO layer  =  1, nlay_s
682            DO ji = kideb , kiut
683               ii = MOD( npb(ji) - 1, jpi ) + 1
684               ij = ( npb(ji) - 1 ) / jpi + 1
685               t_s_b(ji,layer) = MAX(  MIN( t_s_b(ji,layer), rtt ), 190._wp  )
686               zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
687            END DO
688         END DO
689
690         DO layer  =  1, nlay_i
691            DO ji = kideb , kiut
692               ztmelt_i        = -tmut * s_i_b(ji,layer) + rtt 
693               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
694               zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
695            END DO
696         END DO
697
698         ! Compute spatial maximum over all errors
699         ! note that this could be optimized substantially by iterating only the non-converging points
700         zerritmax = 0._wp
701         DO ji = kideb, kiut
702            zerritmax = MAX( zerritmax, zerrit(ji) )   
703         END DO
704         IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
705
706      END DO  ! End of the do while iterative procedure
707
708      IF( ln_nicep .AND. lwp ) THEN
709         WRITE(numout,*) ' zerritmax : ', zerritmax
710         WRITE(numout,*) ' nconv     : ', nconv
711      ENDIF
712
713      !
714      !-------------------------------------------------------------------------!
715      !   11) Fluxes at the interfaces                                          !
716      !-------------------------------------------------------------------------!
717      DO ji = kideb, kiut
718#if ! defined key_coupled
719         ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
720         qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) )
721#endif
722         !                                ! surface ice conduction flux
723         isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) )  )
724         fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji))   &
725            &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_b(ji,1) - t_su_b(ji))
726         !                                ! bottom ice conduction flux
727         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
728      END DO
729
730      !-------------------------!
731      ! Heat conservation       !
732      !-------------------------!
733      IF( con_i .AND. jiindex_1d > 0 ) THEN
734         DO ji = kideb, kiut
735            ! Upper snow value
736            fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 
737            ! Bott. snow value
738            fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 
739         END DO
740         DO ji = kideb, kiut         ! Upper ice layer
741            fc_i(ji,0) = - REAL( isnow(ji) ) * &  ! interface flux if there is snow
742               ( zkappa_i(ji,0)  * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) &
743               - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * & 
744               zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not
745         END DO
746         DO layer = 1, nlay_i - 1         ! Internal ice layers
747            DO ji = kideb, kiut
748               fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) )
749               ii = MOD( npb(ji) - 1, jpi ) + 1
750               ij = ( npb(ji) - 1 ) / jpi + 1
751            END DO
752         END DO
753         DO ji = kideb, kiut         ! Bottom ice layers
754            fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
755         END DO
756      ENDIF
757      !
758   END SUBROUTINE lim_thd_dif
759
760#else
761   !!----------------------------------------------------------------------
762   !!                   Dummy Module                 No LIM-3 sea-ice model
763   !!----------------------------------------------------------------------
764CONTAINS
765   SUBROUTINE lim_thd_dif          ! Empty routine
766   END SUBROUTINE lim_thd_dif
767#endif
768   !!======================================================================
769END MODULE limthd_dif
Note: See TracBrowser for help on using the repository browser.