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.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMO/LIM_SRC_3/limthd.F90 @ 1156

Last change on this file since 1156 was 1156, checked in by rblod, 16 years ago

Update Id and licence information, see ticket #210
  • Property svn:keywords set to Id
File size: 46.2 KB
Line 
1MODULE limthd
2   !!======================================================================
3   !!                  ***  MODULE limthd   ***
4   !!              LIM thermo ice model : ice thermodynamic
5   !!======================================================================
6#if defined key_lim3
7   !!----------------------------------------------------------------------
8   !!   'key_lim3'                                      LIM3 sea-ice model
9   !!----------------------------------------------------------------------
10   !!   lim_thd      : thermodynamic of sea ice
11   !!   lim_thd_init : initialisation of sea-ice thermodynamic
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE phycst          ! physical constants
15   USE dom_oce         ! ocean space and time domain variables
16   USE lbclnk
17   USE in_out_manager  ! I/O manager
18   USE ice             ! LIM sea-ice variables
19   USE ice_oce         ! sea-ice/ocean variables
20   USE sbc_oce         ! Surface boundary condition: ocean fields
21   USE sbc_ice         ! Surface boundary condition: ice fields
22   USE thd_ice         ! LIM thermodynamic sea-ice variables
23   USE dom_ice         ! LIM sea-ice domain
24   USE iceini
25   USE limthd_dif
26   USE limthd_dh
27   USE limthd_sal
28   USE limthd_ent
29   USE limtab
30   USE par_ice
31   USE limvar
32   USE prtctl          ! Print control
33   USE lib_mpp
34
35   IMPLICIT NONE
36   PRIVATE
37
38   !! * Routine accessibility
39   PUBLIC lim_thd         ! called by lim_step
40
41   !! * Module variables
42   REAL(wp)  ::            &  ! constant values
43      epsi20 = 1e-20   ,  &
44      epsi16 = 1e-16   ,  &
45      epsi06 = 1e-06   ,  &
46      epsi04 = 1e-04   ,  &
47      zzero  = 0.e0     , &
48      zone   = 1.e0
49
50   !! * Substitutions
51#  include "domzgr_substitute.h90"
52#  include "vectopt_loop_substitute.h90"
53   !!----------------------------------------------------------------------
54   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2005)
55   !! $Id$
56   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58
59CONTAINS
60
61   SUBROUTINE lim_thd( kt )
62      !!-------------------------------------------------------------------
63      !!                ***  ROUTINE lim_thd  ***       
64      !! 
65      !! ** Purpose : This routine manages the ice thermodynamic.
66      !!         
67      !! ** Action : - Initialisation of some variables
68      !!             - Some preliminary computation (oceanic heat flux
69      !!               at the ice base, snow acc.,heat budget of the leads)
70      !!             - selection of the icy points and put them in an array
71      !!             - call lim_vert_ther for vert ice thermodynamic
72      !!             - back to the geographic grid
73      !!             - selection of points for lateral accretion
74      !!             - call lim_lat_acc  for the ice accretion
75      !!             - back to the geographic grid
76      !!     
77      !! ** References :
78      !!       H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
79      !!
80      !! History :
81      !!   1.0  !  00-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet)
82      !!   2.0  !  02-07 (C. Ethe, G. Madec) F90
83      !!   3.0  !  05-11 (M. Vancoppenolle ) Multi-layer thermodynamics,
84      !!                                     salinity variations
85      !!---------------------------------------------------------------------
86      INTEGER, INTENT(in) ::   kt     ! number of iteration
87      !! * Local variables
88      INTEGER  ::  ji, jj, jk, jl, nbpb   ! nb of icy pts for thermo. cal.
89
90      REAL(wp) ::  &
91         zfric_umin = 5e-03 ,  &   ! lower bound for the friction velocity
92         zfric_umax = 2e-02        ! upper bound for the friction velocity
93
94      REAL(wp) ::   &
95         zinda              ,  &   ! switch for test. the val. of concen.
96         zindb,                &   ! switches for test. the val of arg
97         zthsnice           ,  &
98         zfric_u            ,  &   ! friction velocity
99         zfnsol             ,  &   ! total non solar heat
100         zfontn             ,  &   ! heat flux from snow thickness
101         zfntlat, zpareff   ,  &   ! test. the val. of lead heat budget
102         zeps
103
104      REAL(wp) ::   &
105         zareamin
106
107      REAL(wp), DIMENSION(jpi,jpj) :: &
108         zhicifp            ,  &   ! ice thickness for outputs
109         zqlbsbq                   ! link with lead energy budget qldif
110
111      !!-------------------------------------------------------------------
112
113      IF( numit == nstart  )   CALL lim_thd_init  ! Initialization (first time-step only)
114
115      IF( kt == nit000 .AND. lwp ) THEN
116         WRITE(numout,*) 'limthd : Ice Thermodynamics'
117         WRITE(numout,*) '~~~~~~'
118      ENDIF
119
120      IF( numit == nstart  )   CALL lim_thd_sal_init  ! Initialization (first time-step only)
121      !------------------------------------------------------------------------------!
122      ! 1) Initialization of diagnostic variables                                    !
123      !------------------------------------------------------------------------------!
124      zeps = 1.0e-10
125      tatm_ice(:,:) = tatm_ice(:,:) + 273.15 ! convert C to K
126
127      !--------------------
128      ! 1.2) Heat content   
129      !--------------------
130      ! Change the units of heat content; from global units to
131      ! J.m3
132
133      DO jl = 1, jpl
134         DO jk = 1, nlay_i
135            DO jj = 1, jpj
136               DO ji = 1, jpi
137                  !Energy of melting q(S,T) [J.m-3]
138                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / area(ji,jj) / & 
139                     MAX( v_i(ji,jj,jl) , epsi06 ) * nlay_i
140                  !0 if no ice and 1 if yes
141                  zindb            = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) ) 
142                  !convert units ! very important that this line is here
143                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
144               END DO
145            END DO
146         END DO
147      END DO
148
149      DO jl = 1, jpl
150         DO jk = 1, nlay_s
151            DO jj = 1, jpj
152               DO ji = 1, jpi
153                  !Energy of melting q(S,T) [J.m-3]
154                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / area(ji,jj) / & 
155                     MAX( v_s(ji,jj,jl) , epsi06 ) * nlay_s
156                  !0 if no ice and 1 if yes
157                  zindb            = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) ) 
158                  !convert units ! very important that this line is here
159                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb 
160               END DO
161            END DO
162         END DO
163      END DO
164
165      !-----------------------------
166      ! 1.3) Set some dummies to 0
167      !-----------------------------
168      rdvosif(:,:) = 0.e0   ! variation of ice volume at surface
169      rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom
170      fdvolif(:,:) = 0.e0   ! total variation of ice volume
171      rdvonif(:,:) = 0.e0   ! lateral variation of ice volume
172      fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice
173      ffltbif(:,:) = 0.e0   ! linked with fstric
174      qfvbq  (:,:) = 0.e0   ! linked with fstric
175      rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area
176      rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area
177      hicifp (:,:) = 0.e0   ! daily thermodynamic ice production.
178      fsbri  (:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean
179      fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean
180      fseqv  (:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay
181
182      !-----------------------------------
183      ! 1.4) Compute global heat content
184      !-----------------------------------
185      qt_i_in(:,:)  = 0.e0
186      qt_s_in(:,:)  = 0.e0
187      qt_i_fin(:,:)  = 0.e0
188      qt_s_fin(:,:)  = 0.e0
189      sum_fluxq(:,:) = 0.e0
190      fatm(:,:) = 0.e0
191
192      ! 2) Partial computation of forcing for the thermodynamic sea ice model.      !
193      !-----------------------------------------------------------------------------!
194
195!CDIR NOVERRCHK
196      DO jj = 1, jpj
197!CDIR NOVERRCHK
198         DO ji = 1, jpi
199            zthsnice       = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) )
200            zindb          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) 
201            phicif(ji,jj)  = vt_i(ji,jj)
202            pfrld(ji,jj)   = 1.0 - at_i(ji,jj)
203            zinda          = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )
204
205            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
206            !           !  practically no "direct lateral ablation"
207            !           
208            !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean
209            !           !  temperature and turbulent mixing (McPhee, 1992)
210            ! friction velocity
211            zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 
212
213            ! here the drag will depend on ice thickness and type (0.006)
214            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 
215            ! also category dependent
216            !           !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead
217            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice
218            !                       
219
220            ! still need to be updated : fdtcn !!!!
221            !           !-- Lead heat budget (part 1, next one is in limthd_dh
222            !           !-- qldif -- (or qldif_1d in 1d routines)
223            zfontn         = sprecip(ji,jj) * lfus              ! energy of melting
224            zfnsol         = qns(ji,jj)                         ! total non solar flux
225            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj)                              &
226               &                               + zfnsol + fdtcn(ji,jj) - zfontn     &
227               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   &
228               &                               * ( 1.0 - at_i(ji,jj) ) * rdt_ice   
229
230            ! Positive heat budget is used for bottom ablation
231            zfntlat        = 1.0 - MAX( zzero , SIGN( zone ,  - qldif(ji,jj) ) )
232            != 1 if positive heat budget
233            zpareff        = 1.0 - zinda * zfntlat
234            != 0 if ice and positive heat budget and 1 if one of those two is
235            !false
236            zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / &
237               MAX( at_i(ji,jj) * rdt_ice , epsi16 )
238
239            ! Heat budget of the lead, energy transferred from ice to ocean
240            qldif  (ji,jj) = zpareff * qldif(ji,jj)
241            qdtcn  (ji,jj) = zpareff * qdtcn(ji,jj)
242
243            ! Energy needed to bring ocean surface layer until its freezing
244            ! qcmif, limflx
245            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda )
246
247            !  calculate oceanic heat flux (limthd_dh)
248            fbif   (ji,jj) = zindb * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )
249
250            ! computation of the daily thermodynamic ice production (only needed for output)
251            zhicifp(ji,jj) = ht_i(ji,jj,1) * at_i(ji,jj) 
252            zhicifp(ji,jj) = ht_i(ji,jj,1) * at_i(ji,jj)
253         END DO
254      END DO
255
256      !------------------------------------------------------------------------------!
257      ! 3) Select icy points and fulfill arrays for the vectorial grid.           
258      !------------------------------------------------------------------------------!
259
260      DO jl = 1, jpl      !loop over ice categories
261
262         IF( kt == nit000 .AND. lwp ) THEN
263            WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl 
264            WRITE(numout,*) ' ~~~~~~~~'
265         ENDIF
266
267         zareamin = 1.0e-10
268         nbpb = 0
269         DO jj = 1, jpj
270            DO ji = 1, jpi
271               IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN     
272                  nbpb      = nbpb  + 1
273                  npb(nbpb) = (jj - 1) * jpi + ji
274               ENDIF
275               ! debug point to follow
276               IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN
277                  jiindex_1d = nbpb
278               ENDIF
279            END DO
280         END DO
281
282         !------------------------------------------------------------------------------!
283         ! 4) Thermodynamic computation
284         !------------------------------------------------------------------------------!
285
286         IF( lk_mpp ) CALL mpp_ini_ice(nbpb)
287
288         IF (nbpb > 0) THEN  ! If there is no ice, do nothing.
289
290            !-------------------------
291            ! 4.1 Move to 1D arrays
292            !-------------------------
293
294            CALL tab_2d_1d( nbpb, at_i_b     (1:nbpb)     , at_i            , jpi, jpj, npb(1:nbpb) )
295            CALL tab_2d_1d( nbpb, a_i_b      (1:nbpb)     , a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) )
296            CALL tab_2d_1d( nbpb, ht_i_b     (1:nbpb)     , ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
297            CALL tab_2d_1d( nbpb, ht_s_b     (1:nbpb)     , ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
298
299            CALL tab_2d_1d( nbpb, t_su_b     (1:nbpb)     , t_su(:,:,jl), jpi, jpj, npb(1:nbpb) )
300            CALL tab_2d_1d( nbpb, sm_i_b     (1:nbpb)     , sm_i(:,:,jl), jpi, jpj, npb(1:nbpb) )
301            DO jk = 1, nlay_s
302               CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk)     , t_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
303               CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk)     , e_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
304            END DO
305            DO jk = 1, nlay_i
306               CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk)     , t_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
307               CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk)     , e_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
308               CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk)     , s_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
309            END DO
310
311            CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb)     , tatm_ice(:,:)   , jpi, jpj, npb(1:nbpb) )
312            CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb)     , qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
313            CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb)     , fr1_i0     , jpi, jpj, npb(1:nbpb) )
314            CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb)     , fr2_i0     , jpi, jpj, npb(1:nbpb) )
315            CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb)     , qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
316
317#if ! defined key_coupled
318            CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb)     , qla_ice(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
319            CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb)     , dqla_ice(:,:,jl)   , jpi, jpj, npb(1:nbpb) )
320#endif
321
322            CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb)     , dqns_ice(:,:,jl)   , jpi, jpj, npb(1:nbpb) )
323            CALL tab_2d_1d( nbpb, t_bo_b     (1:nbpb)     , t_bo       , jpi, jpj, npb(1:nbpb) )
324            CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb)     , sprecip    , jpi, jpj, npb(1:nbpb) ) 
325            CALL tab_2d_1d( nbpb, fbif_1d    (1:nbpb)     , fbif       , jpi, jpj, npb(1:nbpb) )
326            CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) )
327            CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) )
328            CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) )
329            CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) )
330
331            CALL tab_2d_1d( nbpb, fseqv_1d   (1:nbpb)     , fseqv      , jpi, jpj, npb(1:nbpb) )
332            CALL tab_2d_1d( nbpb, fsbri_1d   (1:nbpb)     , fsbri      , jpi, jpj, npb(1:nbpb) )
333            CALL tab_2d_1d( nbpb, fhbri_1d   (1:nbpb)     , fhbri      , jpi, jpj, npb(1:nbpb) )
334            CALL tab_2d_1d( nbpb, fstbif_1d  (1:nbpb)     , fstric     , jpi, jpj, npb(1:nbpb) )
335            CALL tab_2d_1d( nbpb, qfvbq_1d   (1:nbpb)     , qfvbq      , jpi, jpj, npb(1:nbpb) )
336
337            !--------------------------------
338            ! 4.3) Thermodynamic processes
339            !--------------------------------
340
341            IF ( con_i ) CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
342            IF ( con_i ) CALL lim_thd_glohec( qt_i_in , qt_s_in ,             &
343               q_i_layer_in , 1 , nbpb , jl )
344
345            !---------------------------------!
346            CALL lim_thd_dif(1,nbpb,jl)   ! Ice/Snow Temperature profile    !
347            !---------------------------------!
348
349            CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
350            ! compulsory for limthd_dh
351
352            IF ( con_i ) CALL lim_thd_glohec( qt_i_fin , qt_s_fin ,           &
353               q_i_layer_fin , 1 , nbpb , jl ) 
354            IF ( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl )
355
356            !---------------------------------!
357            CALL lim_thd_dh(1,nbpb,jl)    ! Ice/Snow thickness              !
358            !---------------------------------!
359
360            !---------------------------------!
361            CALL lim_thd_ent(1,nbpb,jl)   ! Ice/Snow enthalpy remapping     !
362            !---------------------------------!
363
364            !---------------------------------!
365            CALL lim_thd_sal(1,nbpb)      ! Ice salinity computation        !
366            !---------------------------------!
367
368            !           CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
369            IF ( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin,             &
370               q_i_layer_fin , 1 , nbpb , jl ) 
371            IF ( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl )
372
373            !--------------------------------
374            ! 4.4) Move 1D to 2D vectors
375            !--------------------------------
376
377            CALL tab_1d_2d( nbpb, at_i        , npb, at_i_b (1:nbpb), jpi, jpj )
378            CALL tab_1d_2d( nbpb, ht_i(:,:,jl), npb, ht_i_b(1:nbpb), jpi, jpj )
379            CALL tab_1d_2d( nbpb, ht_s(:,:,jl), npb, ht_s_b(1:nbpb), jpi, jpj )
380            CALL tab_1d_2d( nbpb, a_i (:,:,jl), npb, a_i_b(1:nbpb) , jpi, jpj )
381            CALL tab_1d_2d( nbpb, t_su(:,:,jl), npb, t_su_b(1:nbpb), jpi, jpj )
382            CALL tab_1d_2d( nbpb, sm_i(:,:,jl), npb, sm_i_b(1:nbpb), jpi, jpj )
383
384            DO jk = 1, nlay_s
385               CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b(1:nbpb,jk), jpi, jpj)
386               CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b(1:nbpb,jk), jpi, jpj)
387            END DO
388
389            DO jk = 1, nlay_i
390               CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b(1:nbpb,jk), jpi, jpj)
391               CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b(1:nbpb,jk), jpi, jpj)
392               CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b(1:nbpb,jk), jpi, jpj)
393            END DO
394
395            CALL tab_1d_2d( nbpb, fstric     , npb, fstbif_1d (1:nbpb)     , jpi, jpj )
396            CALL tab_1d_2d( nbpb, qldif      , npb, qldif_1d  (1:nbpb)     , jpi, jpj )
397            CALL tab_1d_2d( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj )
398            CALL tab_1d_2d( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj )
399            CALL tab_1d_2d( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj )
400            CALL tab_1d_2d( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj )
401            CALL tab_1d_2d( nbpb, rdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj )
402            CALL tab_1d_2d( nbpb, rdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj )
403            CALL tab_1d_2d( nbpb, fdvolif    , npb, dvlbq_1d  (1:nbpb)     , jpi, jpj )
404            CALL tab_1d_2d( nbpb, rdvonif    , npb, dvnbq_1d  (1:nbpb)     , jpi, jpj ) 
405            CALL tab_1d_2d( nbpb, fseqv      , npb, fseqv_1d  (1:nbpb)     , jpi, jpj )
406
407            IF ( num_sal .EQ. 2 ) THEN
408               CALL tab_1d_2d( nbpb, fsbri      , npb, fsbri_1d  (1:nbpb)     , jpi, jpj )
409               CALL tab_1d_2d( nbpb, fhbri      , npb, fhbri_1d  (1:nbpb)     , jpi, jpj )
410            ENDIF
411
412            !+++++
413            !temporary stuff for a dummyversion
414            CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb)      , jpi, jpj )
415            CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb)      , jpi, jpj )
416            CALL tab_1d_2d( nbpb, fsup2D     , npb, fsup     (1:nbpb)      , jpi, jpj )
417            CALL tab_1d_2d( nbpb, focea2D    , npb, focea    (1:nbpb)      , jpi, jpj )
418            CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new  (1:nbpb)      , jpi, jpj )
419            CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0    (1:nbpb)      , jpi, jpj )
420            CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj)
421            !+++++
422
423            IF( lk_mpp ) CALL mpp_comm_free(ncomm_ice) !RB necessary ??
424         ENDIF ! nbpb
425
426      END DO ! jl
427
428      !------------------------------------------------------------------------------!
429      ! 5) Global variables, diagnostics
430      !------------------------------------------------------------------------------!
431
432      !------------------------
433      ! 5.1) Ice heat content             
434      !------------------------
435
436      ! Enthalpies are global variables we have to readjust the units
437      DO jl = 1, jpl
438         DO jk = 1, nlay_i
439            DO jj = 1, jpj
440               DO ji = 1, jpi
441                  ! Change dimensions
442                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac
443
444                  ! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules
445                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &
446                     area(ji,jj) * a_i(ji,jj,jl) * &
447                     ht_i(ji,jj,jl) / nlay_i
448               END DO !ji
449            END DO !jj
450         END DO !jk
451      END DO !jl
452
453      !------------------------
454      ! 5.2) Snow heat content             
455      !------------------------
456
457      ! Enthalpies are global variables we have to readjust the units
458      DO jl = 1, jpl
459         DO jk = 1, nlay_s
460            DO jj = 1, jpj
461               DO ji = 1, jpi
462                  ! Change dimensions
463                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
464                  ! Multiply by volume, so that heat content in 10^9 Joules
465                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
466                     a_i(ji,jj,jl) * ht_s(ji,jj,jl)  / nlay_s
467               END DO !ji
468            END DO !jj
469         END DO !jk
470      END DO !jl
471
472      !----------------------------------
473      ! 5.3) Change thickness to volume
474      !----------------------------------
475      CALL lim_var_eqv2glo
476
477      !--------------------------------------------
478      ! 5.4) Diagnostic thermodynamic growth rates
479      !--------------------------------------------
480      d_v_i_thd (:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)    ! ice volumes
481      dv_dt_thd(:,:,:)    = d_v_i_thd(:,:,:) / rdt_ice * 86400.0
482
483      IF ( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
484
485      IF(ln_ctl) THEN   ! Control print
486         CALL prt_ctl_info(' ')
487         CALL prt_ctl_info(' - Cell values : ')
488         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
489         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd  : cell area :')
490         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd  : at_i      :')
491         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd  : vt_i      :')
492         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd  : vt_s      :')
493         DO jl = 1, jpl
494            CALL prt_ctl_info(' ')
495            CALL prt_ctl_info(' - Category : ', ivar1=jl)
496            CALL prt_ctl_info('   ~~~~~~~~~~')
497            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_thd  : a_i      : ')
498            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_thd  : ht_i     : ')
499            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_thd  : ht_s     : ')
500            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_thd  : v_i      : ')
501            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_thd  : v_s      : ')
502            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_thd  : e_s      : ')
503            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_thd  : t_su     : ')
504            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_thd  : t_snow   : ')
505            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_thd  : sm_i     : ')
506            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_thd  : smv_i    : ')
507            DO jk = 1, nlay_i
508               CALL prt_ctl_info(' ')
509               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
510               CALL prt_ctl_info('   ~~~~~~~')
511               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd  : t_i      : ')
512               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd  : e_i      : ')
513            END DO
514         END DO
515
516      ENDIF
517
518   END SUBROUTINE lim_thd
519
520   !===============================================================================
521
522   SUBROUTINE lim_thd_glohec(eti,ets,etilayer,kideb,kiut,jl)
523      !!-----------------------------------------------------------------------
524      !!                   ***  ROUTINE lim_thd_glohec ***
525      !!                 
526      !! ** Purpose :  Compute total heat content for each category
527      !!               Works with 1d vectors only
528      !!
529      !! history :
530      !!  9.9  ! 07-04 (M.Vancoppenolle) original code
531      !!-----------------------------------------------------------------------
532      !! * Local variables
533      INTEGER, INTENT(in) :: &
534         kideb, kiut,        &  ! bounds for the spatial loop
535         jl                     ! category number
536
537      REAL(wp), DIMENSION (jpij,jpl), INTENT(out) ::  &
538         eti, ets            ! vertically-summed heat content for ice /snow
539
540      REAL(wp), DIMENSION (jpij,jkmax), INTENT(out) ::  &
541         etilayer            ! heat content for ice layers
542
543      REAL(wp) :: &
544         zdes,    &          ! snow heat content increment (dummy)
545         zeps                ! very small value (1.e-10)
546
547      INTEGER  :: &
548         ji,jk               ! loop indices
549
550      !!-----------------------------------------------------------------------
551      eti(:,:) = 0.0
552      ets(:,:) = 0.0
553      zeps     = 1.0e-10
554
555      ! total q over all layers, ice [J.m-2]
556      DO jk = 1, nlay_i
557         DO ji = kideb, kiut
558            etilayer(ji,jk) = q_i_b(ji,jk) &
559               * ht_i_b(ji) / nlay_i
560            eti(ji,jl) = eti(ji,jl) + etilayer(ji,jk) 
561         END DO
562      END DO
563
564      ! total q over all layers, snow [J.m-2]
565      DO ji = kideb, kiut
566         zdes = q_s_b(ji,1) * ht_s_b(ji) / nlay_s 
567         ets(ji,jl) = ets(ji,jl) + zdes       
568      END DO
569
570      WRITE(numout,*) ' lim_thd_glohec '
571      WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice
572      WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice
573      WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) +         &
574         ets(jiindex_1d,jl) ) / rdt_ice
575
576   END SUBROUTINE lim_thd_glohec
577
578   !===============================================================================
579
580   SUBROUTINE lim_thd_con_dif(kideb,kiut,jl)
581      !!-----------------------------------------------------------------------
582      !!                   ***  ROUTINE lim_thd_con_dif ***
583      !!                 
584      !! ** Purpose :   Test energy conservation after heat diffusion
585      !!
586      !! history :
587      !!  9.9  ! 07-04 (M.Vancoppenolle) original code
588      !!-------------------------------------------------------------------
589      !! * Local variables
590      INTEGER, INTENT(in) ::        &
591         kideb, kiut,               &  !: bounds for the spatial loop
592         jl                            !: category number
593
594      REAL(wp)                 ::   &  !: ! goes to trash
595         meance,                    &  !: mean conservation error
596         max_cons_err,              &  !: maximum tolerated conservation error
597         max_surf_err                  !: maximum tolerated surface error
598
599      INTEGER ::                    &
600         numce                         !: number of points for which conservation
601      !  is violated
602      INTEGER  :: &
603         ji,jk,                     &  !: loop indices
604         zji, zjj
605      !!---------------------------------------------------------------------
606
607      max_cons_err =  1.0
608      max_surf_err =  0.001
609
610      !--------------------------
611      ! Increment of energy
612      !--------------------------
613      ! global
614      DO ji = kideb, kiut
615         dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl)  &
616            + qt_s_fin(ji,jl) - qt_s_in(ji,jl)
617      END DO
618      ! layer by layer
619      dq_i_layer(:,:)    = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
620
621      !----------------------------------------
622      ! Atmospheric heat flux, ice heat budget
623      !----------------------------------------
624
625      DO ji = kideb, kiut
626         zji                 = MOD( npb(ji) - 1, jpi ) + 1
627         zjj                 = ( npb(ji) - 1 ) / jpi + 1
628
629         fatm(ji,jl) = &
630            qnsr_ice_1d(ji)                  + & ! atm non solar
631            (1.0-i0(ji))*qsr_ice_1d(ji)          ! atm solar
632
633         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) &
634            - fstroc(zji,zjj,jl)
635      END DO
636
637      !--------------------
638      ! Conservation error
639      !--------------------
640
641      DO ji = kideb, kiut
642         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
643      END DO
644
645      numce = 0
646      meance = 0.0
647      DO ji = kideb, kiut
648         IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
649            numce = numce + 1
650            meance = meance + cons_error(ji,jl)
651         ENDIF
652      ENDDO
653      IF (numce .GT. 0 ) meance = meance / numce
654
655      WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
656      WRITE(numout,*) ' After lim_thd_dif, category : ', jl
657      WRITE(numout,*) ' Mean conservation error on big error points ', meance, &
658         numit
659      WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit
660
661      !-------------------------------------------------------
662      ! Surface error due to imbalance between Fatm and Fcsu
663      !-------------------------------------------------------
664      numce  = 0.0
665      meance = 0.0
666
667      DO ji = kideb, kiut
668         surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) )
669         IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. &
670            max_surf_err ) ) THEN
671            numce = numce + 1 
672            meance = meance + surf_error(ji,jl)
673         ENDIF
674      ENDDO
675      IF (numce .GT. 0 ) meance = meance / numce
676
677      WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err
678      WRITE(numout,*) ' After lim_thd_dif, category : ', jl
679      WRITE(numout,*) ' Mean surface error on big error points ', meance, numit
680      WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit
681
682      IF (jiindex_1D.GT.0) WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d)
683      IF (jiindex_1D.GT.0) WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl)
684      IF (jiindex_1D.GT.0) WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d)
685
686      !---------------------------------------
687      ! Write ice state in case of big errors
688      !---------------------------------------
689
690      DO ji = kideb, kiut
691         IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. &
692            ( cons_error(ji,jl) .GT. max_cons_err  ) ) THEN
693            zji                 = MOD( npb(ji) - 1, jpi ) + 1
694            zjj                 = ( npb(ji) - 1 ) / jpi + 1
695
696            WRITE(numout,*) ' alerte 1     '
697            WRITE(numout,*) ' Untolerated conservation / surface error after '
698            WRITE(numout,*) ' heat diffusion in the ice '
699            WRITE(numout,*) ' Category   : ', jl
700            WRITE(numout,*) ' zji , zjj  : ', zji, zjj
701            WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj)
702            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
703            WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)
704            WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) / rdt_ice
705            WRITE(numout,*) ' Fdt        : ', sum_fluxq(ji,jl)
706            WRITE(numout,*)
707            !        WRITE(numout,*) ' qt_i_in   : ', qt_i_in(ji,jl)
708            !        WRITE(numout,*) ' qt_s_in   : ', qt_s_in(ji,jl)
709            !        WRITE(numout,*) ' qt_i_fin  : ', qt_i_fin(ji,jl)
710            !        WRITE(numout,*) ' qt_s_fin  : ', qt_s_fin(ji,jl)
711            !        WRITE(numout,*) ' qt        : ', qt_i_fin(ji,jl) + &
712            !                                         qt_s_fin(ji,jl)
713            WRITE(numout,*) ' ht_i       : ', ht_i_b(ji)
714            WRITE(numout,*) ' ht_s       : ', ht_s_b(ji)
715            WRITE(numout,*) ' t_su       : ', t_su_b(ji)
716            WRITE(numout,*) ' t_s        : ', t_s_b(ji,1)
717            WRITE(numout,*) ' t_i        : ', t_i_b(ji,1:nlay_i)
718            WRITE(numout,*) ' t_bo       : ', t_bo_b(ji)
719            WRITE(numout,*) ' q_i        : ', q_i_b(ji,1:nlay_i)
720            WRITE(numout,*) ' s_i        : ', s_i_b(ji,1:nlay_i)
721            WRITE(numout,*) ' tmelts     : ', rtt - tmut*s_i_b(ji,1:nlay_i)
722            WRITE(numout,*)
723            WRITE(numout,*) ' Fluxes '
724            WRITE(numout,*) ' ~~~~~~ '
725            WRITE(numout,*) ' fatm       : ', fatm(ji,jl)
726            WRITE(numout,*) ' fc_su      : ', fc_su    (ji)
727            WRITE(numout,*) ' fstr_inice : ', qsr_ice_1d(ji)*i0(ji)
728            WRITE(numout,*) ' fc_bo      : ', - fc_bo_i  (ji)
729            WRITE(numout,*) ' foc        : ', fbif_1d(ji)
730            WRITE(numout,*) ' fstroc     : ', fstroc   (zji,zjj,jl)
731            WRITE(numout,*) ' i0         : ', i0(ji)
732            WRITE(numout,*) ' qsr_ice    : ', (1.0-i0(ji))*qsr_ice_1d(ji)
733            WRITE(numout,*) ' qns_ice    : ', qnsr_ice_1d(ji)
734            WRITE(numout,*) ' Conduction fluxes : '
735            WRITE(numout,*) ' fc_s      : ', fc_s(ji,0:nlay_s)
736            WRITE(numout,*) ' fc_i      : ', fc_i(ji,0:nlay_i)
737            WRITE(numout,*)
738            WRITE(numout,*) ' Layer by layer ... '
739            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - &
740               qt_s_in(ji,jl) )  & 
741               / rdt_ice
742            WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) -      &
743               fc_s(ji,0)
744            DO jk = 1, nlay_i
745               WRITE(numout,*) ' layer  : ', jk
746               WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice 
747               WRITE(numout,*) ' radab  : ', radab(ji,jk)
748               WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) -               &
749                  fc_i(ji,jk-1)
750               WRITE(numout,*) ' tot f  : ', fc_i(ji,jk) -               &
751                  fc_i(ji,jk-1) - radab(ji,jk)
752            END DO
753
754         ENDIF
755
756      END DO
757
758   END SUBROUTINE lim_thd_con_dif
759
760   !==============================================================================
761
762   SUBROUTINE lim_thd_con_dh(kideb,kiut,jl)
763      !!-----------------------------------------------------------------------
764      !!                   ***  ROUTINE lim_thd_con_dh  ***
765      !!                 
766      !! ** Purpose :   Test energy conservation after enthalpy redistr.
767      !!
768      !! history :
769      !!  9.9  ! 07-04 (M.Vancoppenolle) original code
770      !!-----------------------------------------------------------------------
771      !! * Local variables
772      INTEGER, INTENT(in) ::        &
773         kideb, kiut,               &  !: bounds for the spatial loop
774         jl                            !: category number
775
776      REAL(wp)                 ::   &  !: ! goes to trash
777         meance,                    &  !: mean conservation error
778         max_cons_err                  !: maximum tolerated conservation error
779
780      INTEGER ::                    &
781         numce                         !: number of points for which conservation
782      !  is violated
783      INTEGER  ::  ji, zji, zjj        ! loop indices
784      !!---------------------------------------------------------------------
785
786      max_cons_err = 1.0
787
788      !--------------------------
789      ! Increment of energy
790      !--------------------------
791      ! global
792      DO ji = kideb, kiut
793         dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl)  &
794            + qt_s_fin(ji,jl) - qt_s_in(ji,jl)
795      END DO
796      ! layer by layer
797      dq_i_layer(:,:)    = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
798
799      !----------------------------------------
800      ! Atmospheric heat flux, ice heat budget
801      !----------------------------------------
802
803      DO ji = kideb, kiut
804         zji                 = MOD( npb(ji) - 1, jpi ) + 1
805         zjj                 = ( npb(ji) - 1 ) / jpi + 1
806
807         fatm(ji,jl) = &
808            qnsr_ice_1d(ji)                  + & ! atm non solar
809            !         (1.0-i0(ji))*qsr_ice_1d(ji)          ! atm solar
810            qsr_ice_1d(ji)                       ! atm solar
811
812         sum_fluxq(ji,jl)     = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) & 
813            - fstroc(zji,zjj,jl) 
814         cons_error(ji,jl)   = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
815      END DO
816
817      !--------------------
818      ! Conservation error
819      !--------------------
820
821      DO ji = kideb, kiut
822         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
823      END DO
824
825      numce = 0
826      meance = 0.0
827      DO ji = kideb, kiut
828         IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
829            numce = numce + 1
830            meance = meance + cons_error(ji,jl)
831         ENDIF
832      ENDDO
833      IF (numce .GT. 0 ) meance = meance / numce
834
835      WRITE(numout,*) ' Error report - Category : ', jl
836      WRITE(numout,*) ' ~~~~~~~~~~~~ '
837      WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
838      WRITE(numout,*) ' After lim_thd_ent, category : ', jl
839      WRITE(numout,*) ' Mean conservation error on big error points ', meance, &
840         numit
841      WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit
842
843      !---------------------------------------
844      ! Write ice state in case of big errors
845      !---------------------------------------
846
847      DO ji = kideb, kiut
848         IF ( cons_error(ji,jl) .GT. max_cons_err  ) THEN
849            zji                 = MOD( npb(ji) - 1, jpi ) + 1
850            zjj                 = ( npb(ji) - 1 ) / jpi + 1
851
852            WRITE(numout,*) ' alerte 1 - category : ', jl
853            WRITE(numout,*) ' Untolerated conservation error after limthd_ent '
854            WRITE(numout,*) ' zji , zjj  : ', zji, zjj
855            WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj)
856            WRITE(numout,*) ' * '
857            WRITE(numout,*) ' Ftotal     : ', sum_fluxq(ji,jl)
858            WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) / rdt_ice
859            WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice
860            WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice
861            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
862            WRITE(numout,*) ' * '
863            WRITE(numout,*) ' Fluxes        --- : '
864            WRITE(numout,*) ' fatm       : ', fatm(ji,jl)
865            WRITE(numout,*) ' foce       : ', fbif_1d(ji)
866            WRITE(numout,*) ' fres       : ', ftotal_fin(ji)
867            WRITE(numout,*) ' fhbri      : ', fhbricat(zji,zjj,jl)
868            WRITE(numout,*) ' * '
869            WRITE(numout,*) ' Heat contents --- : '
870            WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) / rdt_ice
871            WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) / rdt_ice
872            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + &
873               qt_s_in(ji,jl) ) / rdt_ice
874            WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) / rdt_ice
875            WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) / rdt_ice
876            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + &
877               qt_s_fin(ji,jl) ) / rdt_ice
878            WRITE(numout,*) ' * '
879            WRITE(numout,*) ' Ice variables --- : '
880            WRITE(numout,*) ' ht_i       : ', ht_i_b(ji)
881            WRITE(numout,*) ' ht_s       : ', ht_s_b(ji)
882            WRITE(numout,*) ' dh_s_tot  : ', dh_s_tot(ji)
883            WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji)
884            WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
885            WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
886
887         ENDIF
888
889      END DO
890
891   END SUBROUTINE lim_thd_con_dh
892   !==============================================================================
893
894   SUBROUTINE lim_thd_enmelt(kideb,kiut)
895      !!-----------------------------------------------------------------------
896      !!                   ***  ROUTINE lim_thd_enmelt ***
897      !!                 
898      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3)
899      !!
900      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
901      !!
902      !! history : Martin Vancoppenolle, May 2007
903      !!-------------------------------------------------------------------
904      INTEGER, INTENT(in) ::        &
905         kideb, kiut                   !: bounds for the spatial loop
906
907      REAL(wp)                 ::   &  !: goes to trash
908         ztmelts               ,    &  !: sea ice freezing point in K
909         zeps 
910
911      INTEGER             ::        &
912         ji,                        &  !: spatial loop index
913         jk                            !: vertical index
914
915      !!-------------------------------------------------------------------
916      zeps = 1.0e-10
917
918      ! Sea ice energy of melting
919      DO jk = 1, nlay_i
920         DO ji = kideb, kiut
921            ztmelts      =   - tmut * s_i_b(ji,jk) + rtt 
922            q_i_b(ji,jk) =   rhoic*( cpic    * ( ztmelts - t_i_b(ji,jk) )  &
923               + lfus  * ( 1.0 - (ztmelts-rtt)/MIN((t_i_b(ji,jk)-rtt),-zeps) )  &
924               - rcp      * ( ztmelts-rtt  ) ) 
925         END DO !ji
926      END DO !jk
927
928      ! Snow energy of melting
929      DO jk = 1, nlay_s
930         DO ji = kideb,kiut
931            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
932         END DO !ji
933      END DO !jk
934
935   END SUBROUTINE lim_thd_enmelt
936
937   !==============================================================================
938
939   SUBROUTINE lim_thd_init
940
941      !!-----------------------------------------------------------------------
942      !!                   ***  ROUTINE lim_thd_init ***
943      !!                 
944      !! ** Purpose :   Physical constants and parameters linked to the ice
945      !!      thermodynamics
946      !!
947      !! ** Method  :   Read the namicethd namelist and check the ice-thermo
948      !!       parameter values called at the first timestep (nit000)
949      !!
950      !! ** input   :   Namelist namicether
951      !!
952      !! history :
953      !!  8.5  ! 03-08 (C. Ethe) original code
954      !!-------------------------------------------------------------------
955      NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, &
956         &                hicmin, hiclim, amax  ,    &
957         &                sbeta  , parlat, hakspl, hibspl, exld,  &
958         &                hakdif, hnzst  , thth  , parsub, alphs, betas, & 
959         &                kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi
960      !!-------------------------------------------------------------------
961
962      ! Define the initial parameters
963      ! -------------------------
964      REWIND( numnam_ice )
965      READ  ( numnam_ice , namicethd )
966      IF (lwp) THEN
967         WRITE(numout,*)
968         WRITE(numout,*)'lim_thd_init : ice parameters for ice thermodynamic computation '
969         WRITE(numout,*)'~~~~~~~~~~~~'
970         WRITE(numout,*)'   maximum melting at the bottom                           hmelt        = ', hmelt
971         WRITE(numout,*)'   ice thick. for lateral accretion in NH (SH)             hiccrit(1/2) = ', hiccrit
972         WRITE(numout,*)'   Frazil ice thickness as a function of wind or not       fraz_swi     = ', fraz_swi
973         WRITE(numout,*)'   Maximum proportion of frazil ice collecting at bottom   maxfrazb     = ', maxfrazb
974         WRITE(numout,*)'   Thresold relative drift speed for collection of frazil  vfrazb       = ', vfrazb
975         WRITE(numout,*)'   Squeezing coefficient for collection of frazil          Cfrazb       = ', Cfrazb
976         WRITE(numout,*)'   ice thick. corr. to max. energy stored in brine pocket  hicmin       = ', hicmin 
977         WRITE(numout,*)'   minimum ice thickness                                   hiclim       = ', hiclim 
978         WRITE(numout,*)'   maximum lead fraction                                   amax         = ', amax
979         WRITE(numout,*)'   numerical carac. of the scheme for diffusion in ice '
980         WRITE(numout,*)'   Cranck-Nicholson (=0.5), implicit (=1), explicit (=0)   sbeta        = ', sbeta
981         WRITE(numout,*)'   percentage of energy used for lateral ablation          parlat       = ', parlat
982         WRITE(numout,*)'   slope of distr. for Hakkinen-Mellor lateral melting     hakspl       = ', hakspl 
983         WRITE(numout,*)'   slope of distribution for Hibler lateral melting        hibspl       = ', hibspl
984         WRITE(numout,*)'   exponent for leads-closure rate                         exld         = ', exld
985         WRITE(numout,*)'   coefficient for diffusions of ice and snow              hakdif       = ', hakdif
986         WRITE(numout,*)'   threshold thick. for comp. of eq. thermal conductivity  zhth         = ', thth 
987         WRITE(numout,*)'   thickness of the surf. layer in temp. computation       hnzst        = ', hnzst
988         WRITE(numout,*)'   switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub 
989         WRITE(numout,*)'   coefficient for snow density when snow ice formation    alphs        = ', alphs
990         WRITE(numout,*)'   coefficient for ice-lead partition of snowfall          betas        = ', betas
991         WRITE(numout,*)'   extinction radiation parameter in sea ice (1.0)         kappa_i      = ', kappa_i
992         WRITE(numout,*)'   maximal n. of iter. for heat diffusion computation      nconv_i_thd  = ', nconv_i_thd
993         WRITE(numout,*)'   maximal err. on T for heat diffusion computation        maxer_i_thd  = ', maxer_i_thd
994         WRITE(numout,*)'   switch for comp. of thermal conductivity in the ice     thcon_i_swi  = ', thcon_i_swi
995         WRITE(numout,*)
996      ENDIF
997
998      rcdsn = hakdif * rcdsn 
999      rcdic = hakdif * rcdic
1000
1001
1002   END SUBROUTINE lim_thd_init
1003
1004#else
1005   !!======================================================================
1006   !!                       ***  MODULE limthd   ***
1007   !!                        No sea ice model
1008   !!======================================================================
1009CONTAINS
1010   SUBROUTINE lim_thd         ! Empty routine
1011   END SUBROUTINE lim_thd
1012   SUBROUTINE lim_thd_con_dif
1013   END SUBROUTINE lim_thd_con_dif
1014#endif
1015
1016   !!======================================================================
1017END MODULE limthd
Note: See TracBrowser for help on using the repository browser.