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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/ICE – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/ICE/icedyn.F90 @ 11463

Last change on this file since 11463 was 11053, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Merge in latest changes from main branch and finish conversion of "h" variables. NB. This version still doesn't work!
  • Property svn:keywords set to Id
File size: 22.0 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   USE icectl         ! sea-ice: control prints
24   !
25   USE in_out_manager ! I/O manager
26   USE iom            ! I/O manager library
27   USE lib_mpp        ! MPP library
28   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
29   USE lbclnk         ! lateral boundary conditions (or mpp links)
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   ice_dyn        ! called by icestp.F90
36   PUBLIC   ice_dyn_init   ! called by icestp.F90
37   
38   INTEGER ::              nice_dyn   ! choice of the type of dynamics
39   !                                        ! associated indices:
40   INTEGER, PARAMETER ::   np_dynALL     = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
41   INTEGER, PARAMETER ::   np_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection)
42   INTEGER, PARAMETER ::   np_dynADV1D   = 3   ! only advection 1D - test case from Schar & Smolarkiewicz 1996
43   INTEGER, PARAMETER ::   np_dynADV2D   = 4   ! only advection 2D w prescribed vel.(rn_uvice + advection)
44   !
45   ! ** namelist (namdyn) **
46   LOGICAL  ::   ln_dynALL        ! full ice dynamics                      (rheology + advection + ridging/rafting + correction)
47   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections         (rheology + advection)
48   LOGICAL  ::   ln_dynADV1D      ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996)
49   LOGICAL  ::   ln_dynADV2D      ! only advection in 2D w prescribed vel. (rn_uvice + advection)
50   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV1D & np_dynADV2D)
51   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV1D & np_dynADV2D)
52   
53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
57   !! $Id$
58   !! Software governed by the CeCILL licence     (./LICENSE)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE ice_dyn( kt, Kmm )
63      !!-------------------------------------------------------------------
64      !!               ***  ROUTINE ice_dyn  ***
65      !!               
66      !! ** Purpose :   this routine manages sea ice dynamics
67      !!
68      !! ** Action : - calculation of friction in case of landfast ice
69      !!             - call ice_dyn_rhg    = rheology
70      !!             - call ice_dyn_adv    = advection
71      !!             - call ice_dyn_rdgrft = ridging/rafting
72      !!             - call ice_cor        = corrections if fields are out of bounds
73      !!--------------------------------------------------------------------
74      INTEGER, INTENT(in) ::   kt     ! ice time step
75      INTEGER, INTENT(in) ::   Kmm    ! ocean time level index
76      !!
77      INTEGER  ::   ji, jj, jl        ! dummy loop indices
78      REAL(wp) ::   zcoefu, zcoefv
79      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max
80      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i
81      !!--------------------------------------------------------------------
82      !
83      ! controls
84      IF( ln_timing )   CALL timing_start('icedyn')
85      !
86      IF( kt == nit000 .AND. lwp ) THEN
87         WRITE(numout,*)
88         WRITE(numout,*)'ice_dyn: sea-ice dynamics'
89         WRITE(numout,*)'~~~~~~~'
90      ENDIF
91      !                     
92      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
93         tau_icebfr(:,:) = 0._wp
94         DO jl = 1, jpl
95            WHERE( h_i_b(:,:,jl) > ht(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
96         END DO
97      ENDIF
98      !
99      !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
100      DO jl = 1, jpl
101         DO jj = 2, jpjm1
102            DO ji = fs_2, fs_jpim1
103               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), &
104                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), &
105                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), &
106                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) )
107               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), &
108                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), &
109                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), &
110                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) )
111               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), &
112                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), &
113                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), &
114                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) )
115            END DO
116         END DO
117      END DO
118      CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. )
119      !
120      !
121      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
122
123      CASE ( np_dynALL )           !==  all dynamical processes  ==!
124         CALL ice_dyn_rhg   ( kt, Kmm )                                            ! -- rheology 
125         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
126         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting
127         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections
128
129      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
130         CALL ice_dyn_rhg   ( kt, Kmm )                                            ! -- rheology 
131         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness
132         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting)
133
134      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
135         ALLOCATE( zdivu_i(jpi,jpj) )
136         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
137         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
138         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s
139         DO jj = 1, jpj
140            DO ji = 1, jpi
141               zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
142               zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
143               u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
144               v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
145            END DO
146         END DO
147         ! ---
148         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
149
150         ! diagnostics: divergence at T points
151         DO jj = 2, jpjm1
152            DO ji = 2, jpim1
153               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
154                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
155            END DO
156         END DO
157         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
158         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
159
160         DEALLOCATE( zdivu_i )
161
162      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
163         ALLOCATE( zdivu_i(jpi,jpj) )
164         u_ice(:,:) = rn_uice * umask(:,:,1)
165         v_ice(:,:) = rn_vice * vmask(:,:,1)
166         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
167         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
168         ! ---
169         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice
170
171         ! diagnostics: divergence at T points
172         DO jj = 2, jpjm1
173            DO ji = 2, jpim1
174               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
175                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
176            END DO
177         END DO
178         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
179         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) )
180
181         DEALLOCATE( zdivu_i )
182
183      END SELECT
184      !
185      ! controls
186      IF( ln_timing )   CALL timing_stop ('icedyn')
187      !
188   END SUBROUTINE ice_dyn
189
190
191   SUBROUTINE Hbig( phi_max, phs_max, phip_max )
192      !!-------------------------------------------------------------------
193      !!                  ***  ROUTINE Hbig  ***
194      !!
195      !! ** Purpose : Thickness correction in case advection scheme creates
196      !!              abnormally tick ice or snow
197      !!
198      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
199      !!                 (before advection) and reduce it by adapting ice concentration
200      !!              2- check whether snow thickness is larger than the surrounding 9-points
201      !!                 (before advection) and reduce it by sending the excess in the ocean
202      !!              3- check whether snow load deplets the snow-ice interface below sea level$
203      !!                 and reduce it by sending the excess in the ocean
204      !!              4- correct pond fraction to avoid a_ip > a_i
205      !!
206      !! ** input   : Max thickness of the surrounding 9-points
207      !!-------------------------------------------------------------------
208      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts
209      !
210      INTEGER  ::   ji, jj, jl         ! dummy loop indices
211      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra
212      !!-------------------------------------------------------------------
213      ! controls
214      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
215      !
216      CALL ice_var_zapsmall                       !-- zap small areas
217      !
218      DO jl = 1, jpl
219         DO jj = 1, jpj
220            DO ji = 1, jpi
221               IF ( v_i(ji,jj,jl) > 0._wp ) THEN
222                  !
223                  !                               ! -- check h_ip -- !
224                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
225                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN
226                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )
227                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN
228                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)
229                     ENDIF
230                  ENDIF
231                  !
232                  !                               ! -- check h_i -- !
233                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
234                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
235                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
236                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m)
237                  ENDIF
238                  !
239                  !                               ! -- check h_s -- !
240                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
241                  zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
242                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN
243                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
244                     !
245                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice
246                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
247                     !
248                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
249                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl)
250                  ENDIF           
251                  !
252                  !                               ! -- check snow load -- !
253                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean
254                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)
255                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically)
256                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
257                  IF( zvs_excess > 0._wp ) THEN
258                     zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 )
259                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice
260                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0
261                     !
262                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra
263                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess
264                  ENDIF
265                 
266               ENDIF
267            END DO
268         END DO
269      END DO 
270      !                                           !-- correct pond fraction to avoid a_ip > a_i
271      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
272      !
273      ! controls
274      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
275      !
276   END SUBROUTINE Hbig
277
278
279   SUBROUTINE Hpiling
280      !!-------------------------------------------------------------------
281      !!                  ***  ROUTINE Hpiling  ***
282      !!
283      !! ** Purpose : Simple conservative piling comparable with 1-cat models
284      !!
285      !! ** Method  : pile-up ice when no ridging/rafting
286      !!
287      !! ** input   : a_i
288      !!-------------------------------------------------------------------
289      INTEGER ::   jl         ! dummy loop indices
290      !!-------------------------------------------------------------------
291      ! controls
292      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
293      !
294      CALL ice_var_zapsmall                       !-- zap small areas
295      !
296      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
297      DO jl = 1, jpl
298         WHERE( at_i(:,:) > epsi20 )
299            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
300         END WHERE
301      END DO
302      ! controls
303      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
304      !
305   END SUBROUTINE Hpiling
306
307
308   SUBROUTINE ice_dyn_init
309      !!-------------------------------------------------------------------
310      !!                  ***  ROUTINE ice_dyn_init  ***
311      !!
312      !! ** Purpose : Physical constants and parameters linked to the ice
313      !!      dynamics
314      !!
315      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
316      !!       parameter values called at the first timestep (nit000)
317      !!
318      !! ** input   :   Namelist namdyn
319      !!-------------------------------------------------------------------
320      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
321      !!
322      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
323         &             rn_ishlat ,                                                           &
324         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
325      !!-------------------------------------------------------------------
326      !
327      REWIND( numnam_ice_ref )         ! Namelist namdyn in reference namelist : Ice dynamics
328      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
329901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn in reference namelist', lwp )
330      REWIND( numnam_ice_cfg )         ! Namelist namdyn in configuration namelist : Ice dynamics
331      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
332902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn in configuration namelist', lwp )
333      IF(lwm) WRITE( numoni, namdyn )
334      !
335      IF(lwp) THEN                     ! control print
336         WRITE(numout,*)
337         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
338         WRITE(numout,*) '~~~~~~~~~~~~'
339         WRITE(numout,*) '   Namelist namdyn:'
340         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL
341         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV
342         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D
343         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D
344         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')'
345         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat
346         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16
347         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home
348         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra
349         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr
350         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax
351         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile
352         WRITE(numout,*)
353      ENDIF
354      !                             !== set the choice of ice dynamics ==!
355      ioptio = 0 
356      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
357      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF
358      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
359      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
360      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
361      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF
362      !      !--- advection 2D only with prescribed ice velocities (from namelist)
363      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF
364      !
365      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
366      !
367      !                                      !--- Lateral boundary conditions
368      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
369      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
370      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
371      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
372      ENDIF
373      !                                      !--- Landfast ice
374      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp
375      !
376      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN
377         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' )
378      ENDIF
379      !
380      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
381      CALL ice_dyn_rhg_init             ! set ice rheology parameters
382      CALL ice_dyn_adv_init             ! set ice advection parameters
383      !
384   END SUBROUTINE ice_dyn_init
385
386#else
387   !!----------------------------------------------------------------------
388   !!   Default option         Empty module           NO SI3 sea-ice model
389   !!----------------------------------------------------------------------
390#endif 
391
392   !!======================================================================
393END MODULE icedyn
Note: See TracBrowser for help on using the repository browser.