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, 7 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.