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/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5075

Last change on this file since 5075 was 5075, checked in by timgraham, 9 years ago

Upgraded branch to current head of trunk (r5072) so it can be used with the trunk

  • Property svn:keywords set to Id
File size: 40.5 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   USE sbc_oce, ONLY : lk_cpl
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   lim_thd_dif   ! called by lim_thd
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 )
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_1d
77      !!           ice/snow temperatures   : t_i_1d, t_s_1d
78      !!           ice salinities          : s_i_1d
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_1d, ht_s_1d
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, kiut   ! Start/End point on which the  the computation is applied
94
95      !! * Local variables
96      INTEGER ::   ji          ! spatial loop index
97      INTEGER ::   ii, ij      ! temporary dummy loop index
98      INTEGER ::   numeq       ! current reference number of equation
99      INTEGER ::   jk       ! vertical dummy loop index
100      INTEGER ::   nconv       ! number of iterations in iterative procedure
101      INTEGER ::   minnumeqmin, maxnumeqmax
102      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
103      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation
104      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow
105      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
106      REAL(wp) ::   zg1       =  2._wp        !
107      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
108      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
109      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
110      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
111      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered as 0°C
112      REAL(wp) ::   ztmelt_i    ! ice melting temperature
113      REAL(wp) ::   zerritmax   ! current maximal error on temperature
114      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point
115      REAL(wp), POINTER, DIMENSION(:) ::   ztsub       ! old surface temperature (before the iterative procedure )
116      REAL(wp), POINTER, DIMENSION(:) ::   ztsubit     ! surface temperature at previous iteration
117      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
118      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
119      REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface
120      REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function
121      REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function
122      REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature
123      REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4)
124      REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice
125      REAL(wp), POINTER, DIMENSION(:) ::   zihic, zhsu
126      REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity
127      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice
128      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice
129      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice
130      REAL(wp), POINTER, DIMENSION(:,:) ::   ztib        ! Old temperature in the ice
131      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice
132      REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
133      REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat
134      REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice
135      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow
136      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow
137      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow
138      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s      ! Eta factor in the snow
139      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp     ! Temporary temperature in the snow to check the convergence
140      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsb        ! Temporary temperature in the snow
141      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s         ! Vertical cotes of the layers in the snow
142      REAL(wp), POINTER, DIMENSION(:,:) ::   zswiterm    ! Independent term
143      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitbis    ! temporary independent term
144      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis
145      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid     ! tridiagonal system terms
146      ! diag errors on heat
147      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err
148      !!------------------------------------------------------------------     
149      !
150      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow )
151      CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
152      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
153      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
154      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0)
155      CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis  )
156      CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid )
157
158      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
159
160      ! --- diag error on heat diffusion - PART 1 --- !
161      zdq(:) = 0._wp ; zq_ini(:) = 0._wp     
162      DO ji = kideb, kiut
163         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
164            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
165      END DO
166
167      !------------------------------------------------------------------------------!
168      ! 1) Initialization                                                            !
169      !------------------------------------------------------------------------------!
170      ! clem clean: replace just ztfs by rtt
171      DO ji = kideb , kiut
172         ! is there snow or not
173         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  )
174         ! surface temperature of fusion
175         ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt
176         ! layer thickness
177         zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i )
178         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
179      END DO
180
181      !--------------------
182      ! Ice / snow layers
183      !--------------------
184
185      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
186      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
187
188      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
189         DO ji = kideb , kiut
190            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s )
191         END DO
192      END DO
193
194      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
195         DO ji = kideb , kiut
196            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i )
197         END DO
198      END DO
199      !
200      !------------------------------------------------------------------------------|
201      ! 2) Radiations                                                                |
202      !------------------------------------------------------------------------------|
203      !
204      !-------------------
205      ! Computation of i0
206      !-------------------
207      ! i0 describes the fraction of solar radiation which does not contribute
208      ! to the surface energy budget but rather penetrates inside the ice.
209      ! We assume that no radiation is transmitted through the snow
210      ! If there is no no snow
211      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
212      ! zftrice = io.qsr_ice       is below the surface
213      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
214
215      DO ji = kideb , kiut
216         ! switches
217         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )  ) 
218         ! hs > 0, isnow = 1
219         zhsu (ji) = hnzst  ! threshold for the computation of i0
220         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) )     
221
222         i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
223         !fr1_i0_1d = i0 for a thin ice surface
224         !fr1_i0_2d = i0 for a thick ice surface
225         !            a function of the cloud cover
226         !
227         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0)
228         !formula used in Cice
229      END DO
230
231      !-------------------------------------------------------
232      ! Solar radiation absorbed / transmitted at the surface
233      ! Derivative of the non solar flux
234      !-------------------------------------------------------
235      DO ji = kideb , kiut
236         zfsw   (ji) =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
237         zftrice(ji) =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
238         dzf    (ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
239      END DO
240
241      !---------------------------------------------------------
242      ! Transmission - absorption of solar radiation in the ice
243      !---------------------------------------------------------
244
245      DO ji = kideb, kiut           ! snow initialization
246         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
247      END DO
248
249      DO jk = 1, nlay_s          ! Radiation through snow
250         DO ji = kideb, kiut
251            !                             ! radiation transmitted below the layer-th snow layer
252            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
253            !                             ! radiation absorbed by the layer-th snow layer
254            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
255         END DO
256      END DO
257
258      DO ji = kideb, kiut           ! ice initialization
259         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) )
260      END DO
261
262      DO jk = 1, nlay_i          ! Radiation through ice
263         DO ji = kideb, kiut
264            !                             ! radiation transmitted below the layer-th ice layer
265            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
266            !                             ! radiation absorbed by the layer-th ice layer
267            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
268         END DO
269      END DO
270
271      DO ji = kideb, kiut           ! Radiation transmitted below the ice
272         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
273      END DO
274
275      !
276      !------------------------------------------------------------------------------|
277      !  3) Iterative procedure begins                                               |
278      !------------------------------------------------------------------------------|
279      !
280      DO ji = kideb, kiut        ! Old surface temperature
281         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
282         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
283         t_su_1d   (ji) =  MIN( t_su_1d(ji), ztfs(ji) - ztsu_err )  ! necessary
284         zerrit   (ji) =  1000._wp                                ! initial value of error
285      END DO
286
287      DO jk = 1, nlay_s       ! Old snow temperature
288         DO ji = kideb , kiut
289            ztsb(ji,jk) =  t_s_1d(ji,jk)
290         END DO
291      END DO
292
293      DO jk = 1, nlay_i       ! Old ice temperature
294         DO ji = kideb , kiut
295            ztib(ji,jk) =  t_i_1d(ji,jk)
296         END DO
297      END DO
298
299      nconv     =  0           ! number of iterations
300      zerritmax =  1000._wp    ! maximal value of error on all points
301
302      DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd )
303         !
304         nconv = nconv + 1
305         !
306         !------------------------------------------------------------------------------|
307         ! 4) Sea ice thermal conductivity                                              |
308         !------------------------------------------------------------------------------|
309         !
310         IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula
311            DO ji = kideb , kiut
312               ztcond_i(ji,0)        = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt)
313               ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin)
314            END DO
315            DO jk = 1, nlay_i-1
316               DO ji = kideb , kiut
317                  ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
318                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)
319                  ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin)
320               END DO
321            END DO
322         ENDIF
323
324         IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
325            DO ji = kideb , kiut
326               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt )   &
327                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rtt ) 
328               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
329            END DO
330            DO jk = 1, nlay_i-1
331               DO ji = kideb , kiut
332                  ztcond_i(ji,jk) = rcdic +                                                                     & 
333                     &                 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                          &
334                     &                 / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)   &
335                     &               - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt ) 
336                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
337               END DO
338            END DO
339            DO ji = kideb , kiut
340               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt)   &
341                  &                        - 0.011_wp * ( t_bo_1d(ji) - rtt ) 
342               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
343            END DO
344         ENDIF
345         !
346         !------------------------------------------------------------------------------|
347         !  5) kappa factors                                                            |
348         !------------------------------------------------------------------------------|
349         !
350         DO ji = kideb, kiut
351
352            !-- Snow kappa factors
353            zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji))
354            zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji))
355         END DO
356
357         DO jk = 1, nlay_s-1
358            DO ji = kideb , kiut
359               zkappa_s(ji,jk)  = 2.0 * rcdsn / &
360                  MAX(epsi10,2.0*zh_s(ji))
361            END DO
362         END DO
363
364         DO jk = 1, nlay_i-1
365            DO ji = kideb , kiut
366               !-- Ice kappa factors
367               zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ &
368                  MAX(epsi10,2.0*zh_i(ji)) 
369            END DO
370         END DO
371
372         DO ji = kideb , kiut
373            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))
374            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))
375            !-- Interface
376            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, &
377               (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji)))
378            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) &
379               + zkappa_i(ji,0)*REAL( 1 - isnow(ji) )
380         END DO
381         !
382         !------------------------------------------------------------------------------|
383         ! 6) Sea ice specific heat, eta factors                                        |
384         !------------------------------------------------------------------------------|
385         !
386         DO jk = 1, nlay_i
387            DO ji = kideb , kiut
388               ztitemp(ji,jk)   = t_i_1d(ji,jk)
389               zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ &
390                  MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10)
391               zeta_i(ji,jk)    = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), &
392                  epsi10)
393            END DO
394         END DO
395
396         DO jk = 1, nlay_s
397            DO ji = kideb , kiut
398               ztstemp(ji,jk) = t_s_1d(ji,jk)
399               zeta_s(ji,jk)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)
400            END DO
401         END DO
402         !
403         !------------------------------------------------------------------------------|
404         ! 7) surface flux computation                                                  |
405         !------------------------------------------------------------------------------|
406         !
407         IF( .NOT. lk_cpl ) THEN   !--- forced atmosphere case
408            DO ji = kideb , kiut
409               ! update of the non solar flux according to the update in T_su
410               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
411            END DO
412         ENDIF
413
414         ! Update incoming flux
415         DO ji = kideb , kiut
416            ! update incoming flux
417            zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation
418               + qns_ice_1d(ji)                   ! non solar total flux
419            ! (LWup, LWdw, SH, LH)
420         END DO
421
422         !
423         !------------------------------------------------------------------------------|
424         ! 8) tridiagonal system terms                                                  |
425         !------------------------------------------------------------------------------|
426         !
427         !!layer denotes the number of the layer in the snow or in the ice
428         !!numeq denotes the reference number of the equation in the tridiagonal
429         !!system, terms of tridiagonal system are indexed as following :
430         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
431
432         !!ice interior terms (top equation has the same form as the others)
433
434         DO numeq=1,nlay_i+3
435            DO ji = kideb , kiut
436               ztrid(ji,numeq,1) = 0.
437               ztrid(ji,numeq,2) = 0.
438               ztrid(ji,numeq,3) = 0.
439               zswiterm(ji,numeq)= 0.
440               zswitbis(ji,numeq)= 0.
441               zdiagbis(ji,numeq)= 0.
442            ENDDO
443         ENDDO
444
445         DO numeq = nlay_s + 2, nlay_s + nlay_i 
446            DO ji = kideb , kiut
447               jk              = numeq - nlay_s - 1
448               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk-1)
449               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + &
450                  zkappa_i(ji,jk))
451               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk)
452               zswiterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* &
453                  zradab_i(ji,jk)
454            END DO
455         ENDDO
456
457         numeq =  nlay_s + nlay_i + 1
458         DO ji = kideb , kiut
459            !!ice bottom term
460            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
461            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 &
462               +  zkappa_i(ji,nlay_i-1) )
463            ztrid(ji,numeq,3)  =  0.0
464            zswiterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &
465               ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
466               *  t_bo_1d(ji) ) 
467         ENDDO
468
469
470         DO ji = kideb , kiut
471            IF ( ht_s_1d(ji).gt.0.0 ) THEN
472               !
473               !------------------------------------------------------------------------------|
474               !  snow-covered cells                                                          |
475               !------------------------------------------------------------------------------|
476               !
477               !!snow interior terms (bottom equation has the same form as the others)
478               DO numeq = 3, nlay_s + 1
479                  jk =  numeq - 1
480                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk-1)
481                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + &
482                     zkappa_s(ji,jk) )
483                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
484                  zswiterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* &
485                     zradab_s(ji,jk)
486               END DO
487
488               !!case of only one layer in the ice (ice equation is altered)
489               IF ( nlay_i.eq.1 ) THEN
490                  ztrid(ji,nlay_s+2,3)    =  0.0
491                  zswiterm(ji,nlay_s+2)   =  zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
492                     t_bo_1d(ji) 
493               ENDIF
494
495               IF ( t_su_1d(ji) .LT. rtt ) THEN
496
497                  !------------------------------------------------------------------------------|
498                  !  case 1 : no surface melting - snow present                                  |
499                  !------------------------------------------------------------------------------|
500                  zdifcase(ji)    =  1.0
501                  numeqmin(ji)    =  1
502                  numeqmax(ji)    =  nlay_i + nlay_s + 1
503
504                  !!surface equation
505                  ztrid(ji,1,1) = 0.0
506                  ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
507                  ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
508                  zswiterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji)
509
510                  !!first layer of snow equation
511                  ztrid(ji,2,1)  =  - zkappa_s(ji,0)*zg1s*zeta_s(ji,1)
512                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s)
513                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
514                  zswiterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
515
516               ELSE 
517                  !
518                  !------------------------------------------------------------------------------|
519                  !  case 2 : surface is melting - snow present                                  |
520                  !------------------------------------------------------------------------------|
521                  !
522                  zdifcase(ji)    =  2.0
523                  numeqmin(ji)    =  2
524                  numeqmax(ji)    =  nlay_i + nlay_s + 1
525
526                  !!first layer of snow equation
527                  ztrid(ji,2,1)  =  0.0
528                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + &
529                     zkappa_s(ji,0) * zg1s )
530                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
531                  zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            &
532                     ( zradab_s(ji,1) +                         &
533                     zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
534               ENDIF
535            ELSE
536               !
537               !------------------------------------------------------------------------------|
538               !  cells without snow                                                          |
539               !------------------------------------------------------------------------------|
540               !
541               IF (t_su_1d(ji) .LT. rtt) THEN
542                  !
543                  !------------------------------------------------------------------------------|
544                  !  case 3 : no surface melting - no snow                                       |
545                  !------------------------------------------------------------------------------|
546                  !
547                  zdifcase(ji)      =  3.0
548                  numeqmin(ji)      =  nlay_s + 1
549                  numeqmax(ji)      =  nlay_i + nlay_s + 1
550
551                  !!surface equation   
552                  ztrid(ji,numeqmin(ji),1)   =  0.0
553                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1   
554                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
555                  zswiterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
556
557                  !!first layer of ice equation
558                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
559                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) & 
560                     + zkappa_i(ji,0) * zg1 )
561                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1) 
562                  zswiterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 
563
564                  !!case of only one layer in the ice (surface & ice equations are altered)
565
566                  IF (nlay_i.eq.1) THEN
567                     ztrid(ji,numeqmin(ji),1)    =  0.0
568                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0)*2.0
569                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0)*2.0
570                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0)*2.0*zeta_i(ji,1)
571                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + &
572                        zkappa_i(ji,1))
573                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
574
575                     zswiterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* &
576                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) )
577                  ENDIF
578
579               ELSE
580
581                  !
582                  !------------------------------------------------------------------------------|
583                  ! case 4 : surface is melting - no snow                                        |
584                  !------------------------------------------------------------------------------|
585                  !
586                  zdifcase(ji)    =  4.0
587                  numeqmin(ji)    =  nlay_s + 2
588                  numeqmax(ji)    =  nlay_i + nlay_s + 1
589
590                  !!first layer of ice equation
591                  ztrid(ji,numeqmin(ji),1)      =  0.0
592                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* &
593                     zg1) 
594                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
595                  zswiterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
596                     zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
597
598                  !!case of only one layer in the ice (surface & ice equations are altered)
599                  IF (nlay_i.eq.1) THEN
600                     ztrid(ji,numeqmin(ji),1)  =  0.0
601                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + &
602                        zkappa_i(ji,1))
603                     ztrid(ji,numeqmin(ji),3)  =  0.0
604                     zswiterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* &
605                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) &
606                        + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0
607                  ENDIF
608
609               ENDIF
610            ENDIF
611
612         END DO
613
614         !
615         !------------------------------------------------------------------------------|
616         ! 9) tridiagonal system solving                                                |
617         !------------------------------------------------------------------------------|
618         !
619
620         ! Solve the tridiagonal system with Gauss elimination method.
621         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
622         ! McGraw-Hill 1984. 
623
624         maxnumeqmax = 0
625         minnumeqmin = nlay_i+5
626
627         DO ji = kideb , kiut
628            zswitbis(ji,numeqmin(ji)) =  zswiterm(ji,numeqmin(ji))
629            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
630            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
631            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
632         END DO
633
634         DO jk = minnumeqmin+1, maxnumeqmax
635            DO ji = kideb , kiut
636               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
637               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
638                  ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
639               zswitbis(ji,numeq)  =  zswiterm(ji,numeq) - ztrid(ji,numeq,1)* &
640                  zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
641            END DO
642         END DO
643
644         DO ji = kideb , kiut
645            ! ice temperatures
646            t_i_1d(ji,nlay_i)    =  zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
647         END DO
648
649         DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
650            DO ji = kideb , kiut
651               jk    =  numeq - nlay_s - 1
652               t_i_1d(ji,jk)  =  (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* &
653                  t_i_1d(ji,jk+1))/zdiagbis(ji,numeq)
654            END DO
655         END DO
656
657         DO ji = kideb , kiut
658            ! snow temperatures     
659            IF (ht_s_1d(ji).GT.0._wp) &
660               t_s_1d(ji,nlay_s)     =  (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
661               *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) &
662               *        MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 
663
664            ! surface temperature
665            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) )  )  )
666            ztsubit(ji) = t_su_1d(ji)
667            IF( t_su_1d(ji) < ztfs(ji) ) &
668               t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   &
669               &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
670         END DO
671         !
672         !--------------------------------------------------------------------------
673         !  10) Has the scheme converged ?, end of the iterative procedure         |
674         !--------------------------------------------------------------------------
675         !
676         ! check that nowhere it has started to melt
677         ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
678         DO ji = kideb , kiut
679            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp  )
680            zerrit(ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
681         END DO
682
683         DO jk  =  1, nlay_s
684            DO ji = kideb , kiut
685               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rtt ), 190._wp  )
686               zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk)))
687            END DO
688         END DO
689
690         DO jk  =  1, nlay_i
691            DO ji = kideb , kiut
692               ztmelt_i        = -tmut * s_i_1d(ji,jk) + rtt 
693               t_i_1d(ji,jk) =  MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp)
694               zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk)))
695            END DO
696         END DO
697
698         ! Compute spatial maximum over all errors
699         ! note that this could be optimized substantially by iterating only the non-converging points
700         zerritmax = 0._wp
701         DO ji = kideb, kiut
702            zerritmax = MAX( zerritmax, zerrit(ji) )   
703         END DO
704         IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
705
706      END DO  ! End of the do while iterative procedure
707
708      IF( ln_nicep .AND. lwp ) THEN
709         WRITE(numout,*) ' zerritmax : ', zerritmax
710         WRITE(numout,*) ' nconv     : ', nconv
711      ENDIF
712
713      !
714      !-------------------------------------------------------------------------!
715      !   11) Fluxes at the interfaces                                          !
716      !-------------------------------------------------------------------------!
717      DO ji = kideb, kiut
718         ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
719         IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) )
720         !                                ! surface ice conduction flux
721         isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )  )
722         fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
723            &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
724         !                                ! bottom ice conduction flux
725         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
726      END DO
727
728      !-----------------------------------------
729      ! Heat flux used to warm/cool ice in W.m-2
730      !-----------------------------------------
731      DO ji = kideb, kiut
732         IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC
733            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
734               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
735         ELSE                         ! case T_su = 0degC
736            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
737               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
738         ENDIF
739      END DO
740
741      ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
742      CALL lim_thd_enmelt( kideb, kiut )
743
744      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
745      DO ji = kideb, kiut
746         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
747            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) )
748         zhfx_err(ji)   = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 
749         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
750      END DO 
751
752      ! diagnose external surface (forced case) or bottom (forced case) from heat conservation
753      IF( .NOT. lk_cpl ) THEN   ! --- forced case: qns_ice and fc_su are diagnosed
754         !
755         DO ji = kideb, kiut
756            qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji)
757            fc_su     (ji) = fc_su(ji)      - zhfx_err(ji)
758         END DO
759         !
760      ELSE                      ! --- coupled case: ocean turbulent heat flux is diagnosed
761         !
762         DO ji = kideb, kiut
763            fhtur_1d  (ji) = fhtur_1d(ji)   - zhfx_err(ji)
764         END DO
765         !
766      ENDIF
767
768      ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m2)
769      DO ji = kideb, kiut
770         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
771         hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
772      END DO
773   
774      !
775      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
776      CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
777      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
778      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   &
779         &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
780      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
781      CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis )
782      CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid )
783      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
784
785   END SUBROUTINE lim_thd_dif
786
787   SUBROUTINE lim_thd_enmelt( kideb, kiut )
788      !!-----------------------------------------------------------------------
789      !!                   ***  ROUTINE lim_thd_enmelt ***
790      !!                 
791      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
792      !!
793      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
794      !!-------------------------------------------------------------------
795      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
796      !
797      INTEGER  ::   ji, jk   ! dummy loop indices
798      REAL(wp) ::   ztmelts  ! local scalar
799      !!-------------------------------------------------------------------
800      !
801      DO jk = 1, nlay_i             ! Sea ice energy of melting
802         DO ji = kideb, kiut
803            ztmelts      = - tmut  * s_i_1d(ji,jk) + rtt 
804            rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) )
805            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                             &
806               &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) )   &
807               &                   - rcp  *                 ( ztmelts-rtt )  ) 
808         END DO
809      END DO
810      DO jk = 1, nlay_s             ! Snow energy of melting
811         DO ji = kideb, kiut
812            q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus )
813         END DO
814      END DO
815      !
816   END SUBROUTINE lim_thd_enmelt
817
818#else
819   !!----------------------------------------------------------------------
820   !!                   Dummy Module                 No LIM-3 sea-ice model
821   !!----------------------------------------------------------------------
822CONTAINS
823   SUBROUTINE lim_thd_dif          ! Empty routine
824   END SUBROUTINE lim_thd_dif
825#endif
826   !!======================================================================
827END MODULE limthd_dif
Note: See TracBrowser for help on using the repository browser.