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.
icethd_do.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icethd_do.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 23.5 KB
Line 
1MODULE icethd_do
2   !!======================================================================
3   !!                       ***  MODULE icethd_do   ***
4   !!   sea-ice: sea ice growth in the leads (open water) 
5   !!======================================================================
6   !! History :       !  2005-12  (M. Vancoppenolle) Original code
7   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
8   !!----------------------------------------------------------------------
9#if defined key_si3
10   !!----------------------------------------------------------------------
11   !!   'key_si3'                                       SI3 sea-ice model
12   !!----------------------------------------------------------------------
13   !!   ice_thd_do        : ice growth in open water (=lateral accretion of ice)
14   !!   ice_thd_do_init   : initialization
15   !!----------------------------------------------------------------------
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE sbc_oce , ONLY : sss_m
19   USE sbc_ice , ONLY : utau_ice, vtau_ice
20   USE ice1D          ! sea-ice: thermodynamics variables
21   USE ice            ! sea-ice: variables
22   USE icetab         ! sea-ice: 2D <==> 1D
23   USE icectl         ! sea-ice: conservation
24   USE icethd_ent     ! sea-ice: thermodynamics, enthalpy
25   USE icevar         ! sea-ice: operations
26   USE icethd_sal     ! sea-ice: salinity profiles
27   !
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
31   USE lbclnk         ! lateral boundary conditions (or mpp links)
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ice_thd_do        ! called by ice_thd
37   PUBLIC   ice_thd_do_init   ! called by ice_stp
38
39   !                          !!** namelist (namthd_do) **
40   REAL(wp) ::   rn_hinew      ! thickness for new ice formation (m)
41   LOGICAL  ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F)
42   REAL(wp) ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom
43   REAL(wp) ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice
44   REAL(wp) ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice
45
46   !! * Substitutions
47#  include "do_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE ice_thd_do
56      !!-------------------------------------------------------------------
57      !!               ***   ROUTINE ice_thd_do  ***
58      !! 
59      !! ** Purpose : Computation of the evolution of the ice thickness and
60      !!              concentration as a function of the heat balance in the leads
61      !!       
62      !! ** Method  : Ice is formed in the open water when ocean looses heat
63      !!              (heat budget of open water is negative) following
64      !!
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 a_i after lateral accretion and
77      !!               update h_s_1d, h_i_1d     
78      !!------------------------------------------------------------------------
79      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
80      INTEGER  ::   iter             !   -       -
81      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                               ! local scalars
82      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
83      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
84      !
85      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
86      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
87      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
88      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
89      !
90      REAL(wp) ::   zv_newfra
91      !
92      INTEGER , DIMENSION(jpij) ::   jcat        ! indexes of categories where new ice grows
93      REAL(wp), DIMENSION(jpij) ::   zswinew     ! switch for new ice or not
94      !
95      REAL(wp), DIMENSION(jpij) ::   zv_newice   ! volume of accreted ice
96      REAL(wp), DIMENSION(jpij) ::   za_newice   ! fractional area of accreted ice
97      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! thickness of accreted ice
98      REAL(wp), DIMENSION(jpij) ::   ze_newice   ! heat content of accreted ice
99      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of accreted ice
100      REAL(wp), DIMENSION(jpij) ::   zo_newice   ! age of accreted ice
101      REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget
102      REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
103      REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
104      REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
105      !
106      REAL(wp), DIMENSION(jpij,jpl) ::   zv_b    ! old volume of ice in category jl
107      REAL(wp), DIMENSION(jpij,jpl) ::   za_b    ! old area of ice in category jl
108      !
109      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i
110      !
111      REAL(wp), DIMENSION(jpi,jpj) ::   zvrel    ! relative ice / frazil velocity
112      !
113      REAL(wp) :: zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used)
114      !!-----------------------------------------------------------------------!
115
116      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
117      IF( ln_icediachk )   CALL ice_cons2D  ( 0, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft )
118
119      at_i(:,:) = SUM( a_i, dim=3 )
120      !------------------------------------------------------------------------------!
121      ! 1) Collection thickness of ice formed in leads and polynyas
122      !------------------------------------------------------------------------------!   
123      ! ht_i_new is the thickness of new ice formed in open water
124      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
125      ! Frazil ice forms in open water, is transported by wind
126      ! accumulates at the edge of the consolidated ice edge
127      ! where it forms aggregates of a specific thickness called
128      ! collection thickness.
129
130      zvrel(:,:) = 0._wp
131
132      ! Default new ice thickness
133      WHERE( qlead(:,:) < 0._wp  .AND. tau_icebfr(:,:) == 0._wp )   ;   ht_i_new(:,:) = rn_hinew ! if cooling and no landfast
134      ELSEWHERE                                                     ;   ht_i_new(:,:) = 0._wp
135      END WHERE
136
137      IF( ln_frazil ) THEN
138         !
139         ht_i_new(:,:) = 0._wp
140         !
141         ! Physical constants
142         zhicrit = 0.04                                          ! frazil ice thickness
143         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoi ) )  ! reduced grav
144         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag)
145         zgamafr = 0.03
146         !
147         DO_2D_00_00
148            IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
149               ! -- Wind stress -- !
150               ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
151                  &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
152               ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
153                  &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
154               ! Square root of wind stress
155               ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
156
157               ! -- Frazil ice velocity -- !
158               rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
159               zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
160               zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
161
162               ! -- Pack ice velocity -- !
163               zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
164               zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
165
166               ! -- Relative frazil/pack ice velocity -- !
167               rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
168               zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
169                  &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
170               zvrel(ji,jj) = SQRT( zvrel2 )
171
172               ! -- new ice thickness (iterative loop) -- !
173               ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
174                  &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
175
176               iter = 1
177               DO WHILE ( iter < 20 ) 
178                  zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   &
179                     &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
180                  zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
181
182                  ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
183                  iter = iter + 1
184               END DO
185               !
186               ! bound ht_i_new (though I don't see why it should be necessary)
187               ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
188               !
189            ENDIF
190            !
191         END_2D
192         !
193         CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1., ht_i_new, 'T', 1.  )
194
195      ENDIF
196
197      !------------------------------------------------------------------------------!
198      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice
199      !------------------------------------------------------------------------------!
200      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
201
202      ! Identify grid points where new ice forms
203      npti = 0   ;   nptidx(:) = 0
204      DO_2D_11_11
205         IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
206            npti = npti + 1
207            nptidx( npti ) = (jj - 1) * jpi + ji
208         ENDIF
209      END_2D
210
211      ! Move from 2-D to 1-D vectors
212      IF ( npti > 0 ) THEN
213
214         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)      , at_i        )
215         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
216         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
217         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
218         DO jl = 1, jpl
219            DO jk = 1, nlay_i
220               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
221            END DO
222         END DO
223         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead      )
224         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo       )
225         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw    )
226         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw    )
227         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new   )
228         CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel      )
229
230         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd    )
231         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw    )
232         CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
233         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m      )
234
235         ! Convert units for ice internal energy
236         DO jl = 1, jpl
237            DO jk = 1, nlay_i               
238               WHERE( v_i_2d(1:npti,jl) > 0._wp )
239                  ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
240               ELSEWHERE
241                  ze_i_2d(1:npti,jk,jl) = 0._wp
242               END WHERE
243            END DO
244         END DO
245
246         ! Keep old ice areas and volume in memory
247         zv_b(1:npti,:) = v_i_2d(1:npti,:) 
248         za_b(1:npti,:) = a_i_2d(1:npti,:)
249
250         ! --- Salinity of new ice --- !
251         SELECT CASE ( nn_icesal )
252         CASE ( 1 )                    ! Sice = constant
253            zs_newice(1:npti) = rn_icesal
254         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
255            DO ji = 1, npti
256               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
257            END DO
258         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
259            zs_newice(1:npti) =   2.3
260         END SELECT
261
262         ! --- Heat content of new ice --- !
263         ! We assume that new ice is formed at the seawater freezing point
264         DO ji = 1, npti
265            ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
266            ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
267               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
268               &                      - rcp   *         ztmelts )
269         END DO
270
271         ! --- Age of new ice --- !
272         zo_newice(1:npti) = 0._wp
273
274         ! --- Volume of new ice --- !
275         DO ji = 1, npti
276
277            zEi           = - ze_newice(ji) * r1_rhoi              ! specific enthalpy of forming ice [J/kg]
278
279            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
280                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
281                                                                   
282            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
283                                             
284            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)
285                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point   
286            zv_newice(ji) = - zfmdt * r1_rhoi
287
288            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux 
289
290            ! Contribution to heat flux to the ocean [W.m-2], >0 
291            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
292            ! Total heat flux used in this process [W.m-2] 
293            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
294            ! mass flux
295            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_rdtice
296            ! salt flux
297            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_rdtice
298         END DO
299         
300         zv_frazb(1:npti) = 0._wp
301         IF( ln_frazil ) THEN
302            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
303            DO ji = 1, npti
304               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
305               zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
306               zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
307               zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
308            END DO
309         END IF
310         
311         ! --- Area of new ice --- !
312         DO ji = 1, npti
313            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
314         END DO
315
316         !------------------------------------------------------------------------------!
317         ! 3) Redistribute new ice area and volume into ice categories                  !
318         !------------------------------------------------------------------------------!
319
320         ! --- lateral ice growth --- !
321         ! If lateral ice growth gives an ice concentration > amax, then
322         ! we keep the excessive volume in memory and attribute it later to bottom accretion
323         DO ji = 1, npti
324            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
325               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
326               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
327               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) )
328               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) )
329            ELSE
330               zda_res(ji) = 0._wp
331               zdv_res(ji) = 0._wp
332            ENDIF
333         END DO
334
335         ! find which category to fill
336         DO jl = 1, jpl
337            DO ji = 1, npti
338               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
339                  a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
340                  v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
341                  jcat(ji) = jl
342               ENDIF
343            END DO
344         END DO
345         at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
346
347         ! Heat content
348         DO ji = 1, npti
349            jl = jcat(ji)                                                    ! categroy in which new ice is put
350            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
351         END DO
352
353         DO jk = 1, nlay_i
354            DO ji = 1, npti
355               jl = jcat(ji)
356               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
357               ze_i_2d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
358                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) )  &
359                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
360            END DO
361         END DO
362
363         ! --- bottom ice growth + ice enthalpy remapping --- !
364         DO jl = 1, jpl
365
366            ! for remapping
367            h_i_old (1:npti,0:nlay_i+1) = 0._wp
368            eh_i_old(1:npti,0:nlay_i+1) = 0._wp
369            DO jk = 1, nlay_i
370               DO ji = 1, npti
371                  h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
372                  eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
373               END DO
374            END DO
375
376            ! new volumes including lateral/bottom accretion + residual
377            DO ji = 1, npti
378               rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
379               zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
380               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)               
381               v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
382               ! for remapping
383               h_i_old (ji,nlay_i+1) = zv_newfra
384               eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
385            END DO
386            ! --- Ice enthalpy remapping --- !
387            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 
388         END DO
389
390         ! --- Update salinity --- !
391         DO jl = 1, jpl
392            DO ji = 1, npti
393               sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
394            END DO
395         END DO
396
397         ! Change units for e_i
398         DO jl = 1, jpl
399            DO jk = 1, nlay_i
400               ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i 
401            END DO
402         END DO
403
404         ! Move 2D vectors to 1D vectors
405         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
406         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
407         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
408          DO jl = 1, jpl
409            DO jk = 1, nlay_i
410               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
411            END DO
412         END DO
413         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
414         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
415         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
416         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw )
417         !
418      ENDIF ! npti > 0
419      !
420      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
421      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
422      !
423   END SUBROUTINE ice_thd_do
424
425
426   SUBROUTINE ice_thd_do_init
427      !!-----------------------------------------------------------------------
428      !!                   ***  ROUTINE ice_thd_do_init ***
429      !!                 
430      !! ** Purpose :   Physical constants and parameters associated with
431      !!                ice growth in the leads
432      !!
433      !! ** Method  :   Read the namthd_do namelist and check the parameters
434      !!                called at the first timestep (nit000)
435      !!
436      !! ** input   :   Namelist namthd_do
437      !!-------------------------------------------------------------------
438      INTEGER  ::   ios   ! Local integer
439      !!
440      NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
441      !!-------------------------------------------------------------------
442      !
443      READ  ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901)
444901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_do in reference namelist' )
445      READ  ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 )
446902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_do in configuration namelist' )
447      IF(lwm) WRITE( numoni, namthd_do )
448      !
449      IF(lwp) THEN                          ! control print
450         WRITE(numout,*)
451         WRITE(numout,*) 'ice_thd_do_init: Ice growth in open water'
452         WRITE(numout,*) '~~~~~~~~~~~~~~~'
453         WRITE(numout,*) '   Namelist namthd_do:'
454         WRITE(numout,*) '      ice thickness for lateral accretion                       rn_hinew   = ', rn_hinew
455         WRITE(numout,*) '      Frazil ice thickness as a function of wind or not         ln_frazil  = ', ln_frazil
456         WRITE(numout,*) '      Maximum proportion of frazil ice collecting at bottom     rn_maxfraz = ', rn_maxfraz
457         WRITE(numout,*) '      Threshold relative drift speed for collection of frazil   rn_vfraz   = ', rn_vfraz
458         WRITE(numout,*) '      Squeezing coefficient for collection of frazil            rn_Cfraz   = ', rn_Cfraz
459      ENDIF
460      !
461      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_do_init : rn_hinew should be >= rn_himin' )
462      !
463   END SUBROUTINE ice_thd_do_init
464   
465#else
466   !!----------------------------------------------------------------------
467   !!   Default option                                NO SI3 sea-ice model
468   !!----------------------------------------------------------------------
469#endif
470
471   !!======================================================================
472END MODULE icethd_do
Note: See TracBrowser for help on using the repository browser.