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 branches/dev_003_CPL/NEMO/LIM_SRC_3 – NEMO

source: branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd.F90 @ 991

Last change on this file since 991 was 991, checked in by smasson, 16 years ago

dev_003_CPL: preliminary draft (not working), see ticket #155
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 2.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            IF ( lk_cpl ) THEN
317               CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb)     , qla_ice(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
318               CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb)     , dqla_ice(:,:,jl)   , jpi, jpj, npb(1:nbpb) )
319            ENDIF
320            CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb)     , dqns_ice(:,:,jl)   , jpi, jpj, npb(1:nbpb) )
321            CALL tab_2d_1d( nbpb, t_bo_b     (1:nbpb)     , t_bo       , jpi, jpj, npb(1:nbpb) )
322            CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb)     , sprecip    , jpi, jpj, npb(1:nbpb) ) 
323            CALL tab_2d_1d( nbpb, fbif_1d    (1:nbpb)     , fbif       , jpi, jpj, npb(1:nbpb) )
324            CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) )
325            CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) )
326            CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) )
327            CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) )
328
329            CALL tab_2d_1d( nbpb, fseqv_1d   (1:nbpb)     , fseqv      , jpi, jpj, npb(1:nbpb) )
330            CALL tab_2d_1d( nbpb, fsbri_1d   (1:nbpb)     , fsbri      , jpi, jpj, npb(1:nbpb) )
331            CALL tab_2d_1d( nbpb, fhbri_1d   (1:nbpb)     , fhbri      , jpi, jpj, npb(1:nbpb) )
332            CALL tab_2d_1d( nbpb, fstbif_1d  (1:nbpb)     , fstric     , jpi, jpj, npb(1:nbpb) )
333            CALL tab_2d_1d( nbpb, qfvbq_1d   (1:nbpb)     , qfvbq      , jpi, jpj, npb(1:nbpb) )
334
335            !--------------------------------
336            ! 4.3) Thermodynamic processes
337            !--------------------------------
338
339            IF ( con_i ) CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
340            IF ( con_i ) CALL lim_thd_glohec( qt_i_in , qt_s_in ,             &
341               q_i_layer_in , 1 , nbpb , jl )
342
343            !---------------------------------!
344            CALL lim_thd_dif(1,nbpb,jl)   ! Ice/Snow Temperature profile    !
345            !---------------------------------!
346
347            CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
348            ! compulsory for limthd_dh
349
350            IF ( con_i ) CALL lim_thd_glohec( qt_i_fin , qt_s_fin ,           &
351               q_i_layer_fin , 1 , nbpb , jl ) 
352            IF ( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl )
353
354            !---------------------------------!
355            CALL lim_thd_dh(1,nbpb,jl)    ! Ice/Snow thickness              !
356            !---------------------------------!
357
358            !---------------------------------!
359            CALL lim_thd_ent(1,nbpb,jl)   ! Ice/Snow enthalpy remapping     !
360            !---------------------------------!
361
362            !---------------------------------!
363            CALL lim_thd_sal(1,nbpb)      ! Ice salinity computation        !
364            !---------------------------------!
365
366            !           CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
367            IF ( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin,             &
368               q_i_layer_fin , 1 , nbpb , jl ) 
369            IF ( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl )
370
371            !--------------------------------
372            ! 4.4) Move 1D to 2D vectors
373            !--------------------------------
374
375            CALL tab_1d_2d( nbpb, at_i        , npb, at_i_b (1:nbpb), jpi, jpj )
376            CALL tab_1d_2d( nbpb, ht_i(:,:,jl), npb, ht_i_b(1:nbpb), jpi, jpj )
377            CALL tab_1d_2d( nbpb, ht_s(:,:,jl), npb, ht_s_b(1:nbpb), jpi, jpj )
378            CALL tab_1d_2d( nbpb, a_i (:,:,jl), npb, a_i_b(1:nbpb) , jpi, jpj )
379            CALL tab_1d_2d( nbpb, t_su(:,:,jl), npb, t_su_b(1:nbpb), jpi, jpj )
380            CALL tab_1d_2d( nbpb, sm_i(:,:,jl), npb, sm_i_b(1:nbpb), jpi, jpj )
381
382            DO jk = 1, nlay_s
383               CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b(1:nbpb,jk), jpi, jpj)
384               CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b(1:nbpb,jk), jpi, jpj)
385            END DO
386
387            DO jk = 1, nlay_i
388               CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b(1:nbpb,jk), jpi, jpj)
389               CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b(1:nbpb,jk), jpi, jpj)
390               CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b(1:nbpb,jk), jpi, jpj)
391            END DO
392
393            CALL tab_1d_2d( nbpb, fstric     , npb, fstbif_1d (1:nbpb)     , jpi, jpj )
394            CALL tab_1d_2d( nbpb, qldif      , npb, qldif_1d  (1:nbpb)     , jpi, jpj )
395            CALL tab_1d_2d( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj )
396            CALL tab_1d_2d( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj )
397            CALL tab_1d_2d( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj )
398            CALL tab_1d_2d( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj )
399            CALL tab_1d_2d( nbpb, rdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj )
400            CALL tab_1d_2d( nbpb, rdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj )
401            CALL tab_1d_2d( nbpb, fdvolif    , npb, dvlbq_1d  (1:nbpb)     , jpi, jpj )
402            CALL tab_1d_2d( nbpb, rdvonif    , npb, dvnbq_1d  (1:nbpb)     , jpi, jpj ) 
403            CALL tab_1d_2d( nbpb, fseqv      , npb, fseqv_1d  (1:nbpb)     , jpi, jpj )
404
405            IF ( num_sal .EQ. 2 ) THEN
406               CALL tab_1d_2d( nbpb, fsbri      , npb, fsbri_1d  (1:nbpb)     , jpi, jpj )
407               CALL tab_1d_2d( nbpb, fhbri      , npb, fhbri_1d  (1:nbpb)     , jpi, jpj )
408            ENDIF
409
410            !+++++
411            !temporary stuff for a dummyversion
412            CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb)      , jpi, jpj )
413            CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb)      , jpi, jpj )
414            CALL tab_1d_2d( nbpb, fsup2D     , npb, fsup     (1:nbpb)      , jpi, jpj )
415            CALL tab_1d_2d( nbpb, focea2D    , npb, focea    (1:nbpb)      , jpi, jpj )
416            CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new  (1:nbpb)      , jpi, jpj )
417            CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0    (1:nbpb)      , jpi, jpj )
418            CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj)
419            !+++++
420
421            IF( lk_mpp ) CALL mpp_comm_free(ncomm_ice) !RB necessary ??
422         ENDIF ! nbpb
423
424      END DO ! jl
425
426      !------------------------------------------------------------------------------!
427      ! 5) Global variables, diagnostics
428      !------------------------------------------------------------------------------!
429
430      !------------------------
431      ! 5.1) Ice heat content             
432      !------------------------
433
434      ! Enthalpies are global variables we have to readjust the units
435      DO jl = 1, jpl
436         DO jk = 1, nlay_i
437            DO jj = 1, jpj
438               DO ji = 1, jpi
439                  ! Change dimensions
440                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac
441
442                  ! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules
443                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &
444                     area(ji,jj) * a_i(ji,jj,jl) * &
445                     ht_i(ji,jj,jl) / nlay_i
446               END DO !ji
447            END DO !jj
448         END DO !jk
449      END DO !jl
450
451      !------------------------
452      ! 5.2) Snow heat content             
453      !------------------------
454
455      ! Enthalpies are global variables we have to readjust the units
456      DO jl = 1, jpl
457         DO jk = 1, nlay_s
458            DO jj = 1, jpj
459               DO ji = 1, jpi
460                  ! Change dimensions
461                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
462                  ! Multiply by volume, so that heat content in 10^9 Joules
463                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
464                     a_i(ji,jj,jl) * ht_s(ji,jj,jl)  / nlay_s
465               END DO !ji
466            END DO !jj
467         END DO !jk
468      END DO !jl
469
470      !----------------------------------
471      ! 5.3) Change thickness to volume
472      !----------------------------------
473      CALL lim_var_eqv2glo
474
475      !--------------------------------------------
476      ! 5.4) Diagnostic thermodynamic growth rates
477      !--------------------------------------------
478      d_v_i_thd (:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)    ! ice volumes
479      dv_dt_thd(:,:,:)    = d_v_i_thd(:,:,:) / rdt_ice * 86400.0
480
481      IF ( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
482
483      IF(ln_ctl) THEN   ! Control print
484         CALL prt_ctl_info(' ')
485         CALL prt_ctl_info(' - Cell values : ')
486         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
487         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd  : cell area :')
488         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd  : at_i      :')
489         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd  : vt_i      :')
490         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_thd  : vt_s      :')
491         DO jl = 1, jpl
492            CALL prt_ctl_info(' ')
493            CALL prt_ctl_info(' - Category : ', ivar1=jl)
494            CALL prt_ctl_info('   ~~~~~~~~~~')
495            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_thd  : a_i      : ')
496            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_thd  : ht_i     : ')
497            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_thd  : ht_s     : ')
498            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_thd  : v_i      : ')
499            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_thd  : v_s      : ')
500            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_thd  : e_s      : ')
501            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_thd  : t_su     : ')
502            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_thd  : t_snow   : ')
503            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_thd  : sm_i     : ')
504            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_thd  : smv_i    : ')
505            DO jk = 1, nlay_i
506               CALL prt_ctl_info(' ')
507               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
508               CALL prt_ctl_info('   ~~~~~~~')
509               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_thd  : t_i      : ')
510               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_thd  : e_i      : ')
511            END DO
512         END DO
513
514      ENDIF
515
516   END SUBROUTINE lim_thd
517
518   !===============================================================================
519
520   SUBROUTINE lim_thd_glohec(eti,ets,etilayer,kideb,kiut,jl)
521      !!-----------------------------------------------------------------------
522      !!                   ***  ROUTINE lim_thd_glohec ***
523      !!                 
524      !! ** Purpose :  Compute total heat content for each category
525      !!               Works with 1d vectors only
526      !!
527      !! history :
528      !!  9.9  ! 07-04 (M.Vancoppenolle) original code
529      !!-----------------------------------------------------------------------
530      !! * Local variables
531      INTEGER, INTENT(in) :: &
532         kideb, kiut,        &  ! bounds for the spatial loop
533         jl                     ! category number
534
535      REAL(wp), DIMENSION (jpij,jpl), INTENT(out) ::  &
536         eti, ets            ! vertically-summed heat content for ice /snow
537
538      REAL(wp), DIMENSION (jpij,jkmax), INTENT(out) ::  &
539         etilayer            ! heat content for ice layers
540
541      REAL(wp) :: &
542         zdes,    &          ! snow heat content increment (dummy)
543         zeps                ! very small value (1.e-10)
544
545      INTEGER  :: &
546         ji,jk               ! loop indices
547
548      !!-----------------------------------------------------------------------
549      eti(:,:) = 0.0
550      ets(:,:) = 0.0
551      zeps     = 1.0e-10
552
553      ! total q over all layers, ice [J.m-2]
554      DO jk = 1, nlay_i
555         DO ji = kideb, kiut
556            etilayer(ji,jk) = q_i_b(ji,jk) &
557               * ht_i_b(ji) / nlay_i
558            eti(ji,jl) = eti(ji,jl) + etilayer(ji,jk) 
559         END DO
560      END DO
561
562      ! total q over all layers, snow [J.m-2]
563      DO ji = kideb, kiut
564         zdes = q_s_b(ji,1) * ht_s_b(ji) / nlay_s 
565         ets(ji,jl) = ets(ji,jl) + zdes       
566      END DO
567
568      WRITE(numout,*) ' lim_thd_glohec '
569      WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice
570      WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice
571      WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) +         &
572         ets(jiindex_1d,jl) ) / rdt_ice
573
574   END SUBROUTINE lim_thd_glohec
575
576   !===============================================================================
577
578   SUBROUTINE lim_thd_con_dif(kideb,kiut,jl)
579      !!-----------------------------------------------------------------------
580      !!                   ***  ROUTINE lim_thd_con_dif ***
581      !!                 
582      !! ** Purpose :   Test energy conservation after heat diffusion
583      !!
584      !! history :
585      !!  9.9  ! 07-04 (M.Vancoppenolle) original code
586      !!-------------------------------------------------------------------
587      !! * Local variables
588      INTEGER, INTENT(in) ::        &
589         kideb, kiut,               &  !: bounds for the spatial loop
590         jl                            !: category number
591
592      REAL(wp)                 ::   &  !: ! goes to trash
593         meance,                    &  !: mean conservation error
594         max_cons_err,              &  !: maximum tolerated conservation error
595         max_surf_err                  !: maximum tolerated surface error
596
597      INTEGER ::                    &
598         numce                         !: number of points for which conservation
599      !  is violated
600      INTEGER  :: &
601         ji,jk,                     &  !: loop indices
602         zji, zjj
603      !!---------------------------------------------------------------------
604
605      max_cons_err =  1.0
606      max_surf_err =  0.001
607
608      !--------------------------
609      ! Increment of energy
610      !--------------------------
611      ! global
612      DO ji = kideb, kiut
613         dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl)  &
614            + qt_s_fin(ji,jl) - qt_s_in(ji,jl)
615      END DO
616      ! layer by layer
617      dq_i_layer(:,:)    = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
618
619      !----------------------------------------
620      ! Atmospheric heat flux, ice heat budget
621      !----------------------------------------
622
623      DO ji = kideb, kiut
624         zji                 = MOD( npb(ji) - 1, jpi ) + 1
625         zjj                 = ( npb(ji) - 1 ) / jpi + 1
626
627         fatm(ji,jl) = &
628            qnsr_ice_1d(ji)                  + & ! atm non solar
629            (1.0-i0(ji))*qsr_ice_1d(ji)          ! atm solar
630
631         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) &
632            - fstroc(zji,zjj,jl)
633      END DO
634
635      !--------------------
636      ! Conservation error
637      !--------------------
638
639      DO ji = kideb, kiut
640         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
641      END DO
642
643      numce = 0
644      meance = 0.0
645      DO ji = kideb, kiut
646         IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
647            numce = numce + 1
648            meance = meance + cons_error(ji,jl)
649         ENDIF
650      ENDDO
651      IF (numce .GT. 0 ) meance = meance / numce
652
653      WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
654      WRITE(numout,*) ' After lim_thd_dif, category : ', jl
655      WRITE(numout,*) ' Mean conservation error on big error points ', meance, &
656         numit
657      WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit
658
659      !-------------------------------------------------------
660      ! Surface error due to imbalance between Fatm and Fcsu
661      !-------------------------------------------------------
662      numce  = 0.0
663      meance = 0.0
664
665      DO ji = kideb, kiut
666         surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) )
667         IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. &
668            max_surf_err ) ) THEN
669            numce = numce + 1 
670            meance = meance + surf_error(ji,jl)
671         ENDIF
672      ENDDO
673      IF (numce .GT. 0 ) meance = meance / numce
674
675      WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err
676      WRITE(numout,*) ' After lim_thd_dif, category : ', jl
677      WRITE(numout,*) ' Mean surface error on big error points ', meance, numit
678      WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit
679
680      IF (jiindex_1D.GT.0) WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d)
681      IF (jiindex_1D.GT.0) WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl)
682      IF (jiindex_1D.GT.0) WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d)
683
684      !---------------------------------------
685      ! Write ice state in case of big errors
686      !---------------------------------------
687
688      DO ji = kideb, kiut
689         IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. &
690            ( cons_error(ji,jl) .GT. max_cons_err  ) ) THEN
691            zji                 = MOD( npb(ji) - 1, jpi ) + 1
692            zjj                 = ( npb(ji) - 1 ) / jpi + 1
693
694            WRITE(numout,*) ' alerte 1     '
695            WRITE(numout,*) ' Untolerated conservation / surface error after '
696            WRITE(numout,*) ' heat diffusion in the ice '
697            WRITE(numout,*) ' Category   : ', jl
698            WRITE(numout,*) ' zji , zjj  : ', zji, zjj
699            WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj)
700            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
701            WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)
702            WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) / rdt_ice
703            WRITE(numout,*) ' Fdt        : ', sum_fluxq(ji,jl)
704            WRITE(numout,*)
705            !        WRITE(numout,*) ' qt_i_in   : ', qt_i_in(ji,jl)
706            !        WRITE(numout,*) ' qt_s_in   : ', qt_s_in(ji,jl)
707            !        WRITE(numout,*) ' qt_i_fin  : ', qt_i_fin(ji,jl)
708            !        WRITE(numout,*) ' qt_s_fin  : ', qt_s_fin(ji,jl)
709            !        WRITE(numout,*) ' qt        : ', qt_i_fin(ji,jl) + &
710            !                                         qt_s_fin(ji,jl)
711            WRITE(numout,*) ' ht_i       : ', ht_i_b(ji)
712            WRITE(numout,*) ' ht_s       : ', ht_s_b(ji)
713            WRITE(numout,*) ' t_su       : ', t_su_b(ji)
714            WRITE(numout,*) ' t_s        : ', t_s_b(ji,1)
715            WRITE(numout,*) ' t_i        : ', t_i_b(ji,1:nlay_i)
716            WRITE(numout,*) ' t_bo       : ', t_bo_b(ji)
717            WRITE(numout,*) ' q_i        : ', q_i_b(ji,1:nlay_i)
718            WRITE(numout,*) ' s_i        : ', s_i_b(ji,1:nlay_i)
719            WRITE(numout,*) ' tmelts     : ', rtt - tmut*s_i_b(ji,1:nlay_i)
720            WRITE(numout,*)
721            WRITE(numout,*) ' Fluxes '
722            WRITE(numout,*) ' ~~~~~~ '
723            WRITE(numout,*) ' fatm       : ', fatm(ji,jl)
724            WRITE(numout,*) ' fc_su      : ', fc_su    (ji)
725            WRITE(numout,*) ' fstr_inice : ', qsr_ice_1d(ji)*i0(ji)
726            WRITE(numout,*) ' fc_bo      : ', - fc_bo_i  (ji)
727            WRITE(numout,*) ' foc        : ', fbif_1d(ji)
728            WRITE(numout,*) ' fstroc     : ', fstroc   (zji,zjj,jl)
729            WRITE(numout,*) ' i0         : ', i0(ji)
730            WRITE(numout,*) ' qsr_ice    : ', (1.0-i0(ji))*qsr_ice_1d(ji)
731            WRITE(numout,*) ' qns_ice    : ', qnsr_ice_1d(ji)
732            WRITE(numout,*) ' Conduction fluxes : '
733            WRITE(numout,*) ' fc_s      : ', fc_s(ji,0:nlay_s)
734            WRITE(numout,*) ' fc_i      : ', fc_i(ji,0:nlay_i)
735            WRITE(numout,*)
736            WRITE(numout,*) ' Layer by layer ... '
737            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - &
738               qt_s_in(ji,jl) )  & 
739               / rdt_ice
740            WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) -      &
741               fc_s(ji,0)
742            DO jk = 1, nlay_i
743               WRITE(numout,*) ' layer  : ', jk
744               WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice 
745               WRITE(numout,*) ' radab  : ', radab(ji,jk)
746               WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) -               &
747                  fc_i(ji,jk-1)
748               WRITE(numout,*) ' tot f  : ', fc_i(ji,jk) -               &
749                  fc_i(ji,jk-1) - radab(ji,jk)
750            END DO
751
752         ENDIF
753
754      END DO
755
756   END SUBROUTINE lim_thd_con_dif
757
758   !==============================================================================
759
760   SUBROUTINE lim_thd_con_dh(kideb,kiut,jl)
761      !!-----------------------------------------------------------------------
762      !!                   ***  ROUTINE lim_thd_con_dh  ***
763      !!                 
764      !! ** Purpose :   Test energy conservation after enthalpy redistr.
765      !!
766      !! history :
767      !!  9.9  ! 07-04 (M.Vancoppenolle) original code
768      !!-----------------------------------------------------------------------
769      !! * Local variables
770      INTEGER, INTENT(in) ::        &
771         kideb, kiut,               &  !: bounds for the spatial loop
772         jl                            !: category number
773
774      REAL(wp)                 ::   &  !: ! goes to trash
775         meance,                    &  !: mean conservation error
776         max_cons_err                  !: maximum tolerated conservation error
777
778      INTEGER ::                    &
779         numce                         !: number of points for which conservation
780      !  is violated
781      INTEGER  ::  ji, zji, zjj        ! loop indices
782      !!---------------------------------------------------------------------
783
784      max_cons_err = 1.0
785
786      !--------------------------
787      ! Increment of energy
788      !--------------------------
789      ! global
790      DO ji = kideb, kiut
791         dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl)  &
792            + qt_s_fin(ji,jl) - qt_s_in(ji,jl)
793      END DO
794      ! layer by layer
795      dq_i_layer(:,:)    = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
796
797      !----------------------------------------
798      ! Atmospheric heat flux, ice heat budget
799      !----------------------------------------
800
801      DO ji = kideb, kiut
802         zji                 = MOD( npb(ji) - 1, jpi ) + 1
803         zjj                 = ( npb(ji) - 1 ) / jpi + 1
804
805         fatm(ji,jl) = &
806            qnsr_ice_1d(ji)                  + & ! atm non solar
807            !         (1.0-i0(ji))*qsr_ice_1d(ji)          ! atm solar
808            qsr_ice_1d(ji)                       ! atm solar
809
810         sum_fluxq(ji,jl)     = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) & 
811            - fstroc(zji,zjj,jl) 
812         cons_error(ji,jl)   = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
813      END DO
814
815      !--------------------
816      ! Conservation error
817      !--------------------
818
819      DO ji = kideb, kiut
820         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
821      END DO
822
823      numce = 0
824      meance = 0.0
825      DO ji = kideb, kiut
826         IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
827            numce = numce + 1
828            meance = meance + cons_error(ji,jl)
829         ENDIF
830      ENDDO
831      IF (numce .GT. 0 ) meance = meance / numce
832
833      WRITE(numout,*) ' Error report - Category : ', jl
834      WRITE(numout,*) ' ~~~~~~~~~~~~ '
835      WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
836      WRITE(numout,*) ' After lim_thd_ent, category : ', jl
837      WRITE(numout,*) ' Mean conservation error on big error points ', meance, &
838         numit
839      WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit
840
841      !---------------------------------------
842      ! Write ice state in case of big errors
843      !---------------------------------------
844
845      DO ji = kideb, kiut
846         IF ( cons_error(ji,jl) .GT. max_cons_err  ) THEN
847            zji                 = MOD( npb(ji) - 1, jpi ) + 1
848            zjj                 = ( npb(ji) - 1 ) / jpi + 1
849
850            WRITE(numout,*) ' alerte 1 - category : ', jl
851            WRITE(numout,*) ' Untolerated conservation error after limthd_ent '
852            WRITE(numout,*) ' zji , zjj  : ', zji, zjj
853            WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj)
854            WRITE(numout,*) ' * '
855            WRITE(numout,*) ' Ftotal     : ', sum_fluxq(ji,jl)
856            WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) / rdt_ice
857            WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice
858            WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice
859            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
860            WRITE(numout,*) ' * '
861            WRITE(numout,*) ' Fluxes        --- : '
862            WRITE(numout,*) ' fatm       : ', fatm(ji,jl)
863            WRITE(numout,*) ' foce       : ', fbif_1d(ji)
864            WRITE(numout,*) ' fres       : ', ftotal_fin(ji)
865            WRITE(numout,*) ' fhbri      : ', fhbricat(zji,zjj,jl)
866            WRITE(numout,*) ' * '
867            WRITE(numout,*) ' Heat contents --- : '
868            WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) / rdt_ice
869            WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) / rdt_ice
870            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + &
871               qt_s_in(ji,jl) ) / rdt_ice
872            WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) / rdt_ice
873            WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) / rdt_ice
874            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + &
875               qt_s_fin(ji,jl) ) / rdt_ice
876            WRITE(numout,*) ' * '
877            WRITE(numout,*) ' Ice variables --- : '
878            WRITE(numout,*) ' ht_i       : ', ht_i_b(ji)
879            WRITE(numout,*) ' ht_s       : ', ht_s_b(ji)
880            WRITE(numout,*) ' dh_s_tot  : ', dh_s_tot(ji)
881            WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji)
882            WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
883            WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
884
885         ENDIF
886
887      END DO
888
889   END SUBROUTINE lim_thd_con_dh
890   !==============================================================================
891
892   SUBROUTINE lim_thd_enmelt(kideb,kiut)
893      !!-----------------------------------------------------------------------
894      !!                   ***  ROUTINE lim_thd_enmelt ***
895      !!                 
896      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3)
897      !!
898      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
899      !!
900      !! history : Martin Vancoppenolle, May 2007
901      !!-------------------------------------------------------------------
902      INTEGER, INTENT(in) ::        &
903         kideb, kiut                   !: bounds for the spatial loop
904
905      REAL(wp)                 ::   &  !: goes to trash
906         ztmelts               ,    &  !: sea ice freezing point in K
907         zeps 
908
909      INTEGER             ::        &
910         ji,                        &  !: spatial loop index
911         jk                            !: vertical index
912
913      !!-------------------------------------------------------------------
914      zeps = 1.0e-10
915
916      ! Sea ice energy of melting
917      DO jk = 1, nlay_i
918         DO ji = kideb, kiut
919            ztmelts      =   - tmut * s_i_b(ji,jk) + rtt 
920            q_i_b(ji,jk) =   rhoic*( cpic    * ( ztmelts - t_i_b(ji,jk) )  &
921               + lfus  * ( 1.0 - (ztmelts-rtt)/MIN((t_i_b(ji,jk)-rtt),-zeps) )  &
922               - rcp      * ( ztmelts-rtt  ) ) 
923         END DO !ji
924      END DO !jk
925
926      ! Snow energy of melting
927      DO jk = 1, nlay_s
928         DO ji = kideb,kiut
929            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
930         END DO !ji
931      END DO !jk
932
933   END SUBROUTINE lim_thd_enmelt
934
935   !==============================================================================
936
937   SUBROUTINE lim_thd_init
938
939      !!-----------------------------------------------------------------------
940      !!                   ***  ROUTINE lim_thd_init ***
941      !!                 
942      !! ** Purpose :   Physical constants and parameters linked to the ice
943      !!      thermodynamics
944      !!
945      !! ** Method  :   Read the namicethd namelist and check the ice-thermo
946      !!       parameter values called at the first timestep (nit000)
947      !!
948      !! ** input   :   Namelist namicether
949      !!
950      !! history :
951      !!  8.5  ! 03-08 (C. Ethe) original code
952      !!-------------------------------------------------------------------
953      NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, &
954         &                hicmin, hiclim, amax  ,    &
955         &                sbeta  , parlat, hakspl, hibspl, exld,  &
956         &                hakdif, hnzst  , thth  , parsub, alphs, betas, & 
957         &                kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi
958      !!-------------------------------------------------------------------
959
960      ! Define the initial parameters
961      ! -------------------------
962      REWIND( numnam_ice )
963      READ  ( numnam_ice , namicethd )
964      IF (lwp) THEN
965         WRITE(numout,*)
966         WRITE(numout,*)'lim_thd_init : ice parameters for ice thermodynamic computation '
967         WRITE(numout,*)'~~~~~~~~~~~~'
968         WRITE(numout,*)'   maximum melting at the bottom                           hmelt        = ', hmelt
969         WRITE(numout,*)'   ice thick. for lateral accretion in NH (SH)             hiccrit(1/2) = ', hiccrit
970         WRITE(numout,*)'   Frazil ice thickness as a function of wind or not       fraz_swi     = ', fraz_swi
971         WRITE(numout,*)'   Maximum proportion of frazil ice collecting at bottom   maxfrazb     = ', maxfrazb
972         WRITE(numout,*)'   Thresold relative drift speed for collection of frazil  vfrazb       = ', vfrazb
973         WRITE(numout,*)'   Squeezing coefficient for collection of frazil          Cfrazb       = ', Cfrazb
974         WRITE(numout,*)'   ice thick. corr. to max. energy stored in brine pocket  hicmin       = ', hicmin 
975         WRITE(numout,*)'   minimum ice thickness                                   hiclim       = ', hiclim 
976         WRITE(numout,*)'   maximum lead fraction                                   amax         = ', amax
977         WRITE(numout,*)'   numerical carac. of the scheme for diffusion in ice '
978         WRITE(numout,*)'   Cranck-Nicholson (=0.5), implicit (=1), explicit (=0)   sbeta        = ', sbeta
979         WRITE(numout,*)'   percentage of energy used for lateral ablation          parlat       = ', parlat
980         WRITE(numout,*)'   slope of distr. for Hakkinen-Mellor lateral melting     hakspl       = ', hakspl 
981         WRITE(numout,*)'   slope of distribution for Hibler lateral melting        hibspl       = ', hibspl
982         WRITE(numout,*)'   exponent for leads-closure rate                         exld         = ', exld
983         WRITE(numout,*)'   coefficient for diffusions of ice and snow              hakdif       = ', hakdif
984         WRITE(numout,*)'   threshold thick. for comp. of eq. thermal conductivity  zhth         = ', thth 
985         WRITE(numout,*)'   thickness of the surf. layer in temp. computation       hnzst        = ', hnzst
986         WRITE(numout,*)'   switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub 
987         WRITE(numout,*)'   coefficient for snow density when snow ice formation    alphs        = ', alphs
988         WRITE(numout,*)'   coefficient for ice-lead partition of snowfall          betas        = ', betas
989         WRITE(numout,*)'   extinction radiation parameter in sea ice (1.0)         kappa_i      = ', kappa_i
990         WRITE(numout,*)'   maximal n. of iter. for heat diffusion computation      nconv_i_thd  = ', nconv_i_thd
991         WRITE(numout,*)'   maximal err. on T for heat diffusion computation        maxer_i_thd  = ', maxer_i_thd
992         WRITE(numout,*)'   switch for comp. of thermal conductivity in the ice     thcon_i_swi  = ', thcon_i_swi
993         WRITE(numout,*)
994      ENDIF
995
996      rcdsn = hakdif * rcdsn 
997      rcdic = hakdif * rcdic
998
999
1000   END SUBROUTINE lim_thd_init
1001
1002#else
1003   !!======================================================================
1004   !!                       ***  MODULE limthd   ***
1005   !!                        No sea ice model
1006   !!======================================================================
1007CONTAINS
1008   SUBROUTINE lim_thd         ! Empty routine
1009   END SUBROUTINE lim_thd
1010   SUBROUTINE lim_thd_con_dif
1011   END SUBROUTINE lim_thd_con_dif
1012#endif
1013
1014   !!======================================================================
1015END MODULE limthd
Note: See TracBrowser for help on using the repository browser.