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

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 3517

Last change on this file since 3517 was 3517, checked in by gm, 10 years ago

gm: Branch: dev_r3385_NOCS04_HAMF; #665. update sbccpl ; change LIM3 from equivalent salt flux to salt flux and mass flux

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