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

source: branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 5087

Last change on this file since 5087 was 4972, checked in by clem, 10 years ago

LIM3 cleaning (1st part) with modification in the ice namelist. 2nd part is coming soon
  • Property svn:keywords set to Id
File size: 26.1 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 lbclnk         ! ocean lateral boundary conditions (or mpp link)
32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
33   USE limthd_ent
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC lim_thd_lac     ! called by lim_thd
39
40   !!----------------------------------------------------------------------
41   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE lim_thd_lac
48      !!-------------------------------------------------------------------
49      !!               ***   ROUTINE lim_thd_lac  ***
50      !! 
51      !! ** Purpose : Computation of the evolution of the ice thickness and
52      !!      concentration as a function of the heat balance in the leads.
53      !!      It is only used for lateral accretion
54      !!       
55      !! ** Method  : Ice is formed in the open water when ocean lose heat
56      !!      (heat budget of open water Bl is negative) .
57      !!      Computation of the increase of 1-A (ice concentration) fol-
58      !!      lowing the law :
59      !!      (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
60      !!       where - h0 is the thickness of ice created in the lead
61      !!             - a is a minimum fraction for leads
62      !!             - F is a monotonic non-increasing function defined as:
63      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
64      !!             - exld is the exponent closure rate (=2 default val.)
65      !!
66      !! ** Action : - Adjustment of snow and ice thicknesses and heat
67      !!                content in brine pockets
68      !!             - Updating ice internal temperature
69      !!             - Computation of variation of ice volume and mass
70      !!             - Computation of frldb after lateral accretion and
71      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)     
72      !!------------------------------------------------------------------------
73      INTEGER ::   ji,jj,jk,jl      ! dummy loop indices
74      INTEGER ::   nbpac            ! local integers
75      INTEGER ::   ii, ij, iter     !   -       -
76      REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zde                         ! local scalars
77      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
78      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
79      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
80      CHARACTER (len = 15) :: fieldid
81
82      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
83      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
84      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
85      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
86     
87      REAL(wp) ::   zv_newfra
88 
89      INTEGER , POINTER, DIMENSION(:) ::   jcat        ! 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_1d    ! total ice fraction   
101      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
102      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
103
104      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_b      ! old volume of ice in category jl
105      REAL(wp), POINTER, DIMENSION(:,:) ::   za_b      ! old area of ice in category jl
106      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i
107      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i
108      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d  ! 1-D version of oa_i
109      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i
110
111      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i
112
113      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
114      !!-----------------------------------------------------------------------!
115
116      CALL wrk_alloc( jpij, jcat )   ! integer
117      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
118      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
119      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
120      CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d )
121      CALL wrk_alloc( jpi,jpj, zvrel )
122
123      !------------------------------------------------------------------------------|
124      ! 2) Convert units for ice internal energy
125      !------------------------------------------------------------------------------|
126      DO jl = 1, jpl
127         DO jk = 1, nlay_i
128            DO jj = 1, jpj
129               DO ji = 1, jpi
130                  !Energy of melting q(S,T) [J.m-3]
131                  rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
132                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) &
133                      &   / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp )
134                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac
135               END DO
136            END DO
137         END DO
138      END DO
139
140      !------------------------------------------------------------------------------!
141      ! 3) Collection thickness of ice formed in leads and polynyas
142      !------------------------------------------------------------------------------!   
143      ! hicol is the thickness of new ice formed in open water
144      ! hicol can be either prescribed (frazswi = 0)
145      ! or computed (frazswi = 1)
146      ! Frazil ice forms in open water, is transported by wind
147      ! accumulates at the edge of the consolidated ice edge
148      ! where it forms aggregates of a specific thickness called
149      ! collection thickness.
150
151      ! Note : the following algorithm currently breaks vectorization
152      !
153
154      zvrel(:,:) = 0._wp
155
156      ! Default new ice thickness
157      hicol(:,:) = hiccrit
158
159      IF( fraz_swi == 1 ) THEN
160
161         !--------------------
162         ! Physical constants
163         !--------------------
164         hicol(:,:) = 0._wp
165
166         zhicrit = 0.04 ! frazil ice thickness
167         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
168         zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
169         zgamafr = 0.03
170
171         DO jj = 2, jpj
172            DO ji = 2, jpi
173               IF ( qlead(ji,jj) < 0._wp ) THEN
174                  !-------------
175                  ! Wind stress
176                  !-------------
177                  ! C-grid wind stress components
178                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   &
179                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp
180                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   &
181                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp
182                  ! Square root of wind stress
183                  ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) )
184
185                  !---------------------
186                  ! Frazil ice velocity
187                  !---------------------
188                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
189                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
190                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
191
192                  !-------------------
193                  ! Pack ice velocity
194                  !-------------------
195                  ! C-grid ice velocity
196                  rswitch = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
197                  zvgx    = rswitch * ( u_ice(ji-1,jj  ) * tmu(ji-1,jj  )  + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp
198                  zvgy    = rswitch * ( v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)  + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp
199
200                  !-----------------------------------
201                  ! Relative frazil/pack ice velocity
202                  !-----------------------------------
203                  ! absolute relative velocity
204                  zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
205                     &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
206                  zvrel(ji,jj)  = SQRT( zvrel2 )
207
208                  !---------------------
209                  ! Iterative procedure
210                  !---------------------
211                  hicol(ji,jj) = zhicrit + 0.1 
212                  hicol(ji,jj) = zhicrit +   hicol(ji,jj)    &
213                     &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2
214
215!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1)
216!!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01
217!!gm                therefore the 2 lines with hicol can be replaced by 1 line:
218!!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2
219!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop
220
221                  iter = 1
222                  iterate_frazil = .true.
223
224                  DO WHILE ( iter .LT. 100 .AND. iterate_frazil ) 
225                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
226                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
227                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
228                        - zhicrit * ztwogp * zvrel2
229                     zhicol_new = hicol(ji,jj) - zf/zfp
230                     hicol(ji,jj)   = zhicol_new
231
232                     iter = iter + 1
233
234                  END DO ! do while
235
236               ENDIF ! end of selection of pixels where ice forms
237
238            END DO ! loop on ji ends
239         END DO ! loop on jj ends
240      !
241      CALL lbc_lnk( zvrel(:,:), 'T', 1. )
242      CALL lbc_lnk( hicol(:,:), 'T', 1. )
243
244      ENDIF ! End of computation of frazil ice collection thickness
245
246      !------------------------------------------------------------------------------!
247      ! 4) Identify grid points where new ice forms
248      !------------------------------------------------------------------------------!
249
250      !-------------------------------------
251      ! Select points for new ice formation
252      !-------------------------------------
253      ! This occurs if open water energy budget is negative
254      nbpac = 0
255      npac(:) = 0
256      !
257      DO jj = 1, jpj
258         DO ji = 1, jpi
259            IF ( qlead(ji,jj)  <  0._wp ) THEN
260               nbpac = nbpac + 1
261               npac( nbpac ) = (jj - 1) * jpi + ji
262            ENDIF
263         END DO
264      END DO
265
266      ! debug point to follow
267      jiindex_1d = 0
268      IF( ln_nicep ) THEN
269         DO ji = mi0(jiindx), mi1(jiindx)
270            DO jj = mj0(jjindx), mj1(jjindx)
271               IF ( qlead(ji,jj)  <  0._wp ) THEN
272                  jiindex_1d = (jj - 1) * jpi + ji
273               ENDIF
274            END DO
275         END DO
276      ENDIF
277   
278      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
279
280      !------------------------------
281      ! Move from 2-D to 1-D vectors
282      !------------------------------
283      ! If ocean gains heat do nothing
284      ! 0therwise compute new ice formation
285
286      IF ( nbpac > 0 ) THEN
287
288         CALL tab_2d_1d( nbpac, zat_i_1d  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
289         DO jl = 1, jpl
290            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
291            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
292            CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
293            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
294            DO jk = 1, nlay_i
295               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
296            END DO ! jk
297         END DO ! jl
298
299         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) )
300         CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
301         CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw, jpi, jpj, npac(1:nbpac) )
302         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) )
303         CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
304         CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
305
306         CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd, jpi, jpj, npac(1:nbpac) )
307         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw, jpi, jpj, npac(1:nbpac) )
308
309         !------------------------------------------------------------------------------!
310         ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
311         !------------------------------------------------------------------------------!
312
313         !-----------------------------------------
314         ! Keep old ice areas and volume in memory
315         !-----------------------------------------
316         zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
317         za_b(1:nbpac,:) = za_i_1d(1:nbpac,:)
318         !----------------------
319         ! Thickness of new ice
320         !----------------------
321         DO ji = 1, nbpac
322            zh_newice(ji) = hiccrit
323         END DO
324         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
325
326         !----------------------
327         ! Salinity of new ice
328         !----------------------
329         SELECT CASE ( num_sal )
330         CASE ( 1 )                    ! Sice = constant
331            zs_newice(1:nbpac) = bulk_sal
332         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
333            DO ji = 1, nbpac
334               ii =   MOD( npac(ji) - 1 , jpi ) + 1
335               ij =      ( npac(ji) - 1 ) / jpi + 1
336               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  )
337            END DO
338         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
339            zs_newice(1:nbpac) =   2.3
340         END SELECT
341
342         !-------------------------
343         ! Heat content of new ice
344         !-------------------------
345         ! We assume that new ice is formed at the seawater freezing point
346         DO ji = 1, nbpac
347            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
348            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                             &
349               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) )   &
350               &                       - rcp  *         ( ztmelts - rtt )  )
351         END DO ! ji
352
353         !----------------
354         ! Age of new ice
355         !----------------
356         DO ji = 1, nbpac
357            zo_newice(ji) = 0._wp
358         END DO ! ji
359
360         !-------------------
361         ! Volume of new ice
362         !-------------------
363         DO ji = 1, nbpac
364
365            zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg]
366
367            zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg]
368                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
369                                                                   
370            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
371                                             
372            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)
373                                                                   ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point   
374            zv_newice(ji) = - zfmdt / rhoic
375
376            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux 
377
378            ! Contribution to heat flux to the ocean [W.m-2], >0 
379            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
380            ! Total heat flux used in this process [W.m-2] 
381            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
382            ! mass flux
383            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice
384            ! salt flux
385            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
386
387            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
388            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
389            zfrazb        = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
390            zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
391            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
392         END DO
393
394         !-----------------
395         ! Area of new ice
396         !-----------------
397         DO ji = 1, nbpac
398            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
399         END DO
400
401         !------------------------------------------------------------------------------!
402         ! 6) Redistribute new ice area and volume into ice categories                  !
403         !------------------------------------------------------------------------------!
404
405         !------------------------
406         ! 6.1) lateral ice growth
407         !------------------------
408         ! If lateral ice growth gives an ice concentration gt 1, then
409         ! we keep the excessive volume in memory and attribute it later to bottom accretion
410         DO ji = 1, nbpac
411            IF ( za_newice(ji) >  ( amax - zat_i_1d(ji) ) ) THEN
412               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_1d(ji) )
413               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
414               za_newice(ji) = za_newice(ji) - zda_res  (ji)
415               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
416            ELSE
417               zda_res(ji) = 0._wp
418               zdv_res(ji) = 0._wp
419            ENDIF
420         END DO
421
422         ! find which category to fill
423         zat_i_1d(:) = 0._wp
424         DO jl = 1, jpl
425            DO ji = 1, nbpac
426               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
427                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
428                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
429                  jcat    (ji)    = jl
430               ENDIF
431               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
432            END DO
433         END DO
434
435         ! Heat content
436         DO ji = 1, nbpac
437            jl = jcat(ji)                                                    ! categroy in which new ice is put
438            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
439         END DO
440
441         DO jk = 1, nlay_i
442            DO ji = 1, nbpac
443               jl = jcat(ji)
444               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
445               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      &
446                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  &
447                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
448            END DO
449         END DO
450
451         !------------------------------------------------
452         ! 6.2) bottom ice growth + ice enthalpy remapping
453         !------------------------------------------------
454         DO jl = 1, jpl
455
456            ! for remapping
457            h_i_old (1:nbpac,0:nlay_i+1) = 0._wp
458            qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp
459            DO jk = 1, nlay_i
460               DO ji = 1, nbpac
461                  h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )
462                  qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
463               END DO
464            END DO
465
466            ! new volumes including lateral/bottom accretion + residual
467            DO ji = 1, nbpac
468               rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
469               zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
470               za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)               
471               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
472               ! for remapping
473               h_i_old (ji,nlay_i+1) = zv_newfra
474               qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
475            ENDDO
476            ! --- Ice enthalpy remapping --- !
477            CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
478         ENDDO
479
480         !------------
481         ! Update age
482         !------------
483         DO jl = 1, jpl
484            DO ji = 1, nbpac
485               rswitch          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes
486               zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch   
487            END DO
488         END DO   
489
490         !-----------------
491         ! Update salinity
492         !-----------------
493         DO jl = 1, jpl
494            DO ji = 1, nbpac
495               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
496               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
497            END DO
498         END DO
499
500         !------------------------------------------------------------------------------!
501         ! 7) Change 2D vectors to 1D vectors
502         !------------------------------------------------------------------------------!
503         DO jl = 1, jpl
504            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj )
505            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj )
506            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj )
507            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj )
508            DO jk = 1, nlay_i
509               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj )
510            END DO
511         END DO
512         CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj )
513         CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
514
515         CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj )
516         CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj )
517         !
518      ENDIF ! nbpac > 0
519
520      !------------------------------------------------------------------------------!
521      ! 8) Change units for e_i
522      !------------------------------------------------------------------------------!   
523      DO jl = 1, jpl
524         DO jk = 1, nlay_i
525            DO jj = 1, jpj
526               DO ji = 1, jpi
527                  ! heat content in Joules
528                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac ) 
529               END DO
530            END DO
531         END DO
532      END DO
533
534      !
535      CALL wrk_dealloc( jpij, jcat )   ! integer
536      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
537      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
538      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
539      CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d )
540      CALL wrk_dealloc( jpi,jpj, zvrel )
541      !
542   END SUBROUTINE lim_thd_lac
543
544#else
545   !!----------------------------------------------------------------------
546   !!   Default option                               NO  LIM3 sea-ice model
547   !!----------------------------------------------------------------------
548CONTAINS
549   SUBROUTINE lim_thd_lac           ! Empty routine
550   END SUBROUTINE lim_thd_lac
551#endif
552
553   !!======================================================================
554END MODULE limthd_lac
Note: See TracBrowser for help on using the repository browser.