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 tags/nemo_v3_2/nemo_v3_2/NEMO/LIM_SRC_3 – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/LIM_SRC_3/limthd_dif.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm

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