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

Last change on this file since 1171 was 1171, checked in by ctlod, 16 years ago

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