source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 4205

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