source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 4161

Last change on this file since 4161 was 4161, checked in by cetlod, 7 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

  • Property svn:keywords set to Id
File size: 32.9 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) ::   epsi20 = 1e-20_wp   ! constant values
39   REAL(wp) ::   epsi13 = 1e-13_wp   !
40   REAL(wp) ::   epsi11 = 1e-11_wp   !
41   REAL(wp) ::   epsi10 = 1e-10_wp   !
42   REAL(wp) ::   epsi06 = 1e-06_wp   !
43   REAL(wp) ::   epsi03 = 1e-03_wp   !
44   REAL(wp) ::   zzero  = 0._wp      !
45   REAL(wp) ::   zone   = 1._wp      !
46
47   !!----------------------------------------------------------------------
48   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE lim_thd_lac
55      !!-------------------------------------------------------------------
56      !!               ***   ROUTINE lim_thd_lac  ***
57      !! 
58      !! ** Purpose : Computation of the evolution of the ice thickness and
59      !!      concentration as a function of the heat balance in the leads.
60      !!      It is only used for lateral accretion
61      !!       
62      !! ** Method  : Ice is formed in the open water when ocean lose heat
63      !!      (heat budget of open water Bl is negative) .
64      !!      Computation of the increase of 1-A (ice concentration) fol-
65      !!      lowing the law :
66      !!      (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
67      !!       where - h0 is the thickness of ice created in the lead
68      !!             - a is a minimum fraction for leads
69      !!             - F is a monotonic non-increasing function defined as:
70      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
71      !!             - exld is the exponent closure rate (=2 default val.)
72      !!
73      !! ** Action : - Adjustment of snow and ice thicknesses and heat
74      !!                content in brine pockets
75      !!             - Updating ice internal temperature
76      !!             - Computation of variation of ice volume and mass
77      !!             - Computation of frldb after lateral accretion and
78      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)     
79      !!------------------------------------------------------------------------
80      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices
81      INTEGER ::   layer, nbpac     ! local integers
82      INTEGER ::   ii, ij, iter   !   -       -
83      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars
84      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
85      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
86      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
87      CHARACTER (len = 15) :: fieldid
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) )  )   !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               IF( ji == jiindx .AND. jj == jjindx )   jiindex_1d = nbpac
289            ENDIF
290         END DO
291      END DO
292
293      IF( ln_nicep )   WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
294
295      !------------------------------
296      ! Move from 2-D to 1-D vectors
297      !------------------------------
298      ! If ocean gains heat do nothing
299      ! 0therwise compute new ice formation
300
301      IF ( nbpac > 0 ) THEN
302
303         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
304         DO jl = 1, jpl
305            CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
306            CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
307            CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
308            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
309            DO jk = 1, nlay_i
310               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
311            END DO ! jk
312         END DO ! jl
313
314         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif  , jpi, jpj, npac(1:nbpac) )
315         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif  , jpi, jpj, npac(1:nbpac) )
316         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
317         CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) )
318         CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) )
319         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
320         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
321
322         !------------------------------------------------------------------------------!
323         ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
324         !------------------------------------------------------------------------------!
325
326         !----------------------
327         ! Thickness of new ice
328         !----------------------
329         DO ji = 1, nbpac
330            zh_newice(ji) = hiccrit(1)
331         END DO
332         IF( fraz_swi == 1.0 )   zh_newice(:) = hicol_b(:)
333
334         !----------------------
335         ! Salinity of new ice
336         !----------------------
337
338         SELECT CASE ( num_sal )
339         CASE ( 1 )                    ! Sice = constant
340            zs_newice(:) = bulk_sal
341         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
342            DO ji = 1, nbpac
343               ii =   MOD( npac(ji) - 1 , jpi ) + 1
344               ij =      ( npac(ji) - 1 ) / jpi + 1
345               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  )
346            END DO
347         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
348            zs_newice(:) =   2.3
349         END SELECT
350
351
352         !-------------------------
353         ! Heat content of new ice
354         !-------------------------
355         ! We assume that new ice is formed at the seawater freezing point
356         DO ji = 1, nbpac
357            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
358            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
359               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
360               &                       - rcp  *         ( ztmelts - rtt )  )
361            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    &
362               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus
363         END DO ! ji
364         !----------------
365         ! Age of new ice
366         !----------------
367         DO ji = 1, nbpac
368            zo_newice(ji) = 0._wp
369         END DO ! ji
370
371         !--------------------------
372         ! Open water energy budget
373         !--------------------------
374         DO ji = 1, nbpac
375            zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)     !<0
376         END DO ! ji
377
378         !-------------------
379         ! Volume of new ice
380         !-------------------
381         DO ji = 1, nbpac
382            zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
383
384            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
385            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
386            zdh_frazb(ji) =         zfrazb   * zv_newice(ji)
387            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
388         END DO
389
390         !------------------------------------
391         ! Diags for energy conservation test
392         !------------------------------------
393         DO ji = 1, nbpac
394            ii = MOD( npac(ji) - 1 , jpi ) + 1
395            ij =    ( npac(ji) - 1 ) / jpi + 1
396            !
397            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
398            !
399            vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume
400            et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy
401
402         END DO
403
404         ! keep new ice volume in memory
405         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
406
407         !-----------------
408         ! Area of new ice
409         !-----------------
410         DO ji = 1, nbpac
411            ii = MOD( npac(ji) - 1 , jpi ) + 1
412            ij =    ( npac(ji) - 1 ) / jpi + 1
413            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
414            diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
415         END DO !ji
416
417         !------------------------------------------------------------------------------!
418         ! 6) Redistribute new ice area and volume into ice categories                  !
419         !------------------------------------------------------------------------------!
420
421         !-----------------------------------------
422         ! Keep old ice areas and volume in memory
423         !-----------------------------------------
424         zv_old(:,:) = zv_i_ac(:,:) 
425         za_old(:,:) = za_i_ac(:,:)
426
427         !-------------------------------------------
428         ! Compute excessive new ice area and volume
429         !-------------------------------------------
430         ! If lateral ice growth gives an ice concentration gt 1, then
431         ! we keep the excessive volume in memory and attribute it later to bottom accretion
432         DO ji = 1, nbpac
433            IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN
434               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) )
435               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
436               za_newice(ji) = za_newice(ji) - zda_res  (ji)
437               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
438            ELSE
439               zda_res(ji) = 0._wp
440               zdv_res(ji) = 0._wp
441            ENDIF
442         END DO ! ji
443
444         !------------------------------------------------
445         ! Laterally redistribute new ice volume and area
446         !------------------------------------------------
447         zat_i_ac(:) = 0._wp
448         DO jl = 1, jpl
449            DO ji = 1, nbpac
450               IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   &
451                  & zh_newice(ji)    <=  hi_max   (jl)         ) THEN
452                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
453                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
454                  zat_i_ac(ji)    = zat_i_ac(ji)    + za_i_ac  (ji,jl)
455                  zcatac  (ji)    = jl
456               ENDIF
457            END DO
458         END DO
459
460         !----------------------------------
461         ! Heat content - lateral accretion
462         !----------------------------------
463         DO ji = 1, nbpac
464            jl = zcatac(ji)                                                           ! categroy in which new ice is put
465            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) )             ! zindb=1 if ice =0 otherwise
466            zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb  ! old ice thickness
467            zdhex    (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) )           ! difference in thickness
468            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! ice totally new in jl category
469         END DO
470
471         DO jk = 1, nlay_i
472            DO ji = 1, nbpac
473               jl = zcatac(ji)
474               zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
475               zalphai = MIN( zhice_old(ji,jl) * REAL( jk )     / REAL( nlay_i ), zh_newice(ji) )   &
476                  &    - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
477               ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji)                                     &
478                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / REAL( nlay_i )  &
479                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       &
480                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
481            END DO
482         END DO
483
484         !-----------------------------------------------
485         ! Add excessive volume of new ice at the bottom
486         !-----------------------------------------------
487         ! If the ice concentration exceeds 1, the remaining volume of new ice
488         ! is equally redistributed among all ice categories in which there is
489         ! ice
490
491         ! Fraction of level ice
492         jm = 1
493         zat_i_lev(:) = 0._wp
494
495         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
496            DO ji = 1, nbpac
497               zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 
498            END DO
499         END DO
500
501         IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
502         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
503            DO ji = 1, nbpac
504               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
505               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem
506               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) , epsi06 )
507            END DO
508         END DO
509         IF( ln_nicep )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
510
511         !---------------------------------
512         ! Heat content - bottom accretion
513         !---------------------------------
514         jm = 1
515         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
516            DO ji = 1, nbpac
517               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) )       ! zindb=1 if ice =0 otherwise
518               zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
519               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    &
520                  &             +  zindb * zdh_frazb(ji)                               ! frazil ice may coalesce
521               zdummy(ji,jl)    = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb      ! thickness of residual ice
522            END DO
523         END DO
524
525         ! old layers thicknesses and enthalpies
526         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
527            DO jk = 1, nlay_i
528               DO ji = 1, nbpac
529                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / REAL( nlay_i )
530                  zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
531               END DO
532            END DO
533         END DO
534!!gm ???  why the previous do loop  if ocerwriten by the following one ?
535         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
536            DO ji = 1, nbpac
537               zthick0(ji,nlay_i+1,jl) =  zdhicbot(ji,jl)
538               zqm0   (ji,nlay_i+1,jl) =  ze_newice(ji) * zdhicbot(ji,jl)
539            END DO ! ji
540         END DO ! jl
541
542         ! Redistributing energy on the new grid
543         ze_i_ac(:,:,:) = 0._wp
544         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
545            DO jk = 1, nlay_i
546               DO layer = 1, nlay_i + 1
547                  DO ji = 1, nbpac
548                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 
549                     ! Redistributing energy on the new grid
550                     zweight = MAX (  MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) )   &
551                        &    - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp )   &
552                        &    /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
553                     ze_i_ac(ji,jk,jl) =  ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 
554                  END DO ! ji
555               END DO ! layer
556            END DO ! jk
557         END DO ! jl
558
559         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
560            DO jk = 1, nlay_i
561               DO ji = 1, nbpac
562                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 
563                  ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl)   &
564                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
565               END DO
566            END DO
567         END DO
568
569         !------------
570         ! Update age
571         !------------
572         DO jl = 1, jpl
573            DO ji = 1, nbpac
574               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
575               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb   
576            END DO
577         END DO   
578
579         !-----------------
580         ! Update salinity
581         !-----------------
582         !clem IF(  num_sal == 2  ) THEN
583            DO jl = 1, jpl
584               DO ji = 1, nbpac
585                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
586                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
587                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
588               END DO
589            END DO   
590         !clem ENDIF
591
592         !--------------------------------
593         ! Update mass/salt fluxes (clem)
594         !--------------------------------
595         DO jl = 1, jpl
596            DO ji = 1, nbpac
597               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
598               zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
599               rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
600               sfx_thd_1d(ji)   =   sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
601           END DO
602         END DO
603
604         !------------------------------------------------------------------------------!
605         ! 8) Change 2D vectors to 1D vectors
606         !------------------------------------------------------------------------------!
607         DO jl = 1, jpl
608            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
609            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
610            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
611            !clem IF (  num_sal == 2  )   &
612               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
613            DO jk = 1, nlay_i
614               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
615            END DO
616         END DO
617         CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
618         CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
619         !
620      ENDIF ! nbpac > 0
621
622      !------------------------------------------------------------------------------!
623      ! 9) Change units for e_i
624      !------------------------------------------------------------------------------!   
625      DO jl = 1, jpl
626         DO jk = 1, nlay_i          ! heat content in 10^9 Joules
627            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i )  / unit_fac 
628         END DO
629      END DO
630
631      !------------------------------------------------------------------------------|
632      ! 10) Conservation check and changes in each ice category
633      !------------------------------------------------------------------------------|
634      IF( con_i ) THEN
635         CALL lim_column_sum (jpl,   v_i, vt_i_final)
636         fieldid = 'v_i, limthd_lac'
637         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
638         !
639         CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
640         fieldid = 'e_i, limthd_lac'
641         CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid) 
642         !
643         CALL lim_column_sum (jpl,   v_s, vt_s_final)
644         fieldid = 'v_s, limthd_lac'
645         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
646         !
647         !     CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
648         !     fieldid = 'e_s, limthd_lac'
649         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
650         IF( ln_nicep ) THEN
651            WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindx,jjindx)
652            WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindx,jjindx)
653            WRITE(numout,*) ' et_i_init : ', et_i_init(jiindx,jjindx)
654            WRITE(numout,*) ' et_i_final: ', et_i_final(jiindx,jjindx)
655         ENDIF
656         !
657      ENDIF
658      !
659      CALL wrk_dealloc( jpij, zcatac )   ! integer
660      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
661      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
662      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 )
663      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
664      CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
665      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 )
666      !
667   END SUBROUTINE lim_thd_lac
668
669#else
670   !!----------------------------------------------------------------------
671   !!   Default option                               NO  LIM3 sea-ice model
672   !!----------------------------------------------------------------------
673CONTAINS
674   SUBROUTINE lim_thd_lac           ! Empty routine
675   END SUBROUTINE lim_thd_lac
676#endif
677
678   !!======================================================================
679END MODULE limthd_lac
Note: See TracBrowser for help on using the repository browser.