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/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icethd_do.F90 @ 15374

Last change on this file since 15374 was 15334, checked in by clem, 3 years ago

Some harmless reorganization of SI3: 1) extract the parts where mpi com were needed from inside thermo. 2) code an optional upstream scheme inside rheology to calculate P as the sub time step level. 3) prepare the albedo to scheme to recieve an aditional namelist parameter (pivotal ice thickness)
  • Property svn:keywords set to Id
File size: 18.9 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 ice1D          ! sea-ice: thermodynamics variables
20   USE ice            ! sea-ice: variables
21   USE icetab         ! sea-ice: 2D <==> 1D
22   USE icectl         ! sea-ice: conservation
23   USE icethd_ent     ! sea-ice: thermodynamics, enthalpy
24   USE icevar         ! sea-ice: operations
25   USE icethd_sal     ! sea-ice: salinity profiles
26   !
27   USE in_out_manager ! I/O manager
28   USE lib_mpp        ! MPP library
29   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
30   USE lbclnk         ! lateral boundary conditions (or mpp links)
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   ice_thd_do        ! called by ice_thd
36   PUBLIC   ice_thd_do_init   ! called by ice_stp
37
38   !                          !!** namelist (namthd_do) **
39   REAL(wp), PUBLIC ::   rn_hinew      ! thickness for new ice formation (m)
40   LOGICAL , PUBLIC ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F)
41   REAL(wp), PUBLIC ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom
42   REAL(wp), PUBLIC ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice
43   REAL(wp), PUBLIC ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice
44
45   !! * Substitutions
46#  include "do_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE ice_thd_do
55      !!-------------------------------------------------------------------
56      !!               ***   ROUTINE ice_thd_do  ***
57      !! 
58      !! ** Purpose : Computation of the evolution of the ice thickness and
59      !!              concentration as a function of the heat balance in the leads
60      !!       
61      !! ** Method  : Ice is formed in the open water when ocean looses heat
62      !!              (heat budget of open water is negative) following
63      !!
64      !!       (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
65      !!          where - h0 is the thickness of ice created in the lead
66      !!                - a is a minimum fraction for leads
67      !!                - F is a monotonic non-increasing function defined as:
68      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
69      !!                - exld is the exponent closure rate (=2 default val.)
70      !!
71      !! ** Action : - Adjustment of snow and ice thicknesses and heat
72      !!                content in brine pockets
73      !!             - Updating ice internal temperature
74      !!             - Computation of variation of ice volume and mass
75      !!             - Computation of a_i after lateral accretion and
76      !!               update h_s_1d, h_i_1d     
77      !!------------------------------------------------------------------------
78      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
79      !
80      REAL(wp) ::   ztmelts
81      REAL(wp) ::   zdE
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 , DIMENSION(jpij) ::   jcat        ! indexes of categories where new ice grows
90      REAL(wp), DIMENSION(jpij) ::   zswinew     ! switch for new ice or not
91      !
92      REAL(wp), DIMENSION(jpij) ::   zv_newice   ! volume of accreted ice
93      REAL(wp), DIMENSION(jpij) ::   za_newice   ! fractional area of accreted ice
94      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! thickness of accreted ice
95      REAL(wp), DIMENSION(jpij) ::   ze_newice   ! heat content of accreted ice
96      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of accreted ice
97      REAL(wp), DIMENSION(jpij) ::   zo_newice   ! age of accreted ice
98      REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget
99      REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
100      REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
101      REAL(wp), DIMENSION(jpij) ::   zfraz_frac_1d ! relative ice / frazil velocity (1D vector)
102      !
103      REAL(wp), DIMENSION(jpij,jpl) ::   zv_b    ! old volume of ice in category jl
104      REAL(wp), DIMENSION(jpij,jpl) ::   za_b    ! old area of ice in category jl
105      !
106      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i
107      !
108      !!-----------------------------------------------------------------------!
109
110      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
111      IF( ln_icediachk )   CALL ice_cons2D  ( 0, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft )
112
113      at_i(:,:) = SUM( a_i, dim=3 )
114      !------------------------------------------------------------------------------!
115      ! 1) Compute thickness, salinity, enthalpy, age, area and volume of new ice
116      !------------------------------------------------------------------------------!
117      ! it occurs if cooling
118
119      ! Identify grid points where new ice forms
120      npti = 0   ;   nptidx(:) = 0
121      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
122         IF ( qlead(ji,jj)  <  0._wp ) THEN
123            npti = npti + 1
124            nptidx( npti ) = (jj - 1) * jpi + ji
125         ENDIF
126      END_2D
127
128      ! Move from 2-D to 1-D vectors
129      IF ( npti > 0 ) THEN
130
131         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)      , at_i        )
132         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
133         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
134         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
135         DO jl = 1, jpl
136            DO jk = 1, nlay_i
137               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
138            END DO
139         END DO
140         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d     (1:npti) , qlead      )
141         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d      (1:npti) , t_bo       )
142         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d   (1:npti) , sfx_opw    )
143         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d   (1:npti) , wfx_opw    )
144         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice    (1:npti) , ht_i_new   )
145         CALL tab_2d_1d( npti, nptidx(1:npti), zfraz_frac_1d(1:npti) , fraz_frac  )
146
147         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd    )
148         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw    )
149         CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
150         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m      )
151
152         ! Convert units for ice internal energy
153         DO jl = 1, jpl
154            DO jk = 1, nlay_i               
155               WHERE( v_i_2d(1:npti,jl) > 0._wp )
156                  ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
157               ELSEWHERE
158                  ze_i_2d(1:npti,jk,jl) = 0._wp
159               END WHERE
160            END DO
161         END DO
162
163         ! Keep old ice areas and volume in memory
164         zv_b(1:npti,:) = v_i_2d(1:npti,:) 
165         za_b(1:npti,:) = a_i_2d(1:npti,:)
166
167         ! --- Salinity of new ice --- !
168         SELECT CASE ( nn_icesal )
169         CASE ( 1 )                    ! Sice = constant
170            zs_newice(1:npti) = rn_icesal
171         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
172            DO ji = 1, npti
173               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
174            END DO
175         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
176            zs_newice(1:npti) =   2.3
177         END SELECT
178
179         ! --- Heat content of new ice --- !
180         ! We assume that new ice is formed at the seawater freezing point
181         DO ji = 1, npti
182            ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C)
183            ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
184               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
185               &                      - rcp   *         ztmelts )
186         END DO
187
188         ! --- Age of new ice --- !
189         zo_newice(1:npti) = 0._wp
190
191         ! --- Volume of new ice --- !
192         DO ji = 1, npti
193
194            zEi           = - ze_newice(ji) * r1_rhoi              ! specific enthalpy of forming ice [J/kg]
195
196            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
197                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
198                                                                   
199            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
200                                             
201            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)
202                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point   
203            zv_newice(ji) = - zfmdt * r1_rhoi
204
205            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux 
206
207            ! Contribution to heat flux to the ocean [W.m-2], >0 
208            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_Dt_ice
209            ! Total heat flux used in this process [W.m-2] 
210            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_Dt_ice
211            ! mass flux
212            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_Dt_ice
213            ! salt flux
214            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_Dt_ice
215         END DO
216         
217         ! A fraction fraz_frac of frazil ice is accreted at the ice bottom
218         DO ji = 1, npti
219            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
220            zv_frazb(ji)  =           zfraz_frac_1d(ji) * rswitch   * zv_newice(ji)
221            zv_newice(ji) = ( 1._wp - zfraz_frac_1d(ji) * rswitch ) * zv_newice(ji)
222         END DO
223         
224         ! --- Area of new ice --- !
225         DO ji = 1, npti
226            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
227         END DO
228
229         !------------------------------------------------------------------------------!
230         ! 2) Redistribute new ice area and volume into ice categories                  !
231         !------------------------------------------------------------------------------!
232
233         ! --- lateral ice growth --- !
234         ! If lateral ice growth gives an ice concentration > amax, then
235         ! we keep the excessive volume in memory and attribute it later to bottom accretion
236         DO ji = 1, npti
237            IF ( za_newice(ji) >  MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
238               zda_res(ji)   = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
239               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
240               za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res  (ji) )
241               zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res  (ji) )
242            ELSE
243               zda_res(ji) = 0._wp
244               zdv_res(ji) = 0._wp
245            ENDIF
246         END DO
247
248         ! find which category to fill
249         DO jl = 1, jpl
250            DO ji = 1, npti
251               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
252                  a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
253                  v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
254                  jcat(ji) = jl
255               ENDIF
256            END DO
257         END DO
258         at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
259
260         ! Heat content
261         DO ji = 1, npti
262            jl = jcat(ji)                                                    ! categroy in which new ice is put
263            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
264         END DO
265
266         DO jk = 1, nlay_i
267            DO ji = 1, npti
268               jl = jcat(ji)
269               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
270               ze_i_2d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
271                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) )  &
272                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
273            END DO
274         END DO
275
276         ! --- bottom ice growth + ice enthalpy remapping --- !
277         DO jl = 1, jpl
278
279            ! for remapping
280            h_i_old (1:npti,0:nlay_i+1) = 0._wp
281            eh_i_old(1:npti,0:nlay_i+1) = 0._wp
282            DO jk = 1, nlay_i
283               DO ji = 1, npti
284                  h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
285                  eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
286               END DO
287            END DO
288
289            ! new volumes including lateral/bottom accretion + residual
290            DO ji = 1, npti
291               rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
292               zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
293               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)               
294               v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
295               ! for remapping
296               h_i_old (ji,nlay_i+1) = zv_newfra
297               eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
298            END DO
299            ! --- Ice enthalpy remapping --- !
300            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 
301         END DO
302
303         ! --- Update salinity --- !
304         DO jl = 1, jpl
305            DO ji = 1, npti
306               sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
307            END DO
308         END DO
309
310         ! Change units for e_i
311         DO jl = 1, jpl
312            DO jk = 1, nlay_i
313               ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i 
314            END DO
315         END DO
316
317         ! Move 2D vectors to 1D vectors
318         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
319         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
320         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
321          DO jl = 1, jpl
322            DO jk = 1, nlay_i
323               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
324            END DO
325         END DO
326         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
327         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
328         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
329         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw )
330         !
331      ENDIF ! npti > 0
332      !
333      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
334      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
335      !
336   END SUBROUTINE ice_thd_do
337
338
339   SUBROUTINE ice_thd_do_init
340      !!-----------------------------------------------------------------------
341      !!                   ***  ROUTINE ice_thd_do_init ***
342      !!                 
343      !! ** Purpose :   Physical constants and parameters associated with
344      !!                ice growth in the leads
345      !!
346      !! ** Method  :   Read the namthd_do namelist and check the parameters
347      !!                called at the first timestep (nit000)
348      !!
349      !! ** input   :   Namelist namthd_do
350      !!-------------------------------------------------------------------
351      INTEGER  ::   ios   ! Local integer
352      !!
353      NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
354      !!-------------------------------------------------------------------
355      !
356      READ  ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901)
357901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namthd_do in reference namelist' )
358      READ  ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 )
359902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namthd_do in configuration namelist' )
360      IF(lwm) WRITE( numoni, namthd_do )
361      !
362      IF(lwp) THEN                          ! control print
363         WRITE(numout,*)
364         WRITE(numout,*) 'ice_thd_do_init: Ice growth in open water'
365         WRITE(numout,*) '~~~~~~~~~~~~~~~'
366         WRITE(numout,*) '   Namelist namthd_do:'
367         WRITE(numout,*) '      ice thickness for lateral accretion                       rn_hinew   = ', rn_hinew
368         WRITE(numout,*) '      Frazil ice thickness as a function of wind or not         ln_frazil  = ', ln_frazil
369         WRITE(numout,*) '      Maximum proportion of frazil ice collecting at bottom     rn_maxfraz = ', rn_maxfraz
370         WRITE(numout,*) '      Threshold relative drift speed for collection of frazil   rn_vfraz   = ', rn_vfraz
371         WRITE(numout,*) '      Squeezing coefficient for collection of frazil            rn_Cfraz   = ', rn_Cfraz
372      ENDIF
373      !
374      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_do_init : rn_hinew should be >= rn_himin' )
375      !
376   END SUBROUTINE ice_thd_do_init
377   
378#else
379   !!----------------------------------------------------------------------
380   !!   Default option                                NO SI3 sea-ice model
381   !!----------------------------------------------------------------------
382#endif
383
384   !!======================================================================
385END MODULE icethd_do
Note: See TracBrowser for help on using the repository browser.