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 branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 3610

Last change on this file since 3610 was 3610, checked in by acc, 10 years ago

Branch dev_NOC_2012_r3555. #1006. Step 5: Merge in trunk changes between revision 3337 and 3385

  • Property svn:keywords set to Id
File size: 36.8 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   USE wrk_nemo         ! work arrays
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   lim_thd_dif   ! called by lim_thd
30
31   REAL(wp) ::   epsi20 = 1e-20     ! constant values
32   REAL(wp) ::   epsi13 = 1e-13     ! constant values
33
34   !!----------------------------------------------------------------------
35   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE lim_thd_dif( kideb , kiut , 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) ::   zeps      =  1.e-10_wp    !
108      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
109      REAL(wp) ::   zg1       =  2._wp        !
110      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
111      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
112      REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow
113      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
114      REAL(wp) ::   zht_smin  =  1.e-4_wp     ! minimum snow depth
115      REAL(wp) ::   ztmelt_i    ! ice melting temperature
116      REAL(wp) ::   zerritmax   ! current maximal error on temperature
117      REAL(wp), DIMENSION(kiut) ::   ztfs        ! ice melting point
118      REAL(wp), DIMENSION(kiut) ::   ztsuold     ! old surface temperature (before the iterative procedure )
119      REAL(wp), DIMENSION(kiut) ::   ztsuoldit   ! surface temperature at previous iteration
120      REAL(wp), DIMENSION(kiut) ::   zh_i        ! ice layer thickness
121      REAL(wp), DIMENSION(kiut) ::   zh_s        ! snow layer thickness
122      REAL(wp), DIMENSION(kiut) ::   zfsw        ! solar radiation absorbed at the surface
123      REAL(wp), DIMENSION(kiut) ::   zf          ! surface flux function
124      REAL(wp), DIMENSION(kiut) ::   dzf         ! derivative of the surface flux function
125      REAL(wp), DIMENSION(kiut) ::   zerrit      ! current error on temperature
126      REAL(wp), DIMENSION(kiut) ::   zdifcase    ! case of the equation resolution (1->4)
127      REAL(wp), DIMENSION(kiut) ::   zftrice     ! solar radiation transmitted through the ice
128      REAL(wp), DIMENSION(kiut) ::   zihic, zhsu
129      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity
130      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice
131      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice
132      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice
133      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztiold      ! Old temperature in the ice
134      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zeta_i      ! Eta factor in the ice
135      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
136      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zspeche_i   ! Ice specific heat
137      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   z_i         ! Vertical cotes of the layers in the ice
138      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow
139      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow
140      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow
141      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zeta_s       ! Eta factor in the snow
142      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztstemp      ! Temporary temperature in the snow to check the convergence
143      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztsold       ! Temporary temperature in the snow
144      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   z_s          ! Vertical cotes of the layers in the snow
145      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindterm   ! Independent term
146      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindtbis   ! temporary independent term
147      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis
148      REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms
149      !!------------------------------------------------------------------
150     
151      !
152      !------------------------------------------------------------------------------!
153      ! 1) Initialization                                                            !
154      !------------------------------------------------------------------------------!
155      !
156      DO ji = kideb , kiut
157         ! is there snow or not
158         isnow(ji)= INT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  )
159         ! surface temperature of fusion
160!!gm ???  ztfs(ji) = rtt !!!????
161         ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt
162         ! layer thickness
163         zh_i(ji) = ht_i_b(ji) / nlay_i
164         zh_s(ji) = ht_s_b(ji) / nlay_s
165      END DO
166
167      !--------------------
168      ! Ice / snow layers
169      !--------------------
170
171      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
172      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
173
174      DO layer = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
175         DO ji = kideb , kiut
176            z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s
177         END DO
178      END DO
179
180      DO layer = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
181         DO ji = kideb , kiut
182            z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i
183         END DO
184      END DO
185      !
186      !------------------------------------------------------------------------------|
187      ! 2) Radiations                                                                |
188      !------------------------------------------------------------------------------|
189      !
190      !-------------------
191      ! Computation of i0
192      !-------------------
193      ! i0 describes the fraction of solar radiation which does not contribute
194      ! to the surface energy budget but rather penetrates inside the ice.
195      ! We assume that no radiation is transmitted through the snow
196      ! If there is no no snow
197      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
198      ! zftrice = io.qsr_ice       is below the surface
199      ! fstbif  = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
200
201      DO ji = kideb , kiut
202         ! switches
203         isnow(ji) = INT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  ) 
204         ! hs > 0, isnow = 1
205         zhsu (ji) = hnzst  ! threshold for the computation of i0
206         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )     
207
208         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
209         !fr1_i0_1d = i0 for a thin ice surface
210         !fr1_i0_2d = i0 for a thick ice surface
211         !            a function of the cloud cover
212         !
213         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0)
214         !formula used in Cice
215      END DO
216
217      !-------------------------------------------------------
218      ! Solar radiation absorbed / transmitted at the surface
219      ! Derivative of the non solar flux
220      !-------------------------------------------------------
221      DO ji = kideb , kiut
222         zfsw   (ji) =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
223         zftrice(ji) =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
224         dzf    (ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
225      END DO
226
227      !---------------------------------------------------------
228      ! Transmission - absorption of solar radiation in the ice
229      !---------------------------------------------------------
230
231      DO ji = kideb, kiut           ! snow initialization
232         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
233      END DO
234
235      DO layer = 1, nlay_s          ! Radiation through snow
236         DO ji = kideb, kiut
237            !                             ! radiation transmitted below the layer-th snow layer
238            zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )
239            !                             ! radiation absorbed by the layer-th snow layer
240            zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
241         END DO
242      END DO
243
244      DO ji = kideb, kiut           ! ice initialization
245         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
246      END DO
247
248      DO layer = 1, nlay_i          ! Radiation through ice
249         DO ji = kideb, kiut
250            !                             ! radiation transmitted below the layer-th ice layer
251            zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )
252            !                             ! radiation absorbed by the layer-th ice layer
253            zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
254         END DO
255      END DO
256
257      DO ji = kideb, kiut           ! Radiation transmitted below the ice
258         fstbif_1d(ji) = fstbif_1d(ji) + zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji)
259      END DO
260
261      ! +++++
262      ! just to check energy conservation
263      DO ji = kideb, kiut
264         ii                = MOD( npb(ji) - 1, jpi ) + 1
265         ij                = ( npb(ji) - 1 ) / jpi + 1
266         fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i)
267      END DO
268      ! +++++
269
270      DO layer = 1, nlay_i
271         DO ji = kideb, kiut
272            radab(ji,layer) = zradab_i(ji,layer)
273         END DO
274      END DO
275
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)*isnow(ji) &
380               + zkappa_i(ji,0)*(1.0-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) &
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)-zeps)) 
662
663            ! surface temperature
664            isnow(ji)     = INT(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) .LT. ztfs(ji)) &
667               t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1)   &
668               &          + (1.0-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.0)
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         !                                ! update of latent heat fluxes
720         qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) )
721         !                                ! surface ice conduction flux
722         isnow(ji)       = INT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) )  )
723         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji))   &
724            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_b(ji,1) - t_su_b(ji))
725         !                                ! bottom ice conduction flux
726         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
727      END DO
728
729      !-------------------------!
730      ! Heat conservation       !
731      !-------------------------!
732      IF( con_i ) THEN
733         DO ji = kideb, kiut
734            ! Upper snow value
735            fc_s(ji,0) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 
736            ! Bott. snow value
737            fc_s(ji,1) = - isnow(ji)* zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 
738         END DO
739         DO ji = kideb, kiut         ! Upper ice layer
740            fc_i(ji,0) = - isnow(ji) * &  ! interface flux if there is snow
741               ( zkappa_i(ji,0)  * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) &
742               - ( 1.0 - isnow(ji) ) * ( zkappa_i(ji,0) * & 
743               zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not
744         END DO
745         DO layer = 1, nlay_i - 1         ! Internal ice layers
746            DO ji = kideb, kiut
747               fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) )
748               ii = MOD( npb(ji) - 1, jpi ) + 1
749               ij = ( npb(ji) - 1 ) / jpi + 1
750            END DO
751         END DO
752         DO ji = kideb, kiut         ! Bottom ice layers
753            fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
754         END DO
755      ENDIF
756      !
757   END SUBROUTINE lim_thd_dif
758
759#else
760   !!----------------------------------------------------------------------
761   !!                   Dummy Module                 No LIM-3 sea-ice model
762   !!----------------------------------------------------------------------
763CONTAINS
764   SUBROUTINE lim_thd_dif          ! Empty routine
765   END SUBROUTINE lim_thd_dif
766#endif
767   !!======================================================================
768END MODULE limthd_dif
Note: See TracBrowser for help on using the repository browser.