ROMS
Loading...
Searching...
No Matches
ice_evp.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE ice_evp_mod
3
4#if defined ICE_MOMENTUM && defined ICE_EVP
5!
6!git $Id$
7!=======================================================================
8! Copyright (c) 2002-2025 The ROMS Group !
9! Licensed under a MIT/X style license W. Paul Budgell !
10! See License_ROMS.md Katherine Hedstrom !
11!================================================== Hernan G. Arango ===
12! !
13! This module computes fields for elastic-viscous-plastic rheology. !
14! It load elastic ice velocities (uie,vie) and couputes ice pressure !
15! (icPice), bulk viscosity (icBvis), and shear viscosity (icSvis). !
16! !
17!=======================================================================
18!
19 USE mod_param
20 USE mod_grid
21 USE mod_ice
22 USE mod_scalars
23!
24 USE bc_2d_mod, ONLY : bc_r2d_tile
25# ifdef DISTRIBUTE
26 USE mp_exchange_mod, ONLY : mp_exchange2d
27# endif
28!
29 implicit none
30!
31 PRIVATE
32 PUBLIC :: ice_evp
33!
34 CONTAINS
35!
36!***********************************************************************
37 SUBROUTINE ice_evp (ng, tile, model)
38!***********************************************************************
39!
40 USE mod_stepping
41!
42! Imported variable declarations.
43!
44 integer, intent(in) :: ng, tile, model
45!
46! Local variable declarations.
47!
48 character (len=*), parameter :: MyFile = &
49 & __FILE__
50!
51# include "tile.h"
52!
53# ifdef PROFILE
54 CALL wclock_on (ng, model, 42, __line__, myfile)
55# endif
56 CALL ice_evp_tile (ng, tile, model, &
57 & lbi, ubi, lbj, ubj, &
58 & imins, imaxs, jmins, jmaxs, &
59 & liold(ng), lieol(ng), liuol(ng), &
60# ifdef OUTFLOW_MASK
61 & grid(ng) % mask_outflow, &
62# endif
63 & grid(ng) % pm, &
64 & grid(ng) % pn, &
65 & ice(ng) % Fi, &
66 & ice(ng) % Si)
67# ifdef PROFILE
68 CALL wclock_off (ng, model, 42, __line__, myfile)
69# endif
70!
71 RETURN
72 END SUBROUTINE ice_evp
73!
74!***********************************************************************
75 SUBROUTINE ice_evp_tile (ng, tile, model, &
76 & LBi, UBi, LBj, UBj, &
77 & IminS, ImaxS, JminS, JmaxS, &
78 & liold, lieol, liuol, &
79# ifdef OUTFLOW_MASK
80 & mask_outflow, &
81# endif
82 & pm, pn, &
83 & Fi, Si)
84!***********************************************************************
85!
86! Imported variable declarations.
87!
88 integer, intent(in) :: ng, tile, model
89 integer, intent(in) :: LBi, UBi, LBj, UBj
90 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
91 integer, intent(in) :: liold, lieol, liuol
92!
93# ifdef ASSUMED_SHAPE
94# ifdef OUTFLOW_MASK
95 real(r8), intent(in) :: mask_outflow(LBi:,LBj:)
96# endif
97 real(r8), intent(in) :: pm(LBi:,LBj:)
98 real(r8), intent(in) :: pn(LBi:,LBj:)
99 real(r8), intent(inout) :: Fi(LBi:,LBj:,:)
100 real(r8), intent(inout) :: Si(LBi:,LBj:,:,:)
101# else
102# ifdef OUTFLOW_MASK
103 real(r8), intent(in) :: mask_outflow(LBi:UBi,LBj:UBj)
104# endif
105 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
106 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
107 real(r8), intent(inout) :: Fi(LBi:UBi,LBj:UBj,nIceF)
108 real(r8), intent(inout) :: Si(LBi:UBi,LBj:UBj,2,nIceS)
109# endif
110!
111! Local variable definitions
112!
113 integer :: i, j
114!
115 real(r8) :: delta, e2r, eone, epx, epy, etwos, zmax
116
117 real(r8), parameter :: epso = 1.0e-12_r8
118!
119 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: eps_xx
120 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: eps_yy
121 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: eps_xy
122
123# include "set_bounds.h"
124!
125!-----------------------------------------------------------------------
126! Compute elastic terms as function of the stress tensor and rates of
127! strain.
128!-----------------------------------------------------------------------
129!
130! Load ice velocity components into EVP velocity arrays.
131!
132 IF (ievp(ng).eq.1) THEN
133 DO j=jstr,jend
134 DO i=istrp,iend
135 si(i,j,1,isuevp)=si(i,j,liuol,isuice)
136 si(i,j,2,isuevp)=si(i,j,liuol,isuice)
137 END DO
138 END DO
139 DO j=jstrp,jend
140 DO i=istr,iend
141 si(i,j,1,isvevp)=si(i,j,liuol,isvice)
142 si(i,j,2,isvevp)=si(i,j,liuol,isvice)
143 END DO
144 END DO
145!
146 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
147 IF (domain(ng)%Western_Edge(tile)) THEN
148 DO j=jstr,jend
149 si(istr,j,1,isuevp)=si(istr,j,liuol,isuice)
150 si(istr,j,2,isuevp)=si(istr,j,liuol,isuice)
151 END DO
152 DO j=jstrp,jend
153 si(istr-1,j,1,isvevp)=si(istr-1,j,liuol,isvice)
154 si(istr-1,j,2,isvevp)=si(istr-1,j,liuol,isvice)
155 END DO
156 END IF
157 END IF
158!
159 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
160 IF (domain(ng)%Eastern_Edge(tile)) THEN
161 DO j=jstr,jend
162 si(iend+1,j,1,isuevp)=si(iend+1,j,liuol,isuice)
163 si(iend+1,j,2,isuevp)=si(iend+1,j,liuol,isuice)
164 END DO
165 DO j=jstrp,jend
166 si(iend+1,j,1,isvevp)=si(iend+1,j,liuol,isvice)
167 si(iend+1,j,2,isvevp)=si(iend+1,j,liuol,isvice)
168 END DO
169 END IF
170 END IF
171!
172 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
173 IF (domain(ng)%Southern_Edge(tile)) THEN
174 DO i=istrp,iend
175 si(i,jstr-1,1,isuevp)=si(i,jstr-1,liuol,isuice)
176 si(i,jstr-1,2,isuevp)=si(i,jstr-1,liuol,isuice)
177 END DO
178 DO i=istr,iend
179 si(i,jstr,1,isvevp)=si(i,jstr,liuol,isvice)
180 si(i,jstr,2,isvevp)=si(i,jstr,liuol,isvice)
181 END DO
182 END IF
183 END IF
184!
185 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
186 IF (domain(ng)%Northern_Edge(tile)) THEN
187 DO i=istrp,iend
188 si(i,jend+1,1,isuevp)=si(i,jend+1,liuol,isuice)
189 si(i,jend+1,2,isuevp)=si(i,jend+1,liuol,isuice)
190 END DO
191 DO i=istr,iend
192 si(i,jend+1,1,isvevp)=si(i,jend+1,liuol,isvice)
193 si(i,jend+1,2,isvevp)=si(i,jend+1,liuol,isvice)
194 END DO
195 END IF
196 END IF
197!
198 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
199 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
200 IF (domain(ng)%SouthWest_Corner(tile)) THEN
201 si(istr ,jstr-1,1,isuevp)=si(istr ,jstr-1,liuol,isuice)
202 si(istr ,jstr-1,2,isuevp)=si(istr ,jstr-1,liuol,isuice)
203 si(istr-1,jstr ,1,isvevp)=si(istr-1,jstr ,liuol,isvice)
204 si(istr-1,jstr ,2,isvevp)=si(istr-1,jstr ,liuol,isvice)
205 END IF
206 END IF
207!
208 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
209 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
210 IF (domain(ng)%SouthEast_Corner(tile)) THEN
211 si(iend+1,jstr-1,1,isuevp)=si(iend+1,jstr-1,liuol,isuice)
212 si(iend+1,jstr-1,2,isuevp)=si(iend+1,jstr-1,liuol,isuice)
213 si(iend+1,jstr ,1,isvevp)=si(iend+1,jstr ,liuol,isvice)
214 si(iend+1,jstr ,2,isvevp)=si(iend+1,jstr ,liuol,isvice)
215 END IF
216 END IF
217!
218 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
219 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
220 IF (domain(ng)%NorthWest_Corner(tile)) THEN
221 si(istr ,jend+1,1,isuevp)=si(istr ,jend+1,liuol,isuice)
222 si(istr ,jend+1,2,isuevp)=si(istr ,jend+1,liuol,isuice)
223 si(istr-1,jend+1,1,isvevp)=si(istr-1,jend+1,liuol,isvice)
224 si(istr-1,jend+1,2,isvevp)=si(istr-1,jend+1,liuol,isvice)
225 END IF
226 END IF
227!
228 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
229 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
230 IF (domain(ng)%NorthEast_Corner(tile)) THEN
231 si(iend+1,jend+1,1,isuevp)=si(iend+1,jend+1,liuol,isuice)
232 si(iend+1,jend+1,2,isuevp)=si(iend+1,jend+1,liuol,isuice)
233 si(iend+1,jend+1,1,isvevp)=si(iend+1,jend+1,liuol,isvice)
234 si(iend+1,jend+1,2,isvevp)=si(iend+1,jend+1,liuol,isvice)
235 END IF
236 END IF
237 END IF
238
239# ifdef DISTRIBUTE
240!
241 CALL mp_exchange2d (ng, tile, model, 4, &
242 & lbi, ubi, lbj, ubj, &
243 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
244 & si(:,:,1,isuevp), &
245 & si(:,:,2,isuevp), &
246 & si(:,:,1,isvevp), &
247 & si(:,:,2,isvevp))
248# endif
249!
250! Compute strain rates. Then, compute ice pressure (strength; icPice),
251! bulk viscosity (icBvis), and shear viscosity (icSvis).
252!
253 e2r = 1.0_r8/(ellip_sq(ng))
254!
255 DO j=jstr,jend
256 DO i=istr,iend
257 eps_xx(i,j)=(si(i+1,j,lieol,isuevp)- &
258 & si(i ,j,lieol,isuevp))*pm(i,j)
259 eps_yy(i,j)=(si(i,j+1,lieol,isvevp)- &
260 & si(i,j ,lieol,isvevp))*pn(i,j)
261 epx=0.25_r8*(si(i+1,j+1,lieol,isvevp)+ &
262 & si(i+1,j ,lieol,isvevp)- &
263 & si(i-1,j+1,lieol,isvevp)- &
264 & si(i-1,j ,lieol,isvevp))*pm(i,j)
265 epy=0.25_r8*(si(i+1,j+1,lieol,isuevp)+ &
266 & si(i ,j+1,lieol,isuevp)- &
267 & si(i+1,j-1,lieol,isuevp)- &
268 & si(i ,j-1,lieol,isuevp))*pn(i,j)
269 eps_xy(i,j)=0.5_r8*(epx+epy)
270!
271 eone=eps_xx(i,j)+eps_yy(i,j)
272 etwos=(eps_xx(i,j)-eps_yy(i,j))* &
273 (eps_xx(i,j)-eps_yy(i,j))+ &
274 & 4.0_r8*eps_xy(i,j)*eps_xy(i,j)
275!
276 delta=abs(eone**2+e2r*etwos)
277 delta=max(sqrt(delta),epso)
278# ifdef ICE_STRENGTH_QUAD
279 fi(i,j,icpice)=fi(i,j,icpgrd)* &
280 si(i,j,liold,ishice)*si(i,j,liold,ishice)* &
281 & exp(-astrength(ng)* &
282 & (1.0_r8-si(i,j,liold,isaice)))
283# else
284 fi(i,j,icpice)=pstar(ng)*si(i,j,liold,ishice)* &
285 & exp(-astrength(ng)* &
286 & (1.0_r8-si(i,j,liold,isaice)))
287# endif
288 fi(i,j,icbvis)=fi(i,j,icpice)/(2.0_r8*delta)
289 zmax=2.5e+8_r8*fi(i,j,icpice)
290 fi(i,j,icbvis)=min(fi(i,j,icbvis), zetamax(ng))
291 fi(i,j,icbvis)=max(fi(i,j,icbvis), zetamin(ng))
292 fi(i,j,icsvis)=e2r*fi(i,j,icbvis)
293 END DO
294 END DO
295!
296! Set lateral boundary conditions.
297!
298 CALL bc_r2d_tile (ng, tile, &
299 & lbi, ubi, lbj, ubj, &
300 & fi(:,:,icpice))
301
302 CALL bc_r2d_tile (ng, tile, &
303 & lbi, ubi, lbj, ubj, &
304 & fi(:,:,icbvis))
305
306 CALL bc_r2d_tile (ng, tile, &
307 & lbi, ubi, lbj, ubj, &
308 & fi(:,:,icsvis))
309
310# ifdef DISTRIBUTE
311 CALL mp_exchange2d (ng, tile, model, 3, &
312 & lbi, ubi, lbj, ubj, &
313 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
314 & fi(:,:,icpice), &
315 & fi(:,:,icbvis), &
316 & fi(:,:,icsvis))
317# endif
318
319# ifdef OUTFLOW_MASK
320!
321! Apply outflow mask.
322!
323 DO j=jstrt,jendt
324 DO i=istrt,iendt
325 fi(i,j,icpice)=fi(i,j,icpice)*mask_outflow(i,j)
326 fi(i,j,icbvis)=fi(i,j,icbvis)*mask_outflow(i,j)
327 fi(i,j,icsvis)=fi(i,j,icsvis)*mask_outflow(i,j)
328 END DO
329 END DO
330# endif
331!
332 RETURN
333 END SUBROUTINE ice_evp_tile
334#endif
335 END MODULE ice_evp_mod
bc_2d_mod
Definition bc_2d.F:2
bc_2d_mod::bc_r2d_tile
subroutine bc_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:44
ice_evp_mod
Definition ice_evp.F:2
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ice
Definition ice_mod.h:1
mod_ice::isvice
integer, parameter isvice
Definition ice_mod.h:147
mod_ice::pstar
real(r8), dimension(:), allocatable pstar
Definition ice_mod.h:236
mod_ice::ievp
integer, dimension(:), allocatable ievp
Definition ice_mod.h:212
mod_ice::icpice
integer, parameter icpice
Definition ice_mod.h:181
mod_ice::ellip_sq
real(r8), dimension(:), allocatable ellip_sq
Definition ice_mod.h:263
mod_ice::icsvis
integer, parameter icsvis
Definition ice_mod.h:184
mod_ice::icbvis
integer, parameter icbvis
Definition ice_mod.h:173
mod_ice::zetamax
real(r8), dimension(:), allocatable zetamax
Definition ice_mod.h:255
mod_ice::zetamin
real(r8), dimension(:), allocatable zetamin
Definition ice_mod.h:254
mod_ice::ice
type(t_ice), dimension(:), allocatable ice
Definition ice_mod.h:283
mod_ice::isuevp
integer, parameter isuevp
Definition ice_mod.h:150
mod_ice::icpgrd
integer, parameter icpgrd
Definition ice_mod.h:180
mod_ice::isaice
integer, parameter isaice
Definition ice_mod.h:137
mod_ice::isvevp
integer, parameter isvevp
Definition ice_mod.h:151
mod_ice::isuice
integer, parameter isuice
Definition ice_mod.h:146
mod_ice::astrength
real(r8), dimension(:), allocatable astrength
Definition ice_mod.h:235
mod_ice::ishice
integer, parameter ishice
Definition ice_mod.h:138
mod_param
Definition mod_param.F:2
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_scalars
Definition mod_scalars.F:2
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::iwest
integer, parameter iwest
Definition mod_scalars.F:1432
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_stepping
Definition mod_stepping.F:2
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3