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

source: branches/2014/dev_r5134_UKMO4_CF_compliance/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 5136

Last change on this file since 5136 was 5134, checked in by clem, 9 years ago

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