source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 4333

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

remove remaining bugs in LIM3, so that it can run in both regional and global config

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