Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
MODULE icedyn
!!======================================================================
!! *** MODULE icedyn ***
!! Sea-Ice dynamics : master routine for sea ice dynamics
!!======================================================================
!! history : 4.0 ! 2018 (C. Rousset) original code SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_dyn : dynamics of sea ice
!! ice_dyn_init : initialization and namelist read
!!----------------------------------------------------------------------
USE phycst ! physical constants
USE dom_oce ! ocean space and time domain
USE ice ! sea-ice: variables
USE icedyn_rhg ! sea-ice: rheology
USE icedyn_adv ! sea-ice: advection
USE icedyn_rdgrft ! sea-ice: ridging/rafting
USE icecor ! sea-ice: corrections
USE icevar ! sea-ice: operations
USE icectl ! sea-ice: control prints
!
USE in_out_manager ! I/O manager
USE iom ! I/O manager library
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE lbclnk ! lateral boundary conditions (or mpp links)
USE timing ! Timing
USE fldread ! read input fields
IMPLICIT NONE
PRIVATE
PUBLIC ice_dyn ! called by icestp.F90
PUBLIC ice_dyn_init ! called by icestp.F90
INTEGER :: nice_dyn ! choice of the type of dynamics
! ! associated indices:
INTEGER, PARAMETER :: np_dynALL = 1 ! full ice dynamics (rheology + advection + ridging/rafting + correction)
INTEGER, PARAMETER :: np_dynRHGADV = 2 ! pure dynamics (rheology + advection)
INTEGER, PARAMETER :: np_dynADV1D = 3 ! only advection 1D - test case from Schar & Smolarkiewicz 1996
INTEGER, PARAMETER :: np_dynADV2D = 4 ! only advection 2D w prescribed vel.(rn_uvice + advection)
!
! ** namelist (namdyn) **
LOGICAL :: ln_dynALL ! full ice dynamics (rheology + advection + ridging/rafting + correction)
LOGICAL :: ln_dynRHGADV ! no ridge/raft & no corrections (rheology + advection)
LOGICAL :: ln_dynADV1D ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996)
LOGICAL :: ln_dynADV2D ! only advection in 2D w prescribed vel. (rn_uvice + advection)
REAL(wp) :: rn_uice ! prescribed u-vel (case np_dynADV1D & np_dynADV2D)
REAL(wp) :: rn_vice ! prescribed v-vel (case np_dynADV1D & np_dynADV2D)
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_icbmsk ! structure of input grounded icebergs mask (file informations, fields read)
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icedyn.F90 14997 2021-06-16 06:43:57Z smasson $
!! Software governed by the CeCILL licence (./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_dyn( kt, Kmm )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_dyn ***
!!
!! ** Purpose : this routine manages sea ice dynamics
!!
!! ** Action : - calculation of friction in case of landfast ice
!! - call ice_dyn_rhg = rheology
!! - call ice_dyn_adv = advection
!! - call ice_dyn_rdgrft = ridging/rafting
!! - call ice_cor = corrections if fields are out of bounds
!!--------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ice time step
INTEGER, INTENT(in) :: Kmm ! ocean time level index
!!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zcoefu, zcoefv
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i
!!--------------------------------------------------------------------
!
! controls
IF( ln_timing ) CALL timing_start('ice_dyn')
!
IF( kt == nit000 .AND. lwp ) THEN
WRITE(numout,*)
WRITE(numout,*)'ice_dyn: sea-ice dynamics'
WRITE(numout,*)'~~~~~~~'
ENDIF
!
! retrieve thickness from volume for landfast param. and UMx advection scheme
WHERE( a_i(:,:,:) >= epsi20 )
h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:)
h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:)
ELSEWHERE
h_i(:,:,:) = 0._wp
h_s(:,:,:) = 0._wp
END WHERE
!
WHERE( a_ip(:,:,:) >= epsi20 )
h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:)
h_il(:,:,:) = v_il(:,:,:) / a_ip(:,:,:)
ELSEWHERE
h_ip(:,:,:) = 0._wp
h_il(:,:,:) = 0._wp
END WHERE
!
IF( ln_landfast_L16 ) THEN
CALL fld_read( kt, 1, sf_icbmsk )
icb_mask(:,:) = sf_icbmsk(1)%fnow(:,:,1)
ENDIF
!
SELECT CASE( nice_dyn ) !-- Set which dynamics is running
CASE ( np_dynALL ) !== all dynamical processes ==!
!
CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology
CALL ice_dyn_adv ( kt ) ! -- advection of ice
CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting
CALL ice_cor ( kt , 1 ) ! -- Corrections
!
CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==!
!
CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology
CALL ice_dyn_adv ( kt ) ! -- advection of ice
CALL Hpiling ! -- simple pile-up (replaces ridging/rafting)
CALL ice_var_zapsmall ! -- zap small areas
!
CASE ( np_dynADV1D ) !== pure advection ==! (1D)
!
! --- monotonicity test from Schar & Smolarkiewicz 1996 --- !
! CFL = 0.5 at a distance from the bound of 1/6 of the basin length
! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zcoefu = ( REAL(jpiglo+1)*0.5_wp - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5_wp - 1._wp )
zcoefv = ( REAL(jpjglo+1)*0.5_wp - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5_wp - 1._wp )
u_ice(ji,jj) = rn_uice * 1.5_wp * SIGN( 1.0_wp, zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)
v_ice(ji,jj) = rn_vice * 1.5_wp * SIGN( 1.0_wp, zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)
END_2D
! ---
CALL ice_dyn_adv ( kt ) ! -- advection of ice
!
CASE ( np_dynADV2D ) !== pure advection ==! (2D w prescribed velocities)
!
u_ice(:,:) = rn_uice * umask(:,:,1)
v_ice(:,:) = rn_vice * vmask(:,:,1)
!CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1)
!CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1)
! ---
CALL ice_dyn_adv ( kt ) ! -- advection of ice
END SELECT
!
!
! diagnostics: divergence at T points
IF( iom_use('icediv') ) THEN
!
SELECT CASE( nice_dyn )
CASE ( np_dynADV1D , np_dynADV2D )
ALLOCATE( zdivu_i(jpi,jpj) )
DO_2D( 0, 0, 0, 0 )
zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
END_2D
CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1.0_wp )
! output
CALL iom_put( 'icediv' , zdivu_i )
DEALLOCATE( zdivu_i )
END SELECT
!
ENDIF
!
! controls
IF( ln_timing ) CALL timing_stop ('ice_dyn')
!
END SUBROUTINE ice_dyn
SUBROUTINE Hpiling
!!-------------------------------------------------------------------
!! *** ROUTINE Hpiling ***
!!
!! ** Purpose : Simple conservative piling comparable with 1-cat models
!!
!! ** Method : pile-up ice when no ridging/rafting
!!
!! ** input : a_i
!!-------------------------------------------------------------------
INTEGER :: jl ! dummy loop indices
!!-------------------------------------------------------------------
! controls
IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
!
at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
DO jl = 1, jpl
WHERE( at_i(:,:) > epsi20 )
a_i(:,:,jl) = a_i(:,:,jl) * ( 1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:) )
END WHERE
END DO
! controls
IF( ln_icediachk ) CALL ice_cons_hsm(1, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
!
END SUBROUTINE Hpiling
SUBROUTINE ice_dyn_init
!!-------------------------------------------------------------------
!! *** ROUTINE ice_dyn_init ***
!!
!! ** Purpose : Physical constants and parameters linked to the ice
!! dynamics
!!
!! ** Method : Read the namdyn namelist and check the ice-dynamic
!! parameter values called at the first timestep (nit000)
!!
!! ** input : Namelist namdyn
!!-------------------------------------------------------------------
INTEGER :: ios, ioptio, ierror ! Local integer output status for namelist read
!
CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
TYPE(FLD_N) :: sn_icbmsk ! informations about the grounded icebergs field to be read
!!
NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice, &
& rn_ishlat , &
& ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile, &
& sn_icbmsk, cn_dir
!!-------------------------------------------------------------------
!
READ ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn in reference namelist' )
READ ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn in configuration namelist' )
IF(lwm) WRITE( numoni, namdyn )
!
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
WRITE(numout,*) '~~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namdyn:'
WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL
WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV
WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D
WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D
WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')'
WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat
WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16
WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_lf_depfra = ', rn_lf_depfra
WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_lf_bfr = ', rn_lf_bfr
WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lf_relax = ', rn_lf_relax
WRITE(numout,*) ' isotropic tensile strength rn_lf_tensile = ', rn_lf_tensile
WRITE(numout,*)
ENDIF
! !== set the choice of ice dynamics ==!
ioptio = 0
! !--- full dynamics (rheology + advection + ridging/rafting + correction)
IF( ln_dynALL ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynALL ; ENDIF
! !--- dynamics without ridging/rafting and corr (rheology + advection)
IF( ln_dynRHGADV ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynRHGADV ; ENDIF
! !--- advection 1D only - test case from Schar & Smolarkiewicz 1996
IF( ln_dynADV1D ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynADV1D ; ENDIF
! !--- advection 2D only with prescribed ice velocities (from namelist)
IF( ln_dynADV2D ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynADV2D ; ENDIF
!
IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
!
! !--- Lateral boundary conditions
IF ( rn_ishlat == 0. ) THEN ; IF(lwp) WRITE(numout,*) ' ===>>> ice lateral free-slip'
ELSEIF ( rn_ishlat == 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ===>>> ice lateral no-slip'
ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN ; IF(lwp) WRITE(numout,*) ' ===>>> ice lateral partial-slip'
ELSEIF ( 2. < rn_ishlat ) THEN ; IF(lwp) WRITE(numout,*) ' ===>>> ice lateral strong-slip'
ENDIF
! !--- Landfast ice
IF( .NOT.ln_landfast_L16 ) tau_icebfr(:,:) = 0._wp
!
! !--- allocate and fill structure for grounded icebergs mask
IF( ln_landfast_L16 ) THEN
ALLOCATE( sf_icbmsk(1), STAT=ierror )
IF( ierror > 0 ) THEN
CALL ctl_stop( 'ice_dyn_init: unable to allocate sf_icbmsk structure' ) ; RETURN
ENDIF
!
CALL fld_fill( sf_icbmsk, (/ sn_icbmsk /), cn_dir, 'ice_dyn_init', &
& 'landfast ice is a function of read grounded icebergs', 'icedyn' )
!
ALLOCATE( sf_icbmsk(1)%fnow(jpi,jpj,1) )
IF( sf_icbmsk(1)%ln_tint ) ALLOCATE( sf_icbmsk(1)%fdta(jpi,jpj,1,2) )
IF( TRIM(sf_icbmsk(1)%clrootname) == 'NOT USED' ) sf_icbmsk(1)%fnow(:,:,1) = 0._wp ! not used field (set to 0)
ELSE
icb_mask(:,:) = 0._wp
ENDIF
! !--- other init
CALL ice_dyn_rdgrft_init ! set ice ridging/rafting parameters
CALL ice_dyn_rhg_init ! set ice rheology parameters
CALL ice_dyn_adv_init ! set ice advection parameters
!
END SUBROUTINE ice_dyn_init
#else
!!----------------------------------------------------------------------
!! Default option Empty module NO SI3 sea-ice model
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE icedyn