New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limthd.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config

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