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 @ 8554

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

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