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

Last change on this file since 1571 was 1571, checked in by ctlod, 15 years ago

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