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_lac.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_lac.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: 33.2 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      INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows
85      REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
86
87      REAL(wp), POINTER, DIMENSION(:) ::   zv_newice   ! volume of accreted ice
88      REAL(wp), POINTER, DIMENSION(:) ::   za_newice   ! fractional area of accreted ice
89      REAL(wp), POINTER, DIMENSION(:) ::   zh_newice   ! thickness of accreted ice
90      REAL(wp), POINTER, DIMENSION(:) ::   ze_newice   ! heat content of accreted ice
91      REAL(wp), POINTER, DIMENSION(:) ::   zs_newice   ! salinity of accreted ice
92      REAL(wp), POINTER, DIMENSION(:) ::   zo_newice   ! age of accreted ice
93      REAL(wp), POINTER, DIMENSION(:) ::   zdv_res     ! residual volume in case of excessive heat budget
94      REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget
95      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_ac    ! total ice fraction   
96      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)   
97      REAL(wp), POINTER, DIMENSION(:) ::   zdh_frazb   ! accretion of frazil ice at the ice bottom
98      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_ac    ! relative ice / frazil velocity (1D vector)
99
100      REAL(wp), POINTER, DIMENSION(:,:) ::   zhice_old   ! previous ice thickness
101      REAL(wp), POINTER, DIMENSION(:,:) ::   zdummy      ! dummy thickness of new ice
102      REAL(wp), POINTER, DIMENSION(:,:) ::   zdhicbot    ! thickness of new ice which is accreted vertically
103      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl
104      REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl
105      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_ac     ! 1-D version of a_i
106      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_ac     ! 1-D version of v_i
107      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_ac    ! 1-D version of oa_i
108      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_ac   ! 1-D version of smv_i
109
110      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_ac   !: 1-D version of e_i
111
112      REAL(wp), POINTER, DIMENSION(:) ::   zqbgow    ! heat budget of the open water (negative)
113      REAL(wp), POINTER, DIMENSION(:) ::   zdhex     ! excessively thick accreted sea ice (hlead-hice)
114
115      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqm0      ! old layer-system heat content
116      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zthick0   ! old ice thickness
117
118      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
119      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   !  snow volume summed over categories
120      REAL(wp), POINTER, DIMENSION(:,:) ::   et_i_init, et_i_final   !  ice energy summed over categories
121      REAL(wp), POINTER, DIMENSION(:,:) ::   et_s_init               !  snow energy summed over categories
122      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
123      !!-----------------------------------------------------------------------!
124
125      CALL wrk_alloc( jpij, zcatac )   ! integer
126      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
127      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
128      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 )
129      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac )
130      CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
131      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 )
132
133      et_i_init(:,:) = 0._wp
134      et_s_init(:,:) = 0._wp
135      vt_i_init(:,:) = 0._wp
136      vt_s_init(:,:) = 0._wp
137
138      !------------------------------------------------------------------------------!
139      ! 1) Conservation check and changes in each ice category
140      !------------------------------------------------------------------------------!
141      IF( con_i ) THEN
142         CALL lim_column_sum        ( jpl, v_i          , vt_i_init)
143         CALL lim_column_sum        ( jpl, v_s          , vt_s_init)
144         CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init)
145         CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init)
146      ENDIF
147
148      !------------------------------------------------------------------------------|
149      ! 2) Convert units for ice internal energy
150      !------------------------------------------------------------------------------|
151      DO jl = 1, jpl
152         DO jk = 1, nlay_i
153            DO jj = 1, jpj
154               DO ji = 1, jpi
155                  !Energy of melting q(S,T) [J.m-3]
156                  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 )
157                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
158                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
159               END DO
160            END DO
161         END DO
162      END DO
163
164      !------------------------------------------------------------------------------!
165      ! 3) Collection thickness of ice formed in leads and polynyas
166      !------------------------------------------------------------------------------!   
167      ! hicol is the thickness of new ice formed in open water
168      ! hicol can be either prescribed (frazswi = 0)
169      ! or computed (frazswi = 1)
170      ! Frazil ice forms in open water, is transported by wind
171      ! accumulates at the edge of the consolidated ice edge
172      ! where it forms aggregates of a specific thickness called
173      ! collection thickness.
174
175      ! Note : the following algorithm currently breaks vectorization
176      !
177
178      zvrel(:,:) = 0._wp
179
180      ! Default new ice thickness
181      hicol(:,:) = hiccrit(1)
182
183      IF( fraz_swi == 1._wp ) THEN
184
185         !--------------------
186         ! Physical constants
187         !--------------------
188         hicol(:,:) = 0._wp
189
190         zhicrit = 0.04 ! frazil ice thickness
191         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
192         zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
193         zgamafr = 0.03
194
195         DO jj = 1, jpj
196            DO ji = 1, jpi
197
198               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
199                  !-------------
200                  ! Wind stress
201                  !-------------
202                  ! C-grid wind stress components
203                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   &
204                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp
205                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   &
206                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp
207                  ! Square root of wind stress
208                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
209
210                  !---------------------
211                  ! Frazil ice velocity
212                  !---------------------
213                  zvfrx         = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
214                  zvfry         = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
215
216                  !-------------------
217                  ! Pack ice velocity
218                  !-------------------
219                  ! C-grid ice velocity
220                  zindb = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
221                  zvgx  = zindb * (  u_ice(ji-1,jj  ) * tmu(ji-1,jj  )    &
222                     &             + u_ice(ji,jj    ) * tmu(ji  ,jj  )  ) * 0.5_wp
223                  zvgy  = zindb * (  v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)    &
224                     &             + v_ice(ji,jj    ) * tmv(ji  ,jj  )  ) * 0.5_wp
225
226                  !-----------------------------------
227                  ! Relative frazil/pack ice velocity
228                  !-----------------------------------
229                  ! absolute relative velocity
230                  zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
231                     &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
232                  zvrel(ji,jj)  = SQRT( zvrel2 )
233
234                  !---------------------
235                  ! Iterative procedure
236                  !---------------------
237                  hicol(ji,jj) = zhicrit + 0.1 
238                  hicol(ji,jj) = zhicrit +   hicol(ji,jj)    &
239                     &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2
240
241!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1)
242!!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01
243!!gm                therefore the 2 lines with hicol can be replaced by 1 line:
244!!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2
245!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop
246
247                  iter = 1
248                  iterate_frazil = .true.
249
250                  DO WHILE ( iter .LT. 100 .AND. iterate_frazil ) 
251                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
252                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
253                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
254                        - zhicrit * ztwogp * zvrel2
255                     zhicol_new = hicol(ji,jj) - zf/zfp
256                     hicol(ji,jj)   = zhicol_new
257
258                     iter = iter + 1
259
260                  END DO ! do while
261
262               ENDIF ! end of selection of pixels where ice forms
263
264            END DO ! loop on ji ends
265         END DO ! loop on jj ends
266
267      ENDIF ! End of computation of frazil ice collection thickness
268
269      !------------------------------------------------------------------------------!
270      ! 4) Identify grid points where new ice forms
271      !------------------------------------------------------------------------------!
272
273      !-------------------------------------
274      ! Select points for new ice formation
275      !-------------------------------------
276      ! This occurs if open water energy budget is negative
277      nbpac = 0
278      DO jj = 1, jpj
279         DO ji = 1, jpi
280            IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN
281               nbpac = nbpac + 1
282               npac( nbpac ) = (jj - 1) * jpi + ji
283            ENDIF
284         END DO
285      END DO
286
287      ! debug point to follow
288      jiindex_1d = 0
289      IF( ln_nicep ) THEN
290         DO ji = mi0(jiindx), mi1(jiindx)
291            DO jj = mj0(jjindx), mj1(jjindx)
292               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN
293                  jiindex_1d = (jj - 1) * jpi + ji
294               ENDIF
295            END DO
296         END DO
297      ENDIF
298   
299      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
300
301      !------------------------------
302      ! Move from 2-D to 1-D vectors
303      !------------------------------
304      ! If ocean gains heat do nothing
305      ! 0therwise compute new ice formation
306
307      IF ( nbpac > 0 ) THEN
308
309         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
310         DO jl = 1, jpl
311            CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
312            CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
313            CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
314            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
315            DO jk = 1, nlay_i
316               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
317            END DO ! jk
318         END DO ! jl
319
320         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif  , jpi, jpj, npac(1:nbpac) )
321         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif  , jpi, jpj, npac(1:nbpac) )
322         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
323         CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) )
324         CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) )
325         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
326         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
327
328         !------------------------------------------------------------------------------!
329         ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
330         !------------------------------------------------------------------------------!
331
332         !----------------------
333         ! Thickness of new ice
334         !----------------------
335         DO ji = 1, nbpac
336            zh_newice(ji) = hiccrit(1)
337         END DO
338         IF( fraz_swi == 1.0 )   zh_newice(:) = hicol_b(:)
339
340         !----------------------
341         ! Salinity of new ice
342         !----------------------
343
344         SELECT CASE ( num_sal )
345         CASE ( 1 )                    ! Sice = constant
346            zs_newice(:) = bulk_sal
347         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
348            DO ji = 1, nbpac
349               ii =   MOD( npac(ji) - 1 , jpi ) + 1
350               ij =      ( npac(ji) - 1 ) / jpi + 1
351               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  )
352            END DO
353         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
354            zs_newice(:) =   2.3
355         END SELECT
356
357
358         !-------------------------
359         ! Heat content of new ice
360         !-------------------------
361         ! We assume that new ice is formed at the seawater freezing point
362         DO ji = 1, nbpac
363            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
364            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
365               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
366               &                       - rcp  *         ( ztmelts - rtt )  )
367            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    &
368               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus
369         END DO ! ji
370         !----------------
371         ! Age of new ice
372         !----------------
373         DO ji = 1, nbpac
374            zo_newice(ji) = 0._wp
375         END DO ! ji
376
377         !--------------------------
378         ! Open water energy budget
379         !--------------------------
380         DO ji = 1, nbpac
381            zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)     !<0
382         END DO ! ji
383
384         !-------------------
385         ! Volume of new ice
386         !-------------------
387         DO ji = 1, nbpac
388            zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
389
390            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
391            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
392            zdh_frazb(ji) =         zfrazb   * zv_newice(ji)
393            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
394         END DO
395
396         !------------------------------------
397         ! Diags for energy conservation test
398         !------------------------------------
399         DO ji = 1, nbpac
400            ii = MOD( npac(ji) - 1 , jpi ) + 1
401            ij =    ( npac(ji) - 1 ) / jpi + 1
402            !
403            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
404            !
405            vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume
406            et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy
407
408         END DO
409
410         ! keep new ice volume in memory
411         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
412
413         !-----------------
414         ! Area of new ice
415         !-----------------
416         DO ji = 1, nbpac
417            ii = MOD( npac(ji) - 1 , jpi ) + 1
418            ij =    ( npac(ji) - 1 ) / jpi + 1
419            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
420            diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
421         END DO !ji
422
423         !------------------------------------------------------------------------------!
424         ! 6) Redistribute new ice area and volume into ice categories                  !
425         !------------------------------------------------------------------------------!
426
427         !-----------------------------------------
428         ! Keep old ice areas and volume in memory
429         !-----------------------------------------
430         zv_old(:,:) = zv_i_ac(:,:) 
431         za_old(:,:) = za_i_ac(:,:)
432
433         !-------------------------------------------
434         ! Compute excessive new ice area and volume
435         !-------------------------------------------
436         ! If lateral ice growth gives an ice concentration gt 1, then
437         ! we keep the excessive volume in memory and attribute it later to bottom accretion
438         DO ji = 1, nbpac
439            IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN
440               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) )
441               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
442               za_newice(ji) = za_newice(ji) - zda_res  (ji)
443               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
444            ELSE
445               zda_res(ji) = 0._wp
446               zdv_res(ji) = 0._wp
447            ENDIF
448         END DO ! ji
449
450         !------------------------------------------------
451         ! Laterally redistribute new ice volume and area
452         !------------------------------------------------
453         zat_i_ac(:) = 0._wp
454         DO jl = 1, jpl
455            DO ji = 1, nbpac
456               IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   &
457                  & zh_newice(ji)    <=  hi_max   (jl)         ) THEN
458                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
459                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
460                  zat_i_ac(ji)    = zat_i_ac(ji)    + za_i_ac  (ji,jl)
461                  zcatac  (ji)    = jl
462               ENDIF
463            END DO
464         END DO
465
466         !----------------------------------
467         ! Heat content - lateral accretion
468         !----------------------------------
469         DO ji = 1, nbpac
470            jl = zcatac(ji)                                                           ! categroy in which new ice is put
471            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) )             ! zindb=1 if ice =0 otherwise
472            zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb  ! old ice thickness
473            zdhex    (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) )           ! difference in thickness
474            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! ice totally new in jl category
475         END DO
476
477         DO jk = 1, nlay_i
478            DO ji = 1, nbpac
479               jl = zcatac(ji)
480               zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
481               zalphai = MIN( zhice_old(ji,jl) * REAL( jk )     / REAL( nlay_i ), zh_newice(ji) )   &
482                  &    - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
483               ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji)                                     &
484                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / REAL( nlay_i )  &
485                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       &
486                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
487            END DO
488         END DO
489
490         !-----------------------------------------------
491         ! Add excessive volume of new ice at the bottom
492         !-----------------------------------------------
493         ! If the ice concentration exceeds 1, the remaining volume of new ice
494         ! is equally redistributed among all ice categories in which there is
495         ! ice
496
497         ! Fraction of level ice
498         jm = 1
499         zat_i_lev(:) = 0._wp
500
501         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
502            DO ji = 1, nbpac
503               zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 
504            END DO
505         END DO
506
507         IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
508         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
509            DO ji = 1, nbpac
510               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
511               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) )  ! clem
512               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 )
513            END DO
514         END DO
515         IF( ln_nicep .AND. jiindex_1d > 0 )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
516
517         !---------------------------------
518         ! Heat content - bottom accretion
519         !---------------------------------
520         jm = 1
521         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
522            DO ji = 1, nbpac
523               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) )       ! zindb=1 if ice =0 otherwise
524               zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
525               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    &
526                  &             +  zindb * zdh_frazb(ji)                               ! frazil ice may coalesce
527               zdummy(ji,jl)    = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb      ! thickness of residual ice
528            END DO
529         END DO
530
531         ! old layers thicknesses and enthalpies
532         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
533            DO jk = 1, nlay_i
534               DO ji = 1, nbpac
535                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / REAL( nlay_i )
536                  zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
537               END DO
538            END DO
539         END DO
540!!gm ???  why the previous do loop  if ocerwriten by the following one ?
541         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
542            DO ji = 1, nbpac
543               zthick0(ji,nlay_i+1,jl) =  zdhicbot(ji,jl)
544               zqm0   (ji,nlay_i+1,jl) =  ze_newice(ji) * zdhicbot(ji,jl)
545            END DO ! ji
546         END DO ! jl
547
548         ! Redistributing energy on the new grid
549         ze_i_ac(:,:,:) = 0._wp
550         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
551            DO jk = 1, nlay_i
552               DO layer = 1, nlay_i + 1
553                  DO ji = 1, nbpac
554                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 
555                     ! Redistributing energy on the new grid
556                     zweight = MAX (  MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) )   &
557                        &    - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp )   &
558                        &    /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
559                     ze_i_ac(ji,jk,jl) =  ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 
560                  END DO ! ji
561               END DO ! layer
562            END DO ! jk
563         END DO ! jl
564
565         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
566            DO jk = 1, nlay_i
567               DO ji = 1, nbpac
568                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 
569                  ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl)   &
570                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
571               END DO
572            END DO
573         END DO
574
575         !------------
576         ! Update age
577         !------------
578         DO jl = 1, jpl
579            DO ji = 1, nbpac
580               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
581               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb   
582            END DO
583         END DO   
584
585         !-----------------
586         ! Update salinity
587         !-----------------
588         !clem IF(  num_sal == 2  ) THEN
589            DO jl = 1, jpl
590               DO ji = 1, nbpac
591                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
592                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
593                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
594               END DO
595            END DO   
596         !clem ENDIF
597
598         !--------------------------------
599         ! Update mass/salt fluxes (clem)
600         !--------------------------------
601         DO jl = 1, jpl
602            DO ji = 1, nbpac
603               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
604               zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
605               rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
606               sfx_thd_1d(ji)   =   sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
607           END DO
608         END DO
609
610         !------------------------------------------------------------------------------!
611         ! 8) Change 2D vectors to 1D vectors
612         !------------------------------------------------------------------------------!
613         DO jl = 1, jpl
614            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
615            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
616            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
617            !clem IF (  num_sal == 2  )   &
618               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
619            DO jk = 1, nlay_i
620               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
621            END DO
622         END DO
623         CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
624         CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
625         !
626      ENDIF ! nbpac > 0
627
628      !------------------------------------------------------------------------------!
629      ! 9) Change units for e_i
630      !------------------------------------------------------------------------------!   
631      DO jl = 1, jpl
632         DO jk = 1, nlay_i          ! heat content in 10^9 Joules
633            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i )  / unit_fac 
634         END DO
635      END DO
636
637      !------------------------------------------------------------------------------|
638      ! 10) Conservation check and changes in each ice category
639      !------------------------------------------------------------------------------|
640      IF( con_i ) THEN
641         CALL lim_column_sum (jpl,   v_i, vt_i_final)
642         fieldid = 'v_i, limthd_lac'
643         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
644         !
645         CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
646         fieldid = 'e_i, limthd_lac'
647         CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid) 
648         !
649         CALL lim_column_sum (jpl,   v_s, vt_s_final)
650         fieldid = 'v_s, limthd_lac'
651         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
652         !
653         !     CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
654         !     fieldid = 'e_s, limthd_lac'
655         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
656         IF( ln_nicep ) THEN
657            DO ji = mi0(jiindx), mi1(jiindx)
658               DO jj = mj0(jjindx), mj1(jjindx)
659                  WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
660                  WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
661                  WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
662                  WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
663               END DO
664            END DO
665         ENDIF
666         !
667      ENDIF
668      !
669      CALL wrk_dealloc( jpij, zcatac )   ! integer
670      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
671      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
672      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 )
673      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
674      CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
675      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 )
676      !
677   END SUBROUTINE lim_thd_lac
678
679#else
680   !!----------------------------------------------------------------------
681   !!   Default option                               NO  LIM3 sea-ice model
682   !!----------------------------------------------------------------------
683CONTAINS
684   SUBROUTINE lim_thd_lac           ! Empty routine
685   END SUBROUTINE lim_thd_lac
686#endif
687
688   !!======================================================================
689END MODULE limthd_lac
Note: See TracBrowser for help on using the repository browser.