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.
icethd_zdf.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90 @ 8531

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

changes in style - part6 - more clarity (still not finished)

File size: 39.8 KB
Line 
1MODULE icethd_zdf
2   !!======================================================================
3   !!                       ***  MODULE icethd_zdf ***
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            ! sea-ice: variables
21   USE ice1D          ! sea-ice: thermodynamics
22   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   ice_thd_zdf        ! called by ice_thd
31   PUBLIC   ice_thd_zdf_init   ! called by ice_stp
32
33   !!** namelist (namthd_zdf) **
34   LOGICAL  ::   ln_zdf_Beer      ! Heat diffusion follows a Beer Law
35   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964)
36   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007)
37   REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K]
38   REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
39   LOGICAL  ::   ln_dqns_i        ! change non-solar surface flux with changing surface temperature (T) or not (F)
40
41   !!----------------------------------------------------------------------
42   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
43   !! $Id: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE ice_thd_zdf
49      !!------------------------------------------------------------------
50      !!                ***  ROUTINE ice_thd_zdf  ***
51      !! ** Purpose :
52      !!           This routine determines the time evolution of snow and sea-ice
53      !!           temperature profiles.
54      !! ** Method  :
55      !!           This is done by solving the heat equation diffusion with
56      !!           a Neumann boundary condition at the surface and a Dirichlet one
57      !!           at the bottom. Solar radiation is partially absorbed into the ice.
58      !!           The specific heat and thermal conductivities depend on ice salinity
59      !!           and temperature to take into account brine pocket melting. The
60      !!           numerical
61      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
62      !!           in the ice and snow system.
63      !!
64      !!           The successive steps of this routine are
65      !!           1.  Thermal conductivity at the interfaces of the ice layers
66      !!           2.  Internal absorbed radiation
67      !!           3.  Scale factors due to non-uniform grid
68      !!           4.  Kappa factors
69      !!           Then iterative procedure begins
70      !!           5.  specific heat in the ice
71      !!           6.  eta factors
72      !!           7.  surface flux computation
73      !!           8.  tridiagonal system terms
74      !!           9.  solving the tridiagonal system with Gauss elimination
75      !!           Iterative procedure ends according to a criterion on evolution
76      !!           of temperature
77      !!
78      !! ** Inputs / Ouputs : (global commons)
79      !!           surface temperature : t_su_1d
80      !!           ice/snow temperatures   : t_i_1d, t_s_1d
81      !!           ice salinities          : s_i_1d
82      !!           number of layers in the ice/snow: nlay_i, nlay_s
83      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
84      !!------------------------------------------------------------------
85      INTEGER ::   ji, jk         ! spatial loop index
86      INTEGER ::   numeq          ! current reference number of equation
87      INTEGER ::   minnumeqmin, maxnumeqmax
88      INTEGER ::   iconv          ! number of iterations in iterative procedure
89      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
90     
91      INTEGER, DIMENSION(jpij) ::   numeqmin   ! reference number of top equation
92      INTEGER, DIMENSION(jpij) ::   numeqmax   ! reference number of bottom equation
93     
94      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
95      REAL(wp) ::   zg1       =  2._wp        !
96      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
97      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
98      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
99      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
100      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C
101      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature
102      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
103      REAL(wp) ::   z1_hsu
104      REAL(wp) ::   zdti_max                  ! current maximal error on temperature
105      REAL(wp) ::   zcpi                      ! Ice specific heat
106      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
107      REAL(wp) ::   zfac                      ! dummy factor
108     
109      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
110      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
111      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
112      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
113      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
114      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
115      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function
116      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
117      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
118     
119      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
120      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
121      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
122      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
123      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
124      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
125      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
126      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
127      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
128      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
129      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
130      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
131      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
132      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
133      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
134      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
135      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
136      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
137      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
138     
139      ! Mono-category
140      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
141      REAL(wp) :: zhe           ! dummy factor
142      REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity
143      !!------------------------------------------------------------------     
144
145      ! --- diag error on heat diffusion - PART 1 --- !
146      DO ji = 1, nidx
147         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
148            &           SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
149      END DO
150
151      !------------------------------------------------------------------------------!
152      ! 1) Initialization                                                            !
153      !------------------------------------------------------------------------------!
154      DO ji = 1, nidx
155         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
156         ! layer thickness
157         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
158         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
159      END DO
160      !
161      WHERE( zh_i(1:nidx) >= epsi10 )   ;   z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx)
162      ELSEWHERE                         ;   z1_h_i(1:nidx) = 0._wp
163      END WHERE
164
165      WHERE( zh_s(1:nidx) >= epsi10 )   ;   z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx)
166      ELSEWHERE                         ;   z1_h_s(1:nidx) = 0._wp
167      END WHERE
168      !
169      ! temperatures
170      ztsub  (1:nidx)   = t_su_1d(1:nidx)  ! temperature at the previous iteration
171      ztsold (1:nidx,:) = t_s_1d(1:nidx,:)   ! Old snow temperature
172      ztiold (1:nidx,:) = t_i_1d(1:nidx,:)   ! Old ice temperature
173      t_su_1d(1:nidx)   = MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary
174      !
175      !------------------------------------------------------------------------------|
176      ! 2) Radiation                                                              |
177      !------------------------------------------------------------------------------|
178      !
179      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
180      DO ji = 1, nidx
181         !-------------------
182         ! Computation of i0
183         !-------------------
184         ! i0 describes the fraction of solar radiation which does not contribute
185         ! to the surface energy budget but rather penetrates inside the ice.
186         ! We assume that no radiation is transmitted through the snow
187         ! If there is no no snow
188         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface
189         ! zftrice = io.qsr_ice       is below the surface
190         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
191         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
192         zfac = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )     
193         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )
194
195         !-------------------------------------------------------
196         ! Solar radiation absorbed / transmitted at the surface
197         ! Derivative of the non solar flux
198         !-------------------------------------------------------
199         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
200         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
201         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux
202         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
203      END DO
204
205      !---------------------------------------------------------
206      ! Transmission - absorption of solar radiation in the ice
207      !---------------------------------------------------------
208      zradtr_s(1:nidx,0) = zftrice(1:nidx)
209      DO jk = 1, nlay_s
210         DO ji = 1, nidx
211            !                             ! radiation transmitted below the layer-th snow layer
212            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
213            !                             ! radiation absorbed by the layer-th snow layer
214            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
215         END DO
216      END DO
217         
218      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )
219      DO jk = 1, nlay_i 
220         DO ji = 1, nidx
221            !                             ! radiation transmitted below the layer-th ice layer
222            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
223            !                             ! radiation absorbed by the layer-th ice layer
224            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
225         END DO
226      END DO
227
228      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice
229
230      !------------------------------------------------------------------------------|
231      !  3) Iterative procedure begins                                               |
232      !------------------------------------------------------------------------------|
233      !
234      iconv    =  0          ! number of iterations
235      zdti_max =  1000._wp   ! maximal value of error on all points
236      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )
237         !
238         iconv = iconv + 1
239         !
240         ztib(1:nidx,:) = t_i_1d(1:nidx,:)
241         ztsb(1:nidx,:) = t_s_1d(1:nidx,:)
242         !
243         !------------------------------------------------------------------------------|
244         ! 4) Sea ice thermal conductivity                                              |
245         !------------------------------------------------------------------------------|
246         !
247         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
248            !
249            DO ji = 1, nidx
250               ztcond_i(ji,0)      = rcdic + zbeta * s_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
251               ztcond_i(ji,nlay_i) = rcdic + zbeta * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
252            END DO
253            DO jk = 1, nlay_i-1
254               DO ji = 1, nidx
255                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
256                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
257               END DO
258            END DO
259            !
260         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
261            !
262            DO ji = 1, nidx
263               ztcond_i(ji,0)      = rcdic + 0.09_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
264                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
265               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )  &
266                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
267            END DO
268            DO jk = 1, nlay_i-1
269               DO ji = 1, nidx
270                  ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /          &
271                     &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   &
272                     &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
273               END DO
274            END DO
275            !
276         ENDIF
277         ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )       
278         !
279         !------------------------------------------------------------------------------|
280         !  5) G(he) - enhancement of thermal conductivity in mono-category case        |
281         !------------------------------------------------------------------------------|
282         !
283         ! Computation of effective thermal conductivity G(h)
284         ! Used in mono-category case only to simulate an ITD implicitly
285         ! Fichefet and Morales Maqueda, JGR 1997
286
287         zghe(1:nidx) = 1._wp
288         
289         SELECT CASE ( nn_monocat )
290
291         CASE ( 1 , 3 )
292
293            zepsilon = 0.1_wp
294            DO ji = 1, nidx
295   
296               ! Mean sea ice thermal conductivity
297               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )
298
299               ! Effective thickness he (zhe)
300               zhe  = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )
301
302               ! G(he)
303               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
304                  zghe(ji) =  MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )
305               ENDIF
306   
307            END DO
308
309         END SELECT
310         !
311         !------------------------------------------------------------------------------|
312         !  6) kappa factors                                                            |
313         !------------------------------------------------------------------------------|
314         !
315         !--- Snow
316         DO jk = 0, nlay_s-1
317            DO ji = 1, nidx
318               zkappa_s(ji,jk)  = zghe(ji) * rn_cnd_s * z1_h_s(ji)
319            END DO
320         END DO
321         DO ji = 1, nidx   ! Snow-ice interface
322            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
323            IF( zfac > epsi10 ) THEN
324               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
325            ELSE
326               zkappa_s(ji,nlay_s) = 0._wp
327            ENDIF
328         END DO
329
330         !--- Ice
331         DO jk = 0, nlay_i
332            DO ji = 1, nidx
333               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
334            END DO
335         END DO
336         DO ji = 1, nidx   ! Snow-ice interface
337            zkappa_i(ji,0)     = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
338         END DO
339         !
340         !------------------------------------------------------------------------------|
341         ! 7) Sea ice specific heat, eta factors                                        |
342         !------------------------------------------------------------------------------|
343         !
344         DO jk = 1, nlay_i
345            DO ji = 1, nidx
346               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
347               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
348            END DO
349         END DO
350
351         DO jk = 1, nlay_s
352            DO ji = 1, nidx
353               zeta_s(ji,jk)  = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
354            END DO
355         END DO
356         !
357         !------------------------------------------------------------------------------|
358         ! 8) surface flux computation                                                  |
359         !------------------------------------------------------------------------------|
360         !
361         IF ( ln_dqns_i ) THEN
362            DO ji = 1, nidx
363               ! update of the non solar flux according to the update in T_su
364               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
365            END DO
366         ENDIF
367
368         DO ji = 1, nidx
369            zf(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
370         END DO
371         !
372         !------------------------------------------------------------------------------|
373         ! 9) tridiagonal system terms                                                  |
374         !------------------------------------------------------------------------------|
375         !
376         !!layer denotes the number of the layer in the snow or in the ice
377         !!numeq denotes the reference number of the equation in the tridiagonal
378         !!system, terms of tridiagonal system are indexed as following :
379         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
380
381         !!ice interior terms (top equation has the same form as the others)
382
383         DO numeq=1,nlay_i+3
384            DO ji = 1, nidx
385               ztrid(ji,numeq,1) = 0.
386               ztrid(ji,numeq,2) = 0.
387               ztrid(ji,numeq,3) = 0.
388               zindterm(ji,numeq)= 0.
389               zindtbis(ji,numeq)= 0.
390               zdiagbis(ji,numeq)= 0.
391            ENDDO
392         ENDDO
393
394         DO numeq = nlay_s + 2, nlay_s + nlay_i 
395            DO ji = 1, nidx
396               jk                 = numeq - nlay_s - 1
397               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
398               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
399               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
400               zindterm(ji,numeq) =  ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
401            END DO
402         ENDDO
403
404         numeq =  nlay_s + nlay_i + 1
405         DO ji = 1, nidx
406            !!ice bottom term
407            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
408            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
409            ztrid(ji,numeq,3)  =  0.0
410            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
411               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
412         ENDDO
413
414
415         DO ji = 1, nidx
416            IF ( ht_s_1d(ji) > 0.0 ) THEN
417               !
418               !------------------------------------------------------------------------------|
419               !  snow-covered cells                                                          |
420               !------------------------------------------------------------------------------|
421               !
422               !!snow interior terms (bottom equation has the same form as the others)
423               DO numeq = 3, nlay_s + 1
424                  jk                  =  numeq - 1
425                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
426                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
427                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
428                  zindterm(ji,numeq)  =  ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
429               END DO
430
431               !!case of only one layer in the ice (ice equation is altered)
432               IF ( nlay_i == 1 ) THEN
433                  ztrid(ji,nlay_s+2,3)    =  0.0
434                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
435               ENDIF
436
437               IF ( t_su_1d(ji) < rt0 ) THEN
438
439                  !------------------------------------------------------------------------------|
440                  !  case 1 : no surface melting - snow present                                  |
441                  !------------------------------------------------------------------------------|
442                  numeqmin(ji)    =  1
443                  numeqmax(ji)    =  nlay_i + nlay_s + 1
444
445                  !!surface equation
446                  ztrid(ji,1,1)  = 0.0
447                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
448                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
449                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji)
450
451                  !!first layer of snow equation
452                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
453                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
454                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
455                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
456
457               ELSE 
458                  !
459                  !------------------------------------------------------------------------------|
460                  !  case 2 : surface is melting - snow present                                  |
461                  !------------------------------------------------------------------------------|
462                  !
463                  numeqmin(ji)    =  2
464                  numeqmax(ji)    =  nlay_i + nlay_s + 1
465
466                  !!first layer of snow equation
467                  ztrid(ji,2,1)  =  0.0
468                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
469                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
470                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
471                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
472               ENDIF
473            ELSE
474               !
475               !------------------------------------------------------------------------------|
476               !  cells without snow                                                          |
477               !------------------------------------------------------------------------------|
478               !
479               IF ( t_su_1d(ji) < rt0 ) THEN
480                  !
481                  !------------------------------------------------------------------------------|
482                  !  case 3 : no surface melting - no snow                                       |
483                  !------------------------------------------------------------------------------|
484                  !
485                  numeqmin(ji)      =  nlay_s + 1
486                  numeqmax(ji)      =  nlay_i + nlay_s + 1
487
488                  !!surface equation   
489                  ztrid(ji,numeqmin(ji),1)   =  0.0
490                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1   
491                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
492                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji)
493
494                  !!first layer of ice equation
495                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
496                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
497                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1) 
498                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 
499
500                  !!case of only one layer in the ice (surface & ice equations are altered)
501
502                  IF ( nlay_i == 1 ) THEN
503                     ztrid(ji,numeqmin(ji),1)    =  0.0
504                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
505                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
506                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
507                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
508                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
509
510                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1) *  &
511                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
512                  ENDIF
513
514               ELSE
515
516                  !
517                  !------------------------------------------------------------------------------|
518                  ! case 4 : surface is melting - no snow                                        |
519                  !------------------------------------------------------------------------------|
520                  !
521                  numeqmin(ji)    =  nlay_s + 2
522                  numeqmax(ji)    =  nlay_i + nlay_s + 1
523
524                  !!first layer of ice equation
525                  ztrid(ji,numeqmin(ji),1)      =  0.0
526                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
527                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
528                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1) *  &
529                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
530
531                  !!case of only one layer in the ice (surface & ice equations are altered)
532                  IF ( nlay_i == 1 ) THEN
533                     ztrid(ji,numeqmin(ji),1)  =  0.0
534                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
535                     ztrid(ji,numeqmin(ji),3)  =  0.0
536                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
537                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
538                  ENDIF
539
540               ENDIF
541            ENDIF
542
543         END DO
544         !
545         !------------------------------------------------------------------------------|
546         ! 10) tridiagonal system solving                                               |
547         !------------------------------------------------------------------------------|
548         !
549         ! Solve the tridiagonal system with Gauss elimination method.
550         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
551         ! McGraw-Hill 1984. 
552
553         maxnumeqmax = 0
554         minnumeqmin = nlay_i+5
555
556         DO ji = 1, nidx
557            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
558            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
559            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
560            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
561         END DO
562
563         DO jk = minnumeqmin+1, maxnumeqmax
564            DO ji = 1, nidx
565               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
566               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
567               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
568            END DO
569         END DO
570
571         DO ji = 1, nidx
572            ! ice temperatures
573            t_i_1d(ji,nlay_i) =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
574         END DO
575
576         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
577            DO ji = 1, nidx
578               jk    =  numeq - nlay_s - 1
579               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
580            END DO
581         END DO
582
583         DO ji = 1, nidx
584            ! snow temperatures     
585            IF( ht_s_1d(ji) > 0._wp ) THEN
586               t_s_1d(ji,nlay_s) =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
587                  &                 / zdiagbis(ji,nlay_s+1)
588            ENDIF
589            ! surface temperature
590            ztsub(ji) = t_su_1d(ji)
591            IF( t_su_1d(ji) < rt0 ) THEN
592               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
593                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
594            ENDIF
595         END DO
596         !
597         !--------------------------------------------------------------------------
598         !  11) Has the scheme converged ?, end of the iterative procedure         |
599         !--------------------------------------------------------------------------
600         !
601         ! check that nowhere it has started to melt
602         ! zdti_max is a measure of error, it has to be under zdti_bnd
603         zdti_max = 0._wp
604         DO ji = 1, nidx
605            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
606            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
607         END DO
608
609         DO jk  =  1, nlay_s
610            DO ji = 1, nidx
611               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
612               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
613            END DO
614         END DO
615
616         DO jk  =  1, nlay_i
617            DO ji = 1, nidx
618               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
619               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
620               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
621            END DO
622         END DO
623
624         ! Compute spatial maximum over all errors
625         ! note that this could be optimized substantially by iterating only the non-converging points
626         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
627
628      END DO  ! End of the do while iterative procedure
629
630      IF( ln_icectl .AND. lwp ) THEN
631         WRITE(numout,*) ' zdti_max : ', zdti_max
632         WRITE(numout,*) ' iconv    : ', iconv
633      ENDIF
634
635      !
636      !-------------------------------------------------------------------------!
637      !   12) Fluxes at the interfaces                                          !
638      !-------------------------------------------------------------------------!
639      DO ji = 1, nidx
640         !                                ! surface ice conduction flux
641         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
642            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
643         !                                ! bottom ice conduction flux
644         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
645      END DO
646
647      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
648      CALL ice_thd_enmelt
649
650      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
651      IF ( ln_dqns_i ) THEN
652         DO ji = 1, nidx
653            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
654         END DO
655      END IF
656
657      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
658      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
659      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
660      DO ji = 1, nidx
661         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
662            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
663
664         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
665            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
666         ELSE                          ! case T_su = 0degC
667            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
668         ENDIF
669         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
670
671         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
672         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
673
674      END DO   
675
676      ! --- Diagnostics SIMIP --- !
677      DO ji = 1, nidx
678         !--- Conduction fluxes (positive downwards)
679         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
680         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
681
682         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
683         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
684         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
685            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) + &
686               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
687         ELSE
688            t_si_1d(ji) = t_su_1d(ji)
689         ENDIF
690      END DO
691      !
692   END SUBROUTINE ice_thd_zdf
693
694
695   SUBROUTINE ice_thd_enmelt
696      !!-----------------------------------------------------------------------
697      !!                   ***  ROUTINE ice_thd_enmelt ***
698      !!                 
699      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
700      !!
701      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
702      !!-------------------------------------------------------------------
703      INTEGER  ::   ji, jk   ! dummy loop indices
704      REAL(wp) ::   ztmelts  ! local scalar
705      !!-------------------------------------------------------------------
706      !
707      DO jk = 1, nlay_i             ! Sea ice energy of melting
708         DO ji = 1, nidx
709            ztmelts      = - tmut  * s_i_1d(ji,jk)
710            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
711                                                                !   (sometimes dif scheme produces abnormally high temperatures)   
712            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
713               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
714               &                    - rcp  *   ztmelts )
715         END DO
716      END DO
717      DO jk = 1, nlay_s             ! Snow energy of melting
718         DO ji = 1, nidx
719            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
720         END DO
721      END DO
722      !
723   END SUBROUTINE ice_thd_enmelt
724
725
726   SUBROUTINE ice_thd_zdf_init
727      !!-----------------------------------------------------------------------
728      !!                   ***  ROUTINE ice_thd_zdf_init ***
729      !!                 
730      !! ** Purpose :   Physical constants and parameters associated with
731      !!                ice thermodynamics
732      !!
733      !! ** Method  :   Read the namthd_zdf namelist and check the parameters
734      !!                called at the first timestep (nit000)
735      !!
736      !! ** input   :   Namelist namthd_zdf
737      !!-------------------------------------------------------------------
738      INTEGER  ::   ios   ! Local integer output status for namelist read
739      !!
740      NAMELIST/namthd_zdf/ ln_zdf_Beer, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i
741      !!-------------------------------------------------------------------
742      !
743      REWIND( numnam_ice_ref )              ! Namelist namthd_zdf in reference namelist : Ice thermodynamics
744      READ  ( numnam_ice_ref, namthd_zdf, IOSTAT = ios, ERR = 901)
745901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in reference namelist', lwp )
746
747      REWIND( numnam_ice_cfg )              ! Namelist namthd_zdf in configuration namelist : Ice thermodynamics
748      READ  ( numnam_ice_cfg, namthd_zdf, IOSTAT = ios, ERR = 902 )
749902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in configuration namelist', lwp )
750      IF(lwm) WRITE ( numoni, namthd_zdf )
751      !
752      !
753      IF(lwp) THEN                          ! control print
754         WRITE(numout,*) 'ice_thd_zdf_init: Ice vertical heat diffusion'
755         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
756         WRITE(numout,*) '   Namelist namthd_zdf:'
757         WRITE(numout,*) '      Diffusion follows a Beer Law                            ln_zdf_Beer  = ', ln_zdf_Beer
758         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64
759         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07
760         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
761         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
762         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i
763     ENDIF
764      !
765      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
766         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conductivity (ln_cndi_U64 or ln_cndi_P07)' )
767      ENDIF
768      !
769   END SUBROUTINE ice_thd_zdf_init
770
771#else
772   !!----------------------------------------------------------------------
773   !!                   Dummy Module                 No ESIM sea-ice model
774   !!----------------------------------------------------------------------
775#endif
776
777   !!======================================================================
778END MODULE icethd_zdf
Note: See TracBrowser for help on using the repository browser.