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/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 3625

Last change on this file since 3625 was 3625, checked in by acc, 11 years ago

Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF

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