source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 7 years ago

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]

  • Property svn:keywords set to Id
File size: 34.6 KB
Line 
1MODULE limthd_lac
2   !!======================================================================
3   !!                       ***  MODULE limthd_lac   ***
4   !!                lateral thermodynamic growth of the ice
5   !!======================================================================
6   !! History :  LIM  ! 2005-12 (M. Vancoppenolle)  Original code
7   !!             -   ! 2006-01 (M. Vancoppenolle)  add ITD
8   !!            3.0  ! 2007-07 (M. Vancoppenolle)  Mass and energy conservation tested
9   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
10   !!----------------------------------------------------------------------
11#if defined key_lim3
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'                                      LIM3 sea-ice model
14   !!----------------------------------------------------------------------
15   !!   lim_lat_acr   : lateral accretion of ice
16   !!----------------------------------------------------------------------
17   USE par_oce        ! ocean parameters
18   USE dom_oce        ! domain variables
19   USE phycst         ! physical constants
20   USE sbc_oce        ! Surface boundary condition: ocean fields
21   USE sbc_ice        ! Surface boundary condition: ice fields
22   USE thd_ice        ! LIM thermodynamics
23   USE dom_ice        ! LIM domain
24   USE par_ice        ! LIM parameters
25   USE ice            ! LIM variables
26   USE limtab         ! LIM 2D <==> 1D
27   USE limcons        ! LIM conservation
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE wrk_nemo       ! work arrays
31   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC lim_thd_lac     ! called by lim_thd
37
38   REAL(wp) ::   epsi10 = 1.e-10_wp   !
39   REAL(wp) ::   zzero  = 0._wp      !
40   REAL(wp) ::   zone   = 1._wp      !
41
42   !!----------------------------------------------------------------------
43   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE lim_thd_lac
50      !!-------------------------------------------------------------------
51      !!               ***   ROUTINE lim_thd_lac  ***
52      !! 
53      !! ** Purpose : Computation of the evolution of the ice thickness and
54      !!      concentration as a function of the heat balance in the leads.
55      !!      It is only used for lateral accretion
56      !!       
57      !! ** Method  : Ice is formed in the open water when ocean lose heat
58      !!      (heat budget of open water Bl is negative) .
59      !!      Computation of the increase of 1-A (ice concentration) fol-
60      !!      lowing the law :
61      !!      (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
62      !!       where - h0 is the thickness of ice created in the lead
63      !!             - a is a minimum fraction for leads
64      !!             - F is a monotonic non-increasing function defined as:
65      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
66      !!             - exld is the exponent closure rate (=2 default val.)
67      !!
68      !! ** Action : - Adjustment of snow and ice thicknesses and heat
69      !!                content in brine pockets
70      !!             - Updating ice internal temperature
71      !!             - Computation of variation of ice volume and mass
72      !!             - Computation of frldb after lateral accretion and
73      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)     
74      !!------------------------------------------------------------------------
75      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices
76      INTEGER ::   layer, nbpac     ! local integers
77      INTEGER ::   ii, ij, iter   !   -       -
78      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars
79      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
80      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
81      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
82      CHARACTER (len = 15) :: fieldid
83
84      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
85      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
86      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
87      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
88       
89      INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows
90      REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
91
92      REAL(wp), POINTER, DIMENSION(:) ::   zv_newice   ! volume of accreted ice
93      REAL(wp), POINTER, DIMENSION(:) ::   za_newice   ! fractional area of accreted ice
94      REAL(wp), POINTER, DIMENSION(:) ::   zh_newice   ! thickness of accreted ice
95      REAL(wp), POINTER, DIMENSION(:) ::   ze_newice   ! heat content of accreted ice
96      REAL(wp), POINTER, DIMENSION(:) ::   zs_newice   ! salinity of accreted ice
97      REAL(wp), POINTER, DIMENSION(:) ::   zo_newice   ! age of accreted ice
98      REAL(wp), POINTER, DIMENSION(:) ::   zdv_res     ! residual volume in case of excessive heat budget
99      REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget
100      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_ac    ! total ice fraction   
101      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)   
102      REAL(wp), POINTER, DIMENSION(:) ::   zdh_frazb   ! accretion of frazil ice at the ice bottom
103      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_ac    ! relative ice / frazil velocity (1D vector)
104
105      REAL(wp), POINTER, DIMENSION(:,:) ::   zhice_old   ! previous ice thickness
106      REAL(wp), POINTER, DIMENSION(:,:) ::   zdummy      ! dummy thickness of new ice
107      REAL(wp), POINTER, DIMENSION(:,:) ::   zdhicbot    ! thickness of new ice which is accreted vertically
108      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl
109      REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl
110      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_ac     ! 1-D version of a_i
111      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_ac     ! 1-D version of v_i
112      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_ac    ! 1-D version of oa_i
113      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_ac   ! 1-D version of smv_i
114
115      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_ac   !: 1-D version of e_i
116
117      REAL(wp), POINTER, DIMENSION(:) ::   zqbgow    ! heat budget of the open water (negative)
118      REAL(wp), POINTER, DIMENSION(:) ::   zdhex     ! excessively thick accreted sea ice (hlead-hice)
119
120      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqm0      ! old layer-system heat content
121      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zthick0   ! old ice thickness
122
123      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
124      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   !  snow volume summed over categories
125      REAL(wp), POINTER, DIMENSION(:,:) ::   et_i_init, et_i_final   !  ice energy summed over categories
126      REAL(wp), POINTER, DIMENSION(:,:) ::   et_s_init               !  snow energy summed over categories
127      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
128      !!-----------------------------------------------------------------------!
129
130      CALL wrk_alloc( jpij, zcatac )   ! integer
131      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
132      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
133      CALL wrk_alloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
134      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac )
135      CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
136      CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
137
138      et_i_init(:,:) = 0._wp
139      et_s_init(:,:) = 0._wp
140      vt_i_init(:,:) = 0._wp
141      vt_s_init(:,:) = 0._wp
142
143      !------------------------------------------------------------------------------!
144      ! 1) Conservation check and changes in each ice category
145      !------------------------------------------------------------------------------!
146      IF( con_i ) THEN
147         CALL lim_column_sum        ( jpl, v_i          , vt_i_init)
148         CALL lim_column_sum        ( jpl, v_s          , vt_s_init)
149         CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init)
150         CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init)
151      ENDIF
152
153      !------------------------------------------------------------------------------|
154      ! 2) Convert units for ice internal energy
155      !------------------------------------------------------------------------------|
156      DO jl = 1, jpl
157         DO jk = 1, nlay_i
158            DO jj = 1, jpj
159               DO ji = 1, jpi
160                  !Energy of melting q(S,T) [J.m-3]
161                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i )
162                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
163                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
164               END DO
165            END DO
166         END DO
167      END DO
168
169      !------------------------------------------------------------------------------!
170      ! 3) Collection thickness of ice formed in leads and polynyas
171      !------------------------------------------------------------------------------!   
172      ! hicol is the thickness of new ice formed in open water
173      ! hicol can be either prescribed (frazswi = 0)
174      ! or computed (frazswi = 1)
175      ! Frazil ice forms in open water, is transported by wind
176      ! accumulates at the edge of the consolidated ice edge
177      ! where it forms aggregates of a specific thickness called
178      ! collection thickness.
179
180      ! Note : the following algorithm currently breaks vectorization
181      !
182
183      zvrel(:,:) = 0._wp
184
185      ! Default new ice thickness
186      hicol(:,:) = hiccrit(1)
187
188      IF( fraz_swi == 1._wp ) THEN
189
190         !--------------------
191         ! Physical constants
192         !--------------------
193         hicol(:,:) = 0._wp
194
195         zhicrit = 0.04 ! frazil ice thickness
196         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
197         zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
198         zgamafr = 0.03
199
200         DO jj = 1, jpj
201            DO ji = 1, jpi
202
203               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
204                  !-------------
205                  ! Wind stress
206                  !-------------
207                  ! C-grid wind stress components
208                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   &
209                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp
210                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   &
211                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp
212                  ! Square root of wind stress
213                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
214
215                  !---------------------
216                  ! Frazil ice velocity
217                  !---------------------
218                  zvfrx         = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
219                  zvfry         = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
220
221                  !-------------------
222                  ! Pack ice velocity
223                  !-------------------
224                  ! C-grid ice velocity
225                  zindb = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
226                  zvgx  = zindb * (  u_ice(ji-1,jj  ) * tmu(ji-1,jj  )    &
227                     &             + u_ice(ji,jj    ) * tmu(ji  ,jj  )  ) * 0.5_wp
228                  zvgy  = zindb * (  v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)    &
229                     &             + v_ice(ji,jj    ) * tmv(ji  ,jj  )  ) * 0.5_wp
230
231                  !-----------------------------------
232                  ! Relative frazil/pack ice velocity
233                  !-----------------------------------
234                  ! absolute relative velocity
235                  zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
236                     &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
237                  zvrel(ji,jj)  = SQRT( zvrel2 )
238
239                  !---------------------
240                  ! Iterative procedure
241                  !---------------------
242                  hicol(ji,jj) = zhicrit + 0.1 
243                  hicol(ji,jj) = zhicrit +   hicol(ji,jj)    &
244                     &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2
245
246!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1)
247!!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01
248!!gm                therefore the 2 lines with hicol can be replaced by 1 line:
249!!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2
250!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop
251
252                  iter = 1
253                  iterate_frazil = .true.
254
255                  DO WHILE ( iter .LT. 100 .AND. iterate_frazil ) 
256                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
257                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
258                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
259                        - zhicrit * ztwogp * zvrel2
260                     zhicol_new = hicol(ji,jj) - zf/zfp
261                     hicol(ji,jj)   = zhicol_new
262
263                     iter = iter + 1
264
265                  END DO ! do while
266
267               ENDIF ! end of selection of pixels where ice forms
268
269            END DO ! loop on ji ends
270         END DO ! loop on jj ends
271
272      ENDIF ! End of computation of frazil ice collection thickness
273
274      !------------------------------------------------------------------------------!
275      ! 4) Identify grid points where new ice forms
276      !------------------------------------------------------------------------------!
277
278      !-------------------------------------
279      ! Select points for new ice formation
280      !-------------------------------------
281      ! This occurs if open water energy budget is negative
282      nbpac = 0
283      DO jj = 1, jpj
284         DO ji = 1, jpi
285            IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN
286               nbpac = nbpac + 1
287               npac( nbpac ) = (jj - 1) * jpi + ji
288            ENDIF
289         END DO
290      END DO
291
292      ! debug point to follow
293      jiindex_1d = 0
294      IF( ln_nicep ) THEN
295         DO ji = mi0(jiindx), mi1(jiindx)
296            DO jj = mj0(jjindx), mj1(jjindx)
297               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN
298                  jiindex_1d = (jj - 1) * jpi + ji
299               ENDIF
300            END DO
301         END DO
302      ENDIF
303   
304      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
305
306      !------------------------------
307      ! Move from 2-D to 1-D vectors
308      !------------------------------
309      ! If ocean gains heat do nothing
310      ! 0therwise compute new ice formation
311
312      IF ( nbpac > 0 ) THEN
313
314         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
315         DO jl = 1, jpl
316            CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
317            CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
318            CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
319            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
320            DO jk = 1, nlay_i
321               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
322            END DO ! jk
323         END DO ! jl
324
325         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif  , jpi, jpj, npac(1:nbpac) )
326         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif  , jpi, jpj, npac(1:nbpac) )
327         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
328         CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) )
329         CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) )
330         CALL tab_2d_1d( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice, jpi, jpj, npac(1:nbpac) )
331         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
332         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
333
334         !------------------------------------------------------------------------------!
335         ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
336         !------------------------------------------------------------------------------!
337
338         !----------------------
339         ! Thickness of new ice
340         !----------------------
341         DO ji = 1, nbpac
342            zh_newice(ji) = hiccrit(1)
343         END DO
344         IF( fraz_swi == 1.0 )   zh_newice(:) = hicol_b(:)
345
346         !----------------------
347         ! Salinity of new ice
348         !----------------------
349
350         SELECT CASE ( num_sal )
351         CASE ( 1 )                    ! Sice = constant
352            zs_newice(:) = bulk_sal
353         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
354            DO ji = 1, nbpac
355               ii =   MOD( npac(ji) - 1 , jpi ) + 1
356               ij =      ( npac(ji) - 1 ) / jpi + 1
357               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  )
358            END DO
359         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
360            zs_newice(:) =   2.3
361         END SELECT
362
363
364         !-------------------------
365         ! Heat content of new ice
366         !-------------------------
367         ! We assume that new ice is formed at the seawater freezing point
368         DO ji = 1, nbpac
369            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
370            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
371               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
372               &                       - rcp  *         ( ztmelts - rtt )  )
373            ! MV HC 2014 comment I dont see why this line below is here... ?
374            ! This implies that ze_newice gets to rhoic*Lfus if it was negative, but this should never happen
375            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    &
376               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus
377         END DO ! ji
378         !----------------
379         ! Age of new ice
380         !----------------
381         DO ji = 1, nbpac
382            zo_newice(ji) = 0._wp
383         END DO ! ji
384
385         !--------------------------
386         ! Open water energy budget
387         !--------------------------
388         DO ji = 1, nbpac
389            zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)     !<0
390         END DO ! ji
391
392         !-------------------
393         ! Volume of new ice
394         !-------------------
395         DO ji = 1, nbpac
396
397            zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg]
398
399            zEw           = rcp * ( t_bo_b(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_b [J/kg]
400
401            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
402
403            zfmdt         = - zqbgow(ji) / zdE                     ! Fm.dt [kg/m2] (<0)
404
405            zv_newice(ji) = - zfmdt / rhoic
406
407            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux 
408
409            ! Contribution to energy flux to the ocean [J/m2], >0 
410            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm 
411
412            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
413            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
414            zdh_frazb(ji) =         zfrazb   * zv_newice(ji)
415            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
416         END DO
417
418         !------------------------------------
419         ! Diags for energy conservation test
420         !------------------------------------
421         DO ji = 1, nbpac
422            ii = MOD( npac(ji) - 1 , jpi ) + 1
423            ij =    ( npac(ji) - 1 ) / jpi + 1
424            !
425            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
426            !
427            vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume
428            et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy
429
430         END DO
431
432         ! keep new ice volume in memory
433         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
434
435         !-----------------
436         ! Area of new ice
437         !-----------------
438         DO ji = 1, nbpac
439            ii = MOD( npac(ji) - 1 , jpi ) + 1
440            ij =    ( npac(ji) - 1 ) / jpi + 1
441            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
442            diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
443         END DO !ji
444
445         !------------------------------------------------------------------------------!
446         ! 6) Redistribute new ice area and volume into ice categories                  !
447         !------------------------------------------------------------------------------!
448
449         !-----------------------------------------
450         ! Keep old ice areas and volume in memory
451         !-----------------------------------------
452         zv_old(:,:) = zv_i_ac(:,:) 
453         za_old(:,:) = za_i_ac(:,:)
454
455         !-------------------------------------------
456         ! Compute excessive new ice area and volume
457         !-------------------------------------------
458         ! If lateral ice growth gives an ice concentration gt 1, then
459         ! we keep the excessive volume in memory and attribute it later to bottom accretion
460         DO ji = 1, nbpac
461            IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN
462               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) )
463               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
464               za_newice(ji) = za_newice(ji) - zda_res  (ji)
465               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
466            ELSE
467               zda_res(ji) = 0._wp
468               zdv_res(ji) = 0._wp
469            ENDIF
470         END DO ! ji
471
472         !------------------------------------------------
473         ! Laterally redistribute new ice volume and area
474         !------------------------------------------------
475         zat_i_ac(:) = 0._wp
476         DO jl = 1, jpl
477            DO ji = 1, nbpac
478               IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   &
479                  & zh_newice(ji)    <=  hi_max   (jl)         ) THEN
480                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
481                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
482                  zat_i_ac(ji)    = zat_i_ac(ji)    + za_i_ac  (ji,jl)
483                  zcatac  (ji)    = jl
484               ENDIF
485            END DO
486         END DO
487
488         !----------------------------------
489         ! Heat content - lateral accretion
490         !----------------------------------
491         DO ji = 1, nbpac
492            jl = zcatac(ji)                                                           ! categroy in which new ice is put
493            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) )             ! zindb=1 if ice =0 otherwise
494            zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb  ! old ice thickness
495            zdhex    (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) )           ! difference in thickness
496            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! ice totally new in jl category
497         END DO
498
499         DO jk = 1, nlay_i
500            DO ji = 1, nbpac
501               jl = zcatac(ji)
502               zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
503               zalphai = MIN( zhice_old(ji,jl) * REAL( jk )     / REAL( nlay_i ), zh_newice(ji) )   &
504                  &    - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
505               ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji)                                     &
506                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / REAL( nlay_i )  &
507                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       &
508                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
509            END DO
510         END DO
511
512         !-----------------------------------------------
513         ! Add excessive volume of new ice at the bottom
514         !-----------------------------------------------
515         ! If the ice concentration exceeds 1, the remaining volume of new ice
516         ! is equally redistributed among all ice categories in which there is
517         ! ice
518
519         ! Fraction of level ice
520         jm = 1
521         zat_i_lev(:) = 0._wp
522
523         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
524            DO ji = 1, nbpac
525               zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 
526            END DO
527         END DO
528
529         IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
530         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
531            DO ji = 1, nbpac
532               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
533               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) )  ! clem
534               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 )
535            END DO
536         END DO
537         IF( ln_nicep .AND. jiindex_1d > 0 )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
538
539         !---------------------------------
540         ! Heat content - bottom accretion
541         !---------------------------------
542         jm = 1
543         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
544            DO ji = 1, nbpac
545               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) )       ! zindb=1 if ice =0 otherwise
546               zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
547               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    &
548                  &             +  zindb * zdh_frazb(ji)                               ! frazil ice may coalesce
549               zdummy(ji,jl)    = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb      ! thickness of residual ice
550            END DO
551         END DO
552
553         ! old layers thicknesses and enthalpies
554         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
555            DO jk = 1, nlay_i
556               DO ji = 1, nbpac
557                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / REAL( nlay_i )
558                  zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
559               END DO
560            END DO
561         END DO
562!!gm ???  why the previous do loop  if ocerwriten by the following one ?
563!! MV HC 2014 just because it'snot the same index and the expression is different
564         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
565            DO ji = 1, nbpac
566               zthick0(ji,nlay_i+1,jl) =  zdhicbot(ji,jl)
567               zqm0   (ji,nlay_i+1,jl) =  ze_newice(ji) * zdhicbot(ji,jl)
568            END DO ! ji
569         END DO ! jl
570
571         ! Redistributing energy on the new grid
572         ze_i_ac(:,:,:) = 0._wp
573         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
574            DO jk = 1, nlay_i
575               DO layer = 1, nlay_i + 1
576                  DO ji = 1, nbpac
577                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 
578                     ! Redistributing energy on the new grid
579                     zweight = MAX (  MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) )   &
580                        &    - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp )   &
581                        &    /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
582                     ze_i_ac(ji,jk,jl) =  ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 
583                  END DO ! ji
584               END DO ! layer
585            END DO ! jk
586         END DO ! jl
587
588         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
589            DO jk = 1, nlay_i
590               DO ji = 1, nbpac
591                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 
592                  ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl)   &
593                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
594               END DO
595            END DO
596         END DO
597
598         !------------
599         ! Update age
600         !------------
601         DO jl = 1, jpl
602            DO ji = 1, nbpac
603               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
604               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb   
605            END DO
606         END DO   
607
608         !-----------------
609         ! Update salinity
610         !-----------------
611         !clem IF(  num_sal == 2  ) THEN
612            DO jl = 1, jpl
613               DO ji = 1, nbpac
614                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
615                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
616                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
617               END DO
618            END DO   
619         !clem ENDIF
620
621         !--------------------------------
622         ! Update mass/salt fluxes (clem)
623         !--------------------------------
624         DO jl = 1, jpl
625            DO ji = 1, nbpac
626               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
627               zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
628               rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
629               sfx_thd_1d(ji)   =   sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
630           END DO
631         END DO
632
633         !------------------------------------------------------------------------------!
634         ! 8) Change 2D vectors to 1D vectors
635         !------------------------------------------------------------------------------!
636         DO jl = 1, jpl
637            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
638            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
639            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
640            !clem IF (  num_sal == 2  )   &
641               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
642            DO jk = 1, nlay_i
643               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
644            END DO
645         END DO
646         CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
647         CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
648         CALL tab_1d_2d( nbpac, rdq_ice, npac(1:nbpac), rdq_ice_1d(1:nbpac), jpi, jpj )
649         !
650      ENDIF ! nbpac > 0
651
652      !------------------------------------------------------------------------------!
653      ! 9) Change units for e_i
654      !------------------------------------------------------------------------------!   
655      DO jl = 1, jpl
656         DO jk = 1, nlay_i          ! heat content in 10^9 Joules
657            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i )  / unit_fac 
658         END DO
659      END DO
660
661      !------------------------------------------------------------------------------|
662      ! 10) Conservation check and changes in each ice category
663      !------------------------------------------------------------------------------|
664      IF( con_i ) THEN
665         CALL lim_column_sum (jpl,   v_i, vt_i_final)
666         fieldid = 'v_i, limthd_lac'
667         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
668         !
669         CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
670         fieldid = 'e_i, limthd_lac'
671         CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid) 
672         !
673         CALL lim_column_sum (jpl,   v_s, vt_s_final)
674         fieldid = 'v_s, limthd_lac'
675         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
676         !
677         !     CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
678         !     fieldid = 'e_s, limthd_lac'
679         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
680         IF( ln_nicep ) THEN
681            DO ji = mi0(jiindx), mi1(jiindx)
682               DO jj = mj0(jjindx), mj1(jjindx)
683                  WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
684                  WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
685                  WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
686                  WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
687               END DO
688            END DO
689         ENDIF
690         !
691      ENDIF
692      !
693      CALL wrk_dealloc( jpij, zcatac )   ! integer
694      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
695      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
696      CALL wrk_dealloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
697      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
698      CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
699      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
700      !
701   END SUBROUTINE lim_thd_lac
702
703#else
704   !!----------------------------------------------------------------------
705   !!   Default option                               NO  LIM3 sea-ice model
706   !!----------------------------------------------------------------------
707CONTAINS
708   SUBROUTINE lim_thd_lac           ! Empty routine
709   END SUBROUTINE lim_thd_lac
710#endif
711
712   !!======================================================================
713END MODULE limthd_lac
Note: See TracBrowser for help on using the repository browser.