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.
limthd_dif.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 2760

Last change on this file since 2760 was 2715, checked in by rblod, 13 years ago

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