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/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE – NEMO

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/ICE/icedyn.F90 @ 13159

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

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 14.3 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 "do_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        ! dummy loop indices
78      REAL(wp) ::   zcoefu, zcoefv
79      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdivu_i
80      !!--------------------------------------------------------------------
81      !
82      ! controls
83      IF( ln_timing )   CALL timing_start('icedyn')
84      !
85      IF( kt == nit000 .AND. lwp ) THEN
86         WRITE(numout,*)
87         WRITE(numout,*)'ice_dyn: sea-ice dynamics'
88         WRITE(numout,*)'~~~~~~~'
89      ENDIF
90      !                     
91      ! retrieve thickness from volume for landfast param. and UMx advection scheme
92      WHERE( a_i(:,:,:) >= epsi20 )
93         h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
94         h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
95      ELSEWHERE
96         h_i(:,:,:) = 0._wp
97         h_s(:,:,:) = 0._wp
98      END WHERE
99      !
100      WHERE( a_ip(:,:,:) >= epsi20 )
101         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
102      ELSEWHERE
103         h_ip(:,:,:) = 0._wp
104      END WHERE
105      !
106      !
107      SELECT CASE( nice_dyn )          !-- Set which dynamics is running
108
109      CASE ( np_dynALL )           !==  all dynamical processes  ==!
110         !
111         CALL ice_dyn_rhg   ( kt, Kmm )                                     ! -- rheology 
112         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
113         CALL ice_dyn_rdgrft( kt )                                          ! -- ridging/rafting
114         CALL ice_cor       ( kt , 1 )                                      ! -- Corrections
115         !
116      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==!
117         !
118         CALL ice_dyn_rhg   ( kt, Kmm )                                     ! -- rheology 
119         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
120         CALL Hpiling                                                       ! -- simple pile-up (replaces ridging/rafting)
121         CALL ice_var_zapsmall                                              ! -- zap small areas
122         !
123      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D)
124         !
125         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
126         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
127         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s
128         DO_2D_11_11
129            zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )
130            zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )
131            u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
132            v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
133         END_2D
134         ! ---
135         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
136         !
137      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities)
138         !
139         u_ice(:,:) = rn_uice * umask(:,:,1)
140         v_ice(:,:) = rn_vice * vmask(:,:,1)
141         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
142         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
143         ! ---
144         CALL ice_dyn_adv   ( kt )                                          ! -- advection of ice
145
146      END SELECT
147      !
148      !
149      ! diagnostics: divergence at T points
150      IF( iom_use('icediv') ) THEN
151         !
152         SELECT CASE( nice_dyn )
153
154         CASE ( np_dynADV1D , np_dynADV2D )
155
156            ALLOCATE( zdivu_i(jpi,jpj) )
157            DO_2D_00_00
158               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
159                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
160            END_2D
161            CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. )
162            ! output
163            CALL iom_put( 'icediv' , zdivu_i )
164
165            DEALLOCATE( zdivu_i )
166
167         END SELECT
168         !
169      ENDIF
170      !
171      ! controls
172      IF( ln_timing )   CALL timing_stop ('icedyn')
173      !
174   END SUBROUTINE ice_dyn
175
176
177   SUBROUTINE Hpiling
178      !!-------------------------------------------------------------------
179      !!                  ***  ROUTINE Hpiling  ***
180      !!
181      !! ** Purpose : Simple conservative piling comparable with 1-cat models
182      !!
183      !! ** Method  : pile-up ice when no ridging/rafting
184      !!
185      !! ** input   : a_i
186      !!-------------------------------------------------------------------
187      INTEGER ::   jl         ! dummy loop indices
188      !!-------------------------------------------------------------------
189      ! controls
190      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
191      !
192      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
193      DO jl = 1, jpl
194         WHERE( at_i(:,:) > epsi20 )
195            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
196         END WHERE
197      END DO
198      ! controls
199      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
200      !
201   END SUBROUTINE Hpiling
202
203
204   SUBROUTINE ice_dyn_init
205      !!-------------------------------------------------------------------
206      !!                  ***  ROUTINE ice_dyn_init  ***
207      !!
208      !! ** Purpose : Physical constants and parameters linked to the ice
209      !!      dynamics
210      !!
211      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
212      !!       parameter values called at the first timestep (nit000)
213      !!
214      !! ** input   :   Namelist namdyn
215      !!-------------------------------------------------------------------
216      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
217      !!
218      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  &
219         &             rn_ishlat ,                                                           &
220         &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
221      !!-------------------------------------------------------------------
222      !
223      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
224901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn in reference namelist' )
225      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
226902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn in configuration namelist' )
227      IF(lwm) WRITE( numoni, namdyn )
228      !
229      IF(lwp) THEN                     ! control print
230         WRITE(numout,*)
231         WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
232         WRITE(numout,*) '~~~~~~~~~~~~'
233         WRITE(numout,*) '   Namelist namdyn:'
234         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr) ln_dynALL       = ', ln_dynALL
235         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                     ln_dynRHGADV    = ', ln_dynRHGADV
236         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)    ln_dynADV1D     = ', ln_dynADV1D
237         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                ln_dynADV2D     = ', ln_dynADV2D
238         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)   = (', rn_uice,',',rn_vice,')'
239         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat
240         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16
241         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra
242         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr
243         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax
244         WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile
245         WRITE(numout,*)
246      ENDIF
247      !                             !== set the choice of ice dynamics ==!
248      ioptio = 0 
249      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
250      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF
251      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
252      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF
253      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
254      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF
255      !      !--- advection 2D only with prescribed ice velocities (from namelist)
256      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF
257      !
258      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
259      !
260      !                                      !--- Lateral boundary conditions
261      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip'
262      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip'
263      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip'
264      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
265      ENDIF
266      !                                      !--- Landfast ice
267      IF( .NOT.ln_landfast_L16 )   tau_icebfr(:,:) = 0._wp
268      !
269      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
270      CALL ice_dyn_rhg_init             ! set ice rheology parameters
271      CALL ice_dyn_adv_init             ! set ice advection parameters
272      !
273   END SUBROUTINE ice_dyn_init
274
275#else
276   !!----------------------------------------------------------------------
277   !!   Default option         Empty module           NO SI3 sea-ice model
278   !!----------------------------------------------------------------------
279#endif 
280
281   !!======================================================================
282END MODULE icedyn
Note: See TracBrowser for help on using the repository browser.