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.
icedyn.F90 in NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE – NEMO

source: NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90 @ 10312

Last change on this file since 10312 was 10312, checked in by clem, 5 years ago

add the parameterization of landfast ice from Lemieux2015 and 2016 with the addition of isotropic tensile strength
File size: 15.7 KB
Line 
1MODULE icedyn
2   !!======================================================================
3   !!                     ***  MODULE  icedyn  ***
4   !!   Sea-Ice dynamics : master routine for sea ice dynamics
5   !!======================================================================
6   !! history :  4.0  ! 2018  (C. Rousset)  original code SI3 [aka Sea Ice cube]
7   !!----------------------------------------------------------------------
8#if defined key_si3
9   !!----------------------------------------------------------------------
10   !!   'key_si3'                                       SI3 sea-ice model
11   !!----------------------------------------------------------------------
12   !!   ice_dyn       : dynamics of sea ice
13   !!   ice_dyn_init  : initialization and namelist read
14   !!----------------------------------------------------------------------
15   USE phycst         ! physical constants
16   USE dom_oce        ! ocean space and time domain
17   USE ice            ! sea-ice: variables
18   USE icedyn_rhg     ! sea-ice: rheology
19   USE icedyn_adv     ! sea-ice: advection
20   USE icedyn_rdgrft  ! sea-ice: ridging/rafting
21   USE icecor         ! sea-ice: corrections
22   USE icevar         ! sea-ice: operations
23   !
24   USE in_out_manager ! I/O manager
25   USE iom            ! I/O manager library
26   USE lib_mpp        ! MPP library
27   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
28   USE lbclnk         ! lateral boundary conditions (or mpp links)
29   USE timing         ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   ice_dyn        ! called by icestp.F90
35   PUBLIC   ice_dyn_init   ! called by icestp.F90
36   
37   INTEGER ::              nice_dyn   ! choice of the type of dynamics
38   !                                        ! associated indices:
39   INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
40   INTEGER, PARAMETER ::   np_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection)
41   INTEGER, PARAMETER ::   np_dynADV     = 3   ! only advection w prescribed vel.(rn_uvice + advection)
42   !
43   ! ** namelist (namdyn) **
44   LOGICAL  ::   ln_dynFULL       ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
45   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections  (rheology + advection)
46   LOGICAL  ::   ln_dynADV        ! only advection w prescribed vel.(rn_uvice + advection)
47   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV)
48   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV)
49   
50   !! * Substitutions
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
54   !! $Id: icedyn.F90 8378 2017-07-26 13:55:59Z clem $
55   !! Software governed by the CeCILL licence     (./LICENSE)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE ice_dyn( kt )
60      !!-------------------------------------------------------------------
61      !!               ***  ROUTINE ice_dyn  ***
62      !!               
63      !! ** Purpose :   this routine manages sea ice dynamics
64      !!
65      !! ** Action : - calculation of friction in case of landfast ice
66      !!             - call ice_dyn_rhg    = rheology
67      !!             - call ice_dyn_adv    = advection
68      !!             - call ice_dyn_rdgrft = ridging/rafting
69      !!             - call ice_cor        = corrections if fields are out of bounds
70      !!--------------------------------------------------------------------
71      INTEGER, INTENT(in) ::   kt     ! ice time step
72      !!
73      INTEGER  ::   ji, jj, jl        ! dummy loop indices
74      REAL(wp) ::   zcoefu, zcoefv
75      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhmax
76      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
77      !!--------------------------------------------------------------------
78      !
79      IF( ln_timing )   CALL timing_start('icedyn')
80      !
81      IF( kt == nit000 .AND. lwp ) THEN
82         WRITE(numout,*)
83         WRITE(numout,*)'ice_dyn: sea-ice dynamics'
84         WRITE(numout,*)'~~~~~~~'
85      ENDIF
86      !                     
87      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
88         tau_icebfr(:,:) = 0._wp
89         DO jl = 1, jpl
90            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
91         END DO
92      ENDIF
93      !
94      zhmax(:,:,:) = h_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
95      DO jl = 1, jpl
96         DO jj = 2, jpjm1
97            DO ji = 2, jpim1
98!!gm use of MAXVAL here is very probably less efficient than expending the 9 values
99               zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( h_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) )
100            END DO
101         END DO
102      END DO
103      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. )
104      !
105      !
106      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
107
108      CASE ( np_dynFULL )          !==  all dynamical processes  ==!
109         CALL ice_dyn_rhg   ( kt )                            ! -- rheology 
110         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness
111         CALL ice_dyn_rdgrft( kt )                            ! -- ridging/rafting
112         CALL ice_cor       ( kt , 1 )                        ! -- Corrections
113
114      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
115         CALL ice_dyn_rhg   ( kt )                            ! -- rheology 
116         CALL ice_dyn_adv   ( kt )                            ! -- advection of ice
117         CALL Hpiling                                         ! -- simple pile-up (replaces ridging/rafting)
118
119      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities)
120         ALLOCATE( zdivu_i(jpi,jpj) )
121
122         u_ice(:,:) = rn_uice * umask(:,:,1)
123         v_ice(:,:) = rn_vice * vmask(:,:,1)
124         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
125         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
126         ! --- monotonicity test from Lipscomb et al 2004 --- !
127         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
128         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s
129         !DO jj = 1, jpj
130         !   DO ji = 1, jpi
131         !      zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
132         !      zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
133         !      u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
134         !      v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
135         !   END DO
136         !END DO
137         ! ---
138         CALL ice_dyn_adv   ( kt )                            ! -- advection of ice
139
140         ! diagnostics: divergence at T points
141         DO jj = 2, jpjm1
142            DO ji = 2, jpim1
143               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
144                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
145            END DO
146         END DO
147         CALL lbc_lnk( zdivu_i, 'T', 1. )
148         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
149
150         DEALLOCATE( zdivu_i )
151
152      END SELECT
153      !
154      IF( ln_timing )   CALL timing_stop('icedyn')
155      !
156   END SUBROUTINE ice_dyn
157
158
159   SUBROUTINE Hbig( phmax )
160      !!-------------------------------------------------------------------
161      !!                  ***  ROUTINE Hbig  ***
162      !!
163      !! ** Purpose : Thickness correction in case advection scheme creates
164      !!              abnormally tick ice
165      !!
166      !! ** Method  : 1- check whether ice thickness resulting from advection is
167      !!                 larger than the surrounding 9-points before advection
168      !!                 and reduce it if a) divergence or b) convergence & at_i>0.8
169      !!              2- bound ice thickness with hi_max (99m)
170      !!
171      !! ** input   : Max thickness of the surrounding 9-points
172      !!-------------------------------------------------------------------
173      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts
174      !
175      INTEGER  ::   ji, jj, jl         ! dummy loop indices
176      REAL(wp) ::   zh, zdv
177      !!-------------------------------------------------------------------
178      !
179      CALL ice_var_zapsmall                       !-- zap small areas
180      !
181      DO jl = 1, jpl
182         DO jj = 1, jpj
183            DO ji = 1, jpi
184               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax
185                  !
186                  zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl)
187                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl) 
188                  !
189                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
190                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
191                     a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
192                  ENDIF
193                  !
194               ENDIF
195            END DO
196         END DO
197      END DO 
198      !                                           !-- correct pond fraction to avoid a_ip > a_i
199      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
200      !
201   END SUBROUTINE Hbig
202
203
204   SUBROUTINE Hpiling
205      !!-------------------------------------------------------------------
206      !!                  ***  ROUTINE Hpiling  ***
207      !!
208      !! ** Purpose : Simple conservative piling comparable with 1-cat models
209      !!
210      !! ** Method  : pile-up ice when no ridging/rafting
211      !!
212      !! ** input   : a_i
213      !!-------------------------------------------------------------------
214      INTEGER ::   jl         ! dummy loop indices
215      !!-------------------------------------------------------------------
216      !
217      CALL ice_var_zapsmall                       !-- zap small areas
218      !
219      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
220      DO jl = 1, jpl
221         WHERE( at_i(:,:) > epsi20 )
222            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
223         END WHERE
224      END DO
225      !
226   END SUBROUTINE Hpiling
227
228
229   SUBROUTINE ice_dyn_init
230      !!-------------------------------------------------------------------
231      !!                  ***  ROUTINE ice_dyn_init  ***
232      !!
233      !! ** Purpose : Physical constants and parameters linked to the ice
234      !!      dynamics
235      !!
236      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
237      !!       parameter values called at the first timestep (nit000)
238      !!
239      !! ** input   :   Namelist namdyn
240      !!-------------------------------------------------------------------
241      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
242      !!
243      NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
244         &             rn_ishlat ,                                             &
245         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
246      !!-------------------------------------------------------------------
247      !
248      REWIND( numnam_ice_ref )         ! Namelist namdyn in reference namelist : Ice dynamics
249      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
250901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn in reference namelist', lwp )
251      REWIND( numnam_ice_cfg )         ! Namelist namdyn in configuration namelist : Ice dynamics
252      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
253902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn in configuration namelist', lwp )
254      IF(lwm) WRITE( numoni, namdyn )
255      !
256      IF(lwp) THEN                     ! control print
257         WRITE(numout,*)
258         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
259         WRITE(numout,*) '~~~~~~~~~~~~'
260         WRITE(numout,*) '   Namelist namdyn:'
261         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL      = ', ln_dynFULL
262         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV
263         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV       = ', ln_dynADV
264         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')'
265         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat
266         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16
267         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home
268         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra
269         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr
270         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax
271         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile
272         WRITE(numout,*)
273      ENDIF
274      !                             !== set the choice of ice dynamics ==!
275      ioptio = 0 
276      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
277      IF( ln_dynFULL   ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF
278      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
279      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
280      !      !--- advection only with prescribed ice velocities (from namelist)
281      IF( ln_dynADV    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF
282      !
283      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
284      !
285      !                                      !--- Lateral boundary conditions
286      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
287      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
288      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
289      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
290      ENDIF
291      !                                      !--- Landfast ice
292      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp
293      !
294      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN
295         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' )
296      ENDIF
297      !
298      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
299      CALL ice_dyn_rhg_init             ! set ice rheology parameters
300      CALL ice_dyn_adv_init             ! set ice advection parameters
301      !
302   END SUBROUTINE ice_dyn_init
303
304#else
305   !!----------------------------------------------------------------------
306   !!   Default option         Empty module           NO SI3 sea-ice model
307   !!----------------------------------------------------------------------
308#endif 
309
310   !!======================================================================
311END MODULE icedyn
Note: See TracBrowser for help on using the repository browser.