ROMS
Loading...
Searching...
No Matches
ice_elastic.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE ice_elastic_mod
3
4#if defined ICE_MOMENTUM && defined ICE_EVP
5!
6!git $Id$
7!================================================== Hernan G. Arango ===
8! Copyright (c) 2002-2025 The ROMS Group Paul Budgell !
9! Licensed under a MIT/X style license Katherine Hedstrom !
10! See License_ROMS.md Scott M. Durski !
11!=======================================================================
12! !
13! This routine timesteps the ice momentum equations. !
14! !
15!=======================================================================
16!
17 USE mod_param
18# ifdef MICLM_NUDGING
19 USE mod_clima
20# endif
21# if defined ICE_SHOREFAST || defined ICE_KEEL
22 USE mod_coupling
23# endif
24# ifdef FASTICE_CLIMATOLOGY
25 USE mod_forces
26# endif
27 USE mod_grid
28 USE mod_ice
29 USE mod_scalars
30!
31 USE exchange_2d_mod, ONLY : exchange_u2d_tile, exchange_v2d_tile
32 USE ice_uibc_mod, ONLY : ice_uibc_tile
33 USE ice_vibc_mod, ONLY : ice_vibc_tile
34# ifdef DISTRIBUTE
35 USE mp_exchange_mod, ONLY : mp_exchange2d
36# endif
37!
38 implicit none
39!
40 PRIVATE
41 PUBLIC :: ice_elastic
42!
43 CONTAINS
44!
45!***********************************************************************
46 SUBROUTINE ice_elastic (ng, tile, model)
47!***********************************************************************
48!
49 USE mod_stepping
50!
51! Imported variable declarations.
52!
53 integer, intent(in) :: ng, tile, model
54!
55! Local variable declarations.
56!
57 character (len=*), parameter :: MyFile = &
58 & __FILE__
59!
60# include "tile.h"
61!
62# ifdef PROFILE
63 CALL wclock_on (ng, model, 42, __line__, myfile)
64# endif
65 CALL ice_elastic_tile (ng, tile, model, &
66 & lbi, ubi, lbj, ubj, &
67 & imins, imaxs, jmins, jmaxs, &
68 & liold(ng), liuol(ng), liunw(ng), &
69 & lieol(ng), lienw(ng), &
70# ifdef MASKING
71 & grid(ng) % rmask, &
72 & grid(ng) % umask, &
73 & grid(ng) % vmask, &
74# endif
75# ifdef WET_DRY
76 & grid(ng) % rmask_wet, &
77 & grid(ng) % umask_wet, &
78 & grid(ng) % vmask_wet, &
79# endif
80# ifdef ICESHELF
81 & grid(ng) % zice, &
82# endif
83 & grid(ng) % pm, &
84 & grid(ng) % pn, &
85 & grid(ng) % om_u, &
86 & grid(ng) % on_u, &
87 & grid(ng) % om_v, &
88 & grid(ng) % on_v, &
89 & grid(ng) % f, &
90# if defined ICE_SHOREFAST || defined ICE_KEEL
91 & grid(ng) % h, &
92 & coupling(ng) % Zt_avg1, &
93# endif
94# ifdef FASTICE_CLIMATOLOGY
95 & forces(ng) % fastice_clm, &
96# endif
97# ifdef MICLM_NUDGING
98 & clima(ng) % uiclm, &
99 & clima(ng) % viclm, &
100 & clima(ng) % MInudgcof, &
101# endif
102 & ice(ng) % Fi, &
103 & ice(ng) % Si)
104# ifdef PROFILE
105 CALL wclock_off (ng, model, 42, __line__, myfile)
106# endif
107!
108 RETURN
109 END SUBROUTINE ice_elastic
110!
111!***********************************************************************
112 SUBROUTINE ice_elastic_tile (ng, tile, model, &
113 & LBi, UBi, LBj, UBj, &
114 & IminS, ImaxS, JminS, JmaxS, &
115 & liold, liuol, liunw, lieol, lienw, &
116# ifdef MASKING
117 & rmask, umask, vmask, &
118# endif
119# ifdef WET_DRY
120 & rmask_wet, umask_wet, vmask_wet, &
121# endif
122# ifdef ICESHELF
123 & zice, &
124# endif
125 & pm, pn, om_u, on_u, om_v, on_v, &
126 & f, &
127# if defined ICE_SHOREFAST || defined ICE_KEEL
128 & h, &
129 & Zt_avg1, &
130# endif
131# ifdef FASTICE_CLIMATOLOGY
132 & fastice_clm, &
133# endif
134# ifdef MICLM_NUDGING
135 & uiclm, viclm, MInudgcof, &
136# endif
137 & Fi, Si)
138!***********************************************************************
139!
140! Imported variable declarations.
141!
142 integer, intent(in) :: ng, tile, model
143 integer, intent(in) :: LBi, UBi, LBj, UBj
144 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
145 integer, intent(in) :: liold, liuol, liunw, lieol, lienw
146!
147# ifdef ASSUMED_SHAPE
148# ifdef MASKING
149 real(r8), intent(in) :: rmask(LBi:,LBj:)
150 real(r8), intent(in) :: umask(LBi:,LBj:)
151 real(r8), intent(in) :: vmask(LBi:,LBj:)
152# endif
153# ifdef WET_DRY
154 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
155 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
156 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
157# endif
158# ifdef ICESHELF
159 real(r8), intent(in) :: zice(LBi:,LBj:)
160# endif
161 real(r8), intent(in) :: pm(LBi:,LBj:)
162 real(r8), intent(in) :: pn(LBi:,LBj:)
163 real(r8), intent(in) :: om_u(LBi:,LBj:)
164 real(r8), intent(in) :: on_u(LBi:,LBj:)
165 real(r8), intent(in) :: om_v(LBi:,LBj:)
166 real(r8), intent(in) :: on_v(LBi:,LBj:)
167 real(r8), intent(in) :: f(LBi:,LBj:)
168# if defined ICE_SHOREFAST || defined ICE_KEEL
169 real(r8), intent(in) :: h(LBi:,LBj:)
170 real(r8), intent(in) :: Zt_avg1(LBi:,LBj:)
171# endif
172# ifdef FASTICE_CLIMATOLOGY
173 real(r8), intent(in) :: fastice_clm(LBi:,LBj:)
174# endif
175# ifdef MICLM_NUDGING
176 real(r8), intent(in) :: uiclm(LBi:,LBj:)
177 real(r8), intent(in) :: viclm(LBi:,LBj:)
178 real(r8), intent(in) :: MInudgcof(LBi:,LBj:)
179# endif
180 real(r8), intent(in) :: Fi(LBi:,LBj:,:)
181 real(r8), intent(inout) :: Si(LBi:,LBj:,:,:)
182# else
183# ifdef MASKING
184 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
185 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
186 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
187# endif
188# ifdef WET_DRY
189 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
190 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
191 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
192# endif
193# ifdef ICESHELF
194 real(r8), intent(in) :: zice(LBi:UBi,LBj:UBj)
195# endif
196 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
197 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
198 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
199 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
200 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
201 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
202 real(r8), intent(in) :: f(LBi:UBi,LBj:UBj)
203# if defined ICE_SHOREFAST || defined ICE_KEEL
204 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
205 real(r8), intent(in) :: Zt_avg1(LBi:UBi,LBj:UBj)
206# endif
207# ifdef FASTICE_CLIMATOLOGY
208 real(r8), intent(in) :: fastice_clm(LBi:UBi,LBj:UBj)
209# endif
210# ifdef MICLM_NUDGING
211 real(r8), intent(in) :: uiclm(LBi:UBi,LBj:UBj)
212 real(r8), intent(in) :: viclm(LBi:UBi,LBj:UBj)
213 real(r8), intent(in) :: MInudgcof(LBi:UBi,LBj:UBj)
214# endif
215 real(r8), intent(in) :: Fi(LBi:UBi,LBj:UBj,nIceF)
216 real(r8), intent(inout) :: Si(LBi:UBi,LBj:UBj,2,nIceS)
217# endif
218!
219! Local variable definitions
220!
221 integer :: i, j
222!
223 real(r8) :: alfa, area, auf, avf, chux, chuy, clear, dsum, fakt
224# ifdef ICE_KEEL
225 real(r8) :: Cbu, Cbv, spd, thic, uvf, vuf
226# endif
227 real(r8) :: cff, f1, f2, s1, s2
228# if defined ICE_SHOREFAST || defined ICE_KEEL
229 real(r8) :: hu, hv
230# endif
231 real(r8) :: pmu, pnu, pmv, pnv
232 real(r8) :: cosstang, sinstang
233 real(r8) :: masu, mfu11, mfu21, mfu12, mfu22
234 real(r8) :: masv, mfv11, mfv21, mfv12, mfv22
235 real(r8) :: uforce, vforce
236
237# include "set_bounds.h"
238!
239!-----------------------------------------------------------------------
240! Timestep ice U-momentum equations.
241!-----------------------------------------------------------------------
242!
243 DO j=jstr,jend
244 DO i=istrp,iend
245 uforce=0.0_r8
246 area=om_u(i,j)*on_u(i,j)
247 chux=0.5_r8*(fi(i-1,j,iciomt)+ &
248 fi(i ,j,iciomt))
249 masu=0.5_r8*(si(i ,j,liold,ishice)+ &
250 & si(i-1,j,liold,ishice))
251 masu=max(masu, 0.1_r8)
252 masu=masu*icerho(ng)
253 auf=max(0.5_r8*(si(i-1,j,liold,isaice)+ &
254 & si(i ,j,liold,isaice)), 0.01_r8)
255 pmu=0.5_r8*(pm(i,j)+pm(i-1,j))
256 pnu=0.5_r8*(pn(i,j)+pn(i-1,j))
257!
258! Forces from ice rheology x-direction.
259!
260 s1=(si(i ,j,lienw,isisxx)- &
261 & si(i-1,j,lienw,isisxx))*pmu
262#ifdef MASKING
263 IF ((umask(i,j ).gt.0.0_r8).and. &
264 (umask(i,j+1).lt.1.0_r8)) THEN
265 f1=0.5_r8*(si(i ,j,lienw,isisxy)+ &
266 & si(i-1,j,lienw,isisxy))
267# ifdef WET_DRY
268 ELSE IF ((umask_wet(i,j ).gt.0.0_r8).and. &
269 & (umask_wet(i,j+1).lt.1.0_r8)) THEN
270 f1=0.5_r8*(si(i ,j,lienw,isisxy)+ &
271 & si(i-1,j,lienw,isisxy))
272# endif
273# ifdef ICESHELF
274 ELSE IF ((zice(i-1,j).eq.0.0_r8).and. &
275 & (zice(i ,j).eq.0.0_r8).and. &
276 & (zice(i-1,j+1)+zice(i,j+1)).ne.0.0_r8) THEN
277 f1=0.5_r8*(si(i ,j,lienw,isisxy)+ &
278 & si(i-1,j,lienw,isisxy))
279# endif
280 ELSE
281 f1=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
282 & si(i ,j+1,lienw,isisxy)+ &
283 & si(i-1,j+1,lienw,isisxy)+ &
284 & si(i-1,j ,lienw,isisxy))
285 END IF
286#else
287 f1=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
288 & si(i ,j+1,lienw,isisxy)+ &
289 & si(i-1,j+1,lienw,isisxy)+ &
290 & si(i-1,j ,lienw,isisxy))
291#endif
292
293#ifdef MASKING
294 IF ((umask(i,j ).gt.0.0_r8).and. &
295 & (umask(i,j-1).lt.1.0_r8)) THEN
296 f2=0.5_r8*(si(i ,j,lienw,isisxy)+ &
297 & si(i-1,j,lienw,isisxy))
298# ifdef WET_DRY
299 ELSE IF ((umask_wet(i,j ).gt.0.0_r8).and. &
300 & (umask_wet(i,j-1).lt.1.0_r8)) THEN
301 f2=0.5_r8*(si(i ,j,lienw,isisxy)+ &
302 & si(i-1,j,lienw,isisxy))
303# endif
304# ifdef ICESHELF
305 ELSE IF ((zice(i-1,j).eq.0.0_r8).and. &
306 & (zice(i ,j).eq.0.0_r8).and. &
307 & zice(i-1,j-1)+zice(i,j-1).ne.0.0_r8) THEN
308 f2=0.5_r8*(si(i ,j,lienw,isisxy)+ &
309 & si(i-1,j,lienw,isisxy))
310# endif
311 ELSE
312 f2=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
313 & si(i-1,j ,lienw,isisxy)+ &
314 & si(i-1,j-1,lienw,isisxy)+ &
315 & si(i ,j-1,lienw,isisxy))
316 END IF
317#else
318 f2=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
319 & si(i-1,j ,lienw,isisxy)+ &
320 & si(i-1,j-1,lienw,isisxy)+ &
321 & si(i ,j-1,lienw,isisxy))
322#endif
323 s2=(f1-f2)*pnu
324 uforce=(s1+s2)*area
325!
326! Add wind forcing.
327!
328 uforce=uforce+ &
329 & icerho(ng)*fi(i,j,icaius)*area
330!
331! Add pressure from tilting ocean surface.
332!
333 uforce=uforce- &
334 & g*masu*on_u(i,j)*(fi(i ,j,ichsse)- &
335 & fi(i-1,j,ichsse))
336!
337! Add Coriolis force.
338!
339 fakt=0.0_r8
340
341#ifdef MASKING
342# ifdef WET_DRY
343 dsum=vmask(i-1,j )*vmask_wet(i-1,j )+ &
344 & vmask(i ,j )*vmask_wet(i ,j )+ &
345 & vmask(i-1,j+1)*vmask_wet(i-1,j+1)+ &
346 & vmask(i ,j+1)*vmask_wet(i ,j+1)
347 IF (dsum.gt.0.0_r8) fakt=1.0_r8/dsum
348
349 mfv11=0.5_r8*(si(i-1,j-1,liold,ishice)*f(i-1,j-1)+ &
350 & si(i-1,j ,liold,ishice)*f(i-1,j ))* &
351 & si(i-1,j,lieol,isvevp)* &
352 & vmask(i-1,j)*vmask_wet(i-1,j)
353 mfv21=0.5_r8*(si(i,j-1,liold,ishice)*f(i,j-1)+ &
354 & si(i,j ,liold,ishice)*f(i,j ))* &
355 & si(i,j,lieol,isvevp)* &
356 & vmask(i,j)*vmask_wet(i,j)
357 mfv12=0.5_r8*(si(i-1,j ,liold,ishice)*f(i-1,j )+ &
358 & si(i-1,j+1,liold,ishice)*f(i-1,j+1))* &
359 & si(i-1,j+1,lieol,isvevp)* &
360 & vmask(i-1,j+1)*vmask_wet(i-1,j+1)
361 mfv22=0.5_r8*(si(i,j ,liold,ishice)*f(i,j )+ &
362 & si(i,j+1,liold,ishice)*f(i,j+1))* &
363 & si(i,j+1,lieol,isvevp)* &
364 & vmask(i,j+1)*vmask_wet(i,j+1)
365# else
366 dsum=vmask(i-1,j )+ &
367 & vmask(i ,j )+ &
368 & vmask(i-1,j+1)+ &
369 & vmask(i ,j+1)
370 IF (dsum.gt.0.0_r8) fakt=1.0_r8/dsum
371
372 mfv11=0.5_r8*(si(i-1,j-1,liold,ishice)*f(i-1,j-1)+ &
373 & si(i-1,j ,liold,ishice)*f(i-1,j ))* &
374 & si(i-1,j,lieol,isvevp)*vmask(i-1,j)
375 mfv21=0.5_r8*(si(i,j-1,liold,ishice)*f(i,j-1)+ &
376 & si(i,j ,liold,ishice)*f(i,j ))* &
377 & si(i,j,lieol,isvevp)*vmask(i,j)
378 mfv12=0.5_r8*(si(i-1,j ,liold,ishice)*f(i-1,j )+ &
379 & si(i-1,j+1,liold,ishice)*f(i-1,j+1))* &
380 & si(i-1,j+1,lieol,isvevp)*vmask(i-1,j+1)
381 mfv22=0.5_r8*(si(i,j ,liold,ishice)*f(i,j )+ &
382 & si(i,j+1,liold,ishice)*f(i,j+1))* &
383 & si(i,j+1,lieol,isvevp)*vmask(i,j+1)
384# endif
385#else
386 fakt=0.25_r8
387 mfv11=0.5_r8*(si(i-1,j-1,liold,ishice)*f(i-1,j-1)+ &
388 & si(i-1,j ,liold,ishice)*f(i-1,j ))* &
389 & si(i-1,j,lieol,isvevp)
390 mfv21=0.5_r8*(si(i,j-1,liold,ishice)*f(i,j-1)+ &
391 & si(i,j ,liold,ishice)*f(i,j ))* &
392 & si(i,j,lieol,isvevp)
393 mfv12=0.5_r8*(si(i-1,j ,liold,ishice)*f(i-1,j )+ &
394 & si(i-1,j+1,liold,ishice)*f(i-1,j+1))* &
395 & si(i-1,j+1,lieol,isvevp)
396 mfv22=0.5_r8*(si(i,j ,liold,ishice)*f(i,j )+ &
397 & si(i,j+1,liold,ishice)*f(i,j+1))* &
398 & si(i,j+1,lieol,isvevp)
399#endif
400!
401 uforce=uforce+ &
402 & fakt*icerho(ng)*area* &
403 & (mfv11+mfv21+mfv12+mfv22)
404!
405! Add stress from ocean current.
406!
407 uforce=uforce/area + &
408 & auf*chux*rho0*fi(i,j,icuavg)
409!
410! Compute basal stress as bottom drag to the deepest ice keels.
411# ifdef ICE_KEEL
412! Use the ice keel basal stress parameterization of Lemieux (2005).
413!
414 hu=0.5_r8*(h(i-1,j)+min(zt_avg1(i-1,j), 0.0_r8) + &
415 & h(i ,j)+min(zt_avg1(i ,j), 0.0_r8))
416 thic=0.5_r8*(si(i ,j,liold,ishice)+ &
417 & si(i-1,j,liold,ishice))/auf
418 vuf=0.25*(si(i ,j ,lieol,isvevp)+ &
419 & si(i ,j+1,lieol,isvevp)+ &
420 & si(i-1,j ,lieol,isvevp)+ &
421 & si(i-1,j+1,lieol,isvevp))
422 spd=sqrt(si(i,j,lieol,isuevp)* &
423 & si(i,j,lieol,isuevp)+vuf*vuf)
424!
425! Determine if the ice keel depth is greater than the water depth
426! (Lemieux k_1 parameter default value is 8.0)
427!
428 clear=max(thic-hu/8.0_r8, 0.0_r8)
429!
430! Compute basal drag coefficient, C_bu. Use Lemieux k_2=20.0.
431!
432 cbu=20.0_r8*clear*exp(-20.0_r8*(1.0_r8-auf))/(spd+5.0e-5_r8)
433!
434! For an implicit formulation, the basal stress just appears in the
435! RHS denominator.
436!
437 alfa=masu+dtevp(ng)*auf*rho0*chux + &
438 & dtevp(ng)*icerho(ng)*cbu
439# else
440 alfa=masu+dtevp(ng)*auf*rho0*chux
441# endif
442!
443! Solving the ice U-momentum equation.
444!
445 si(i,j,lienw,isuevp)=(masu*si(i,j,lieol,isuevp)+ &
446 & dtevp(ng)*uforce)/alfa
447!
448# ifdef MICLM_NUDGING
449 cff = 0.5*(minudgcof(i,j)+minudgcof(i-1,j))
450 si(i,j,lienw,isuevp)=si(i,j,lienw,isuevp)+ &
451 & dtevp(ng)*cff*(uiclm(i,j)- &
452 & si(i,j,lienw,isuevp))
453# endif
454# ifdef ICE_SHOREFAST
455 hu=0.5_r8*(h(i-1,j)+min(zt_avg1(i-1,j), 0.0_r8)+ &
456 & h(i ,j)+min(zt_avg1(i ,j), 0.0_r8))
457 masu=0.5_r8*(si(i ,j,liold,ishice)+ &
458 & si(i-1,j,liold,ishice))
459 clear=hu-0.9_r8*masu
460 clear=max(clear, 0.0_r8)
461 IF (clear.lt.5.0_r8) &
462 & si(i,j,lienw,isuevp)=max(clear-1.0_r8, 0.0_r8)/ &
463 & 4.0_r8*si(i,j,lienw,isuevp)
464# endif
465# ifdef FASTICE_CLIMATOLOGY
466 si(i,j,lienw,isuevp)=0.5_r8*(fastice_clm(i-1,j)+ &
467 & fastice_clm(i ,j))* &
468 & si(i,j,lienw,isuevp)
469# endif
470# ifdef MASKING
471 si(i,j,lienw,isuevp)=si(i,j,lienw,isuevp)*umask(i,j)
472# endif
473# ifdef ICESHELF
474 IF(zice(i-1,j)+zice(i,j).ne.0.0_r8) THEN
475 si(i,j,lienw,isuevp)=0.0_r8
476 END IF
477# endif
478 END DO
479 END DO
480!
481!-----------------------------------------------------------------------
482! Timestep ice V-momentum equations.
483!-----------------------------------------------------------------------
484!
485 DO j=jstrp,jend
486 DO i=istr,iend
487 vforce=0.0_r8
488 area=om_v(i,j)*on_v(i,j)
489 masv=0.5_r8*(si(i,j ,liold,ishice)+ &
490 & si(i,j-1,liold,ishice))
491 masv=max(masv, 0.1_r8)
492 masv=masv*icerho(ng)
493 avf=max(0.5_r8*(si(i,j-1,liold,isaice)+ &
494 & si(i,j ,liold,isaice)), 0.01_r8)
495 chuy=0.5_r8*(fi(i,j-1,iciomt)+ &
496 & fi(i,j ,iciomt))
497 pmv=0.5_r8*(pm(i,j)+pm(i,j-1))
498 pnv=0.5_r8*(pn(i,j)+pn(i,j-1))
499!
500! Forces from ice rheology y-direction.
501!
502 s1=(si(i,j ,lienw,isisyy)- &
503 & si(i,j-1,lienw,isisyy))*pnv
504#ifdef MASKING
505 IF ((vmask(i ,j).gt.0.0_r8).and. &
506 & (vmask(i+1,j).lt.1.0_r8)) THEN
507 f1=0.5_r8*(si(i,j ,lienw,isisxy)+ &
508 & si(i,j-1,lienw,isisxy))
509# ifdef WET_DRY
510 ELSE IF ((vmask_wet(i ,j).gt.0.0_r8).and. &
511 & (vmask_wet(i+1,j).lt.1.0_r8)) THEN
512 f1=0.5_r8*(si(i,j ,lienw,isisxy)+ &
513 & si(i,j-1,lienw,isisxy))
514# endif
515# ifdef ICESHELF
516 ELSE IF ((zice(i,j-1).eq.0.0_r8).and. &
517 & (zice(i,j ).eq.0.0_r8).and. &
518 & (zice(i+1,j-1)+zice(i+1,j)).ne.0.0_r8) THEN
519 f1=0.5_r8*(si(i,j ,lienw,isisxy)+ &
520 & si(i,j-1,lienw,isisxy))
521# endif
522 ELSE
523 f1=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
524 & si(i+1,j ,lienw,isisxy)+ &
525 & si(i+1,j-1,lienw,isisxy)+ &
526 & si(i ,j-1,lienw,isisxy))
527 END IF
528#else
529 f1=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
530 & si(i+1,j ,lienw,isisxy)+ &
531 & si(i+1,j-1,lienw,isisxy)+ &
532 & si(i ,j-1,lienw,isisxy))
533#endif
534
535#ifdef MASKING
536 IF ((vmask(i ,j).gt.0.0_r8).and. &
537 & (vmask(i-1,j).lt.1.0_r8)) THEN
538 f2=0.5_r8*(si(i,j ,lienw,isisxy)+ &
539 & si(i,j-1,lienw,isisxy))
540# ifdef WET_DRY
541 ELSE IF ((vmask_wet(i ,j).gt.0.0_r8).and. &
542 & (vmask_wet(i-1,j).lt.1.0_r8)) THEN
543 f2=0.5_r8*(si(i,j ,lienw,isisxy)+ &
544 & si(i,j-1,lienw,isisxy))
545# endif
546# ifdef ICESHELF
547 ELSE IF ((zice(i,j-1).eq.0.0_r8).and. &
548 & (zice(i,j ).eq.0.0_r8).and. &
549 & (zice(i-1,j-1)+zice(i-1,j)).ne.0.0_r8) THEN
550 f2=0.5_r8*(si(i,j ,lienw,isisxy)+ &
551 & si(i,j-1,lienw,isisxy))
552# endif
553 ELSE
554 f2=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
555 & si(i-1,j ,lienw,isisxy)+ &
556 & si(i-1,j-1,lienw,isisxy)+ &
557 & si(i ,j-1,lienw,isisxy))
558 END IF
559#else
560 f2=0.25_r8*(si(i ,j ,lienw,isisxy)+ &
561 & si(i-1,j ,lienw,isisxy)+ &
562 & si(i-1,j-1,lienw,isisxy)+ &
563 & si(i ,j-1,lienw,isisxy))
564#endif
565 s2=(f1-f2)*pmv
566 vforce=(s1+s2)*area
567!
568! Add wind forcing.
569!
570 vforce=vforce+ &
571 & icerho(ng)*fi(i,j,icaivs)*area
572!
573! Add pressure from tilting ocean surface.
574!
575 vforce=vforce- &
576 & g*masv*om_v(i,j)*(fi(i,j ,ichsse)- &
577 & fi(i,j-1,ichsse))
578!
579! Add Coriolis force.
580!
581 fakt=0.0_r8
582
583#ifdef MASKING
584# ifdef WET_DRY
585 dsum=umask(i ,j-1)*umask_wet(i ,j-1)+ &
586 & umask(i+1,j-1)*umask_wet(i+1,j-1)+ &
587 & umask(i ,j )*umask_wet(i ,j )+ &
588 & umask(i+1,j )*umask_wet(i ,j )
589 IF (dsum.gt.0.0_r8) fakt=1.0_r8/dsum
590
591 mfu11=0.5_r8*(si(i-1,j-1,liold,ishice)*f(i-1,j-1)+ &
592 & si(i ,j-1,liold,ishice)*f(i ,j-1))* &
593 & si(i,j-1,lieol,isuevp)* &
594 & umask(i,j-1)*umask_wet(i,j-1)
595 mfu21=0.5_r8*(si(i ,j-1,liold,ishice)*f(i ,j-1)+ &
596 & si(i+1,j-1,liold,ishice)*f(i+1,j-1))* &
597 & si(i+1,j-1,lieol,isuevp)* &
598 & umask(i+1,j-1)*umask_wet(i+1,j-1)
599 mfu12=0.5_r8*(si(i-1,j,liold,ishice)*f(i-1,j)+ &
600 & si(i ,j,liold,ishice)*f(i ,j))* &
601 & si(i,j,lieol,isuevp)* &
602 & umask(i,j)*umask_wet(i,j)
603 mfu22=0.5_r8*(si(i ,j,liold,ishice)*f(i ,j)+ &
604 & si(i+1,j,liold,ishice)*f(i+1,j))* &
605 & si(i+1,j,lieol,isuevp)* &
606 & umask(i+1,j)*umask_wet(i+1,j)
607# else
608 dsum=umask(i ,j-1)+ &
609 & umask(i+1,j-1)+ &
610 & umask(i ,j )+ &
611 & umask(i+1,j )
612 IF (dsum.gt.0.0_r8) fakt=1.0_r8/dsum
613
614 mfu11=0.5_r8*(si(i-1,j-1,liold,ishice)*f(i-1,j-1)+ &
615 & si(i ,j-1,liold,ishice)*f(i ,j-1))* &
616 & si(i,j-1,lieol,isuevp)*umask(i,j-1)
617 mfu21=0.5_r8*(si(i ,j-1,liold,ishice)*f(i ,j-1)+ &
618 & si(i+1,j-1,liold,ishice)*f(i+1,j-1))* &
619 & si(i+1,j-1,lieol,isuevp)*umask(i+1,j-1)
620 mfu12=0.5_r8*(si(i-1,j,liold,ishice)*f(i-1,j)+ &
621 & si(i ,j,liold,ishice)*f(i ,j))* &
622 & si(i,j,lieol,isuevp)*umask(i,j)
623 mfu22=0.5_r8*(si(i ,j,liold,ishice)*f(i ,j)+ &
624 & si(i+1,j,liold,ishice)*f(i+1,j))* &
625 & si(i+1,j,lieol,isuevp)*umask(i+1,j)
626# endif
627#else
628 fakt=0.25_r8
629 mfu11=0.5_r8*(si(i-1,j-1,liold,ishice)*f(i-1,j-1)+ &
630 & si(i ,j-1,liold,ishice)*f(i ,j-1))* &
631 & si(i,j-1,lieol,isuevp)
632 mfu21=0.5_r8*(si(i ,j-1,liold,ishice)*f(i ,j-1)+ &
633 & si(i+1,j-1,liold,ishice)*f(i+1,j-1))* &
634 & si(i+1,j-1,lieol,isuevp)
635 mfu12=0.5_r8*(si(i-1,j,liold,ishice)*f(i-1,j)+ &
636 & si(i ,j,liold,ishice)*f(i ,j))* &
637 & si(i,j,lieol,isuevp)
638 mfu22=0.5_r8*(si(i ,j,liold,ishice)*f(i ,j)+ &
639 & si(i+1,j,liold,ishice)*f(i+1,j))* &
640 & si(i+1,j,lieol,isuevp)
641#endif
642!
643 vforce=vforce- &
644 & fakt*icerho(ng)*area* &
645 & (mfu11+mfu21+mfu12+mfu22)
646!
647! Add stress from ocean current.
648!
649 vforce=vforce/area + &
650 & avf*chuy*rho0*fi(i,j,icvavg)
651!
652! Compute basal stress as bottom drag to the deepest ice keels.
653# ifdef ICE_KEEL
654! Use the ice keel basal stress parameterization of Lemieux (2005).
655!
656 hv=0.5_r8*(h(i,j-1)+min(zt_avg1(i,j-1), 0.0_r8)+ &
657 & h(i,j )+min(zt_avg1(i,j ), 0.0_r8))
658 thic=0.5_r8*(si(i,j ,liold,ishice)+ &
659 & si(i,j-1,liold,ishice))/avf
660 uvf=0.25*(si(i ,j ,lieol,isuevp)+ &
661 & si(i ,j-1,lieol,isuevp)+ &
662 & si(i+1,j ,lieol,isuevp)+
663 & si(i+1,j-1,lieol,isuevp))
664 spd=sqrt(si(i,j,lieol,isvevp)* &
665 & si(i,j,lieol,isvevp)+uvf*uvf)
666!
667! Determine if the ice keel depth is greater than the water depth
668! (Lemieux k_1 parameter default value is 8.0)
669!
670 clear=max(thic-hv/8.0_r8, 0.0_r8)
671!
672! Compute basal drag coefficient, C_bu. Use Lemieux k_2=20.0.
673!
674 cbv=20.0_r8*clear*exp(-20.0_r8*(1.0_r8-avf))/(spd+5.0e-5_r8)
675!
676! For an implicit formulation, the basal stress just appears in the
677! RHS denominator.
678!
679 alfa=masv+dtevp(ng)*avf*rho0*chuy+ &
680 & dtevp(ng)*icerho(ng)*cbv
681# else
682 alfa=masv+dtevp(ng)*avf*rho0*chuy
683# endif
684!
685! Solving the ice V-momentum equation.
686!
687 si(i,j,lienw,isvevp)=(masv*si(i,j,lieol,isvevp)+ &
688 & dtevp(ng)*vforce)/alfa
689!
690# ifdef MICLM_NUDGING
691 cff=0.5*(minudgcof(i,j)+minudgcof(i,j-1))
692 si(i,j,lienw,isvevp)=si(i,j,lienw,isvevp)+ &
693 & dtevp(ng)*cff*(viclm(i,j)- &
694 & si(i,j,lienw,isvevp))
695# endif
696# ifdef ICE_SHOREFAST
697 hv=0.5_r8*(h(i,j-1)+min(zt_avg1(i,j-1), 0.0_r8) + &
698 & h(i,j )+min(zt_avg1(i,j ), 0.0_r8))
699 masv=0.5_r8*(si(i,j ,liold,ishice)+ &
700 & si(i,j-1,liold,ishice))
701 clear=hv-0.9_r8*masv
702 clear=max(clear, 0.0_r8)
703 IF (clear.lt.5.0_r8) &
704 & si(i,j,lienw,isvevp)=max(clear-1.0_r8, 0.0_r8)/ &
705 & 4.0_r8*si(i,j,lienw,isvevp)
706# endif
707# ifdef FASTICE_CLIMATOLOGY
708 si(i,j,lienw,isvevp)=0.5*(fastice_clm(i,j-1)+ &
709 & fastice_clm(i,j ))*
710 & si(i,j,lienw,isvevp)
711# endif
712# ifdef MASKING
713 si(i,j,lienw,isvevp) = si(i,j,lienw,isvevp)*vmask(i,j)
714# endif
715# ifdef ICESHELF
716 IF (zice(i,j-1)+zice(i,j).ne.0.0_r8) THEN
717 si(i,j,lienw,isvevp)=0.0_r8
718 END IF
719# endif
720 END DO
721 END DO
722!
723!-----------------------------------------------------------------------
724! Apply lateral boundary conditions.
725!-----------------------------------------------------------------------
726!
727 CALL ice_uibc_tile (ng, tile, model, &
728 & lbi, ubi, lbj, ubj, &
729 & imins, imaxs, jmins, jmaxs, &
730 & lieol, lienw, &
731 & si(:,:,:,isuevp))
732
733 CALL ice_vibc_tile (ng, tile, model, &
734 & lbi, ubi, lbj, ubj, &
735 & imins, imaxs, jmins, jmaxs, &
736 & lieol, lienw, &
737 & si(:,:,:,isvevp))
738!
739 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
740 CALL exchange_u2d_tile (ng, tile, &
741 & lbi, ubi, lbj, ubj, &
742 & si(:,:,lienw,isuevp))
743 CALL exchange_v2d_tile (ng, tile, &
744 & lbi, ubi, lbj, ubj, &
745 & si(:,:,lienw,isvevp))
746 END IF
747
748# ifdef DISTRIBUTE
749!
750 CALL mp_exchange2d (ng, tile, model, 2, &
751 & lbi, ubi, lbj, ubj, &
752 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
753 & si(:,:,lienw,isuevp), &
754 & si(:,:,lienw,isvevp))
755# endif
756
757!
758!-----------------------------------------------------------------------
759! Update ice velocity
760!-----------------------------------------------------------------------
761!
762 IF(ievp(ng).eq.nevp(ng)) THEN
763 DO j=jstr,jend
764 DO i=istrp,iend
765 si(i,j,liunw,isuice)=si(i,j,lienw,isuevp)
766 END DO
767 END DO
768 DO j=jstrp,jend
769 DO i=istr,iend
770 si(i,j,liunw,isvice)=si(i,j,lienw,isvevp)
771 END DO
772 END DO
773!
774 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
775 IF (domain(ng)%Western_Edge(tile)) THEN
776 DO j=jstr,jend
777 si(istr,j,liunw,isuice)=si(istr,j,lienw,isuevp)
778 END DO
779 DO j=jstrp,jend
780 si(istr-1,j,liunw,isvice)=si(istr-1,j,lienw,isvevp)
781 END DO
782 END IF
783 END IF
784!
785 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
786 IF (domain(ng)%Eastern_Edge(tile)) THEN
787 DO j=jstr,jend
788 si(iend+1,j,liunw,isuice)=si(iend+1,j,lienw,isuevp)
789 END DO
790 DO j=jstrp,jend
791 si(iend+1,j,liunw,isvice)=si(iend+1,j,lienw,isvevp)
792 END DO
793 END IF
794 END IF
795!
796 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
797 IF (domain(ng)%Southern_Edge(tile)) THEN
798 DO i=istrp,iend
799 si(i,jstr-1,liunw,isuice)=si(i,jstr-1,lienw,isuevp)
800 END DO
801 DO i=istr,iend
802 si(i,jstr,liunw,isvice)=si(i,jstr,lienw,isvevp)
803 END DO
804 END IF
805 END IF
806!
807 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
808 IF (domain(ng)%Northern_Edge(tile)) THEN
809 DO i=istrp,iend
810 si(i,jend+1,liunw,isuice)=si(i,jend+1,lienw,isuevp)
811 END DO
812 DO i=istr,iend
813 si(i,jend+1,liunw,isvice)=si(i,jend+1,lienw,isvevp)
814 END DO
815 END IF
816 END IF
817!
818 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
819 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
820 IF (domain(ng)%SouthWest_Corner(tile)) THEN
821 si(istr ,jstr-1,liunw,isuice)=si(istr ,jstr-1,lienw, &
822 & isuevp)
823 si(istr-1,jstr ,liunw,isvice)=si(istr-1,jstr ,lienw, &
824 & isvevp)
825 END IF
826 END IF
827!
828 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
829 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
830 IF (domain(ng)%SouthEast_Corner(tile)) THEN
831 si(iend+1,jstr-1,liunw,isuice)=si(iend+1,jstr-1,lienw, &
832 & isuevp)
833 si(iend+1,jstr ,liunw,isvice)=si(iend+1,jstr ,lienw, &
834 & isvevp)
835 END IF
836 END IF
837!
838 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
839 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
840 IF (domain(ng)%NorthWest_Corner(tile)) THEN
841 si(istr ,jend+1,liunw,isuice)=si(istr ,jend+1,lienw, &
842 & isuevp)
843 si(istr-1,jend+1,liunw,isvice)=si(istr-1,jend+1,lienw, &
844 & isvevp)
845 END IF
846 END IF
847!
848 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
849 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
850 IF (domain(ng)%NorthEast_Corner(tile)) THEN
851 si(iend+1,jend+1,liunw,isuice)=si(iend+1,jend+1,lienw, &
852 & isuevp)
853 si(iend+1,jend+1,liunw,isvice)=si(iend+1,jend+1,lienw, &
854 & isvevp)
855 END IF
856 END IF
857!
858! Exchange boundary data.
859!
860 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
861 CALL exchange_u2d_tile (ng, tile, &
862 & lbi, ubi, lbj, ubj, &
863 & si(:,:,liunw,isuice))
864
865 CALL exchange_v2d_tile (ng, tile, &
866 & lbi, ubi, lbj, ubj, &
867 & si(:,:,liunw,isvice))
868 END IF
869
870# ifdef DISTRIBUTE
871!
872 CALL mp_exchange2d (ng, tile, inlm, 2, &
873 & lbi, ubi, lbj, ubj, &
874 & nghostpoints, ewperiodic(ng), nsperiodic(ng), &
875 & si(:,:,liunw,isuice), &
876 & si(:,:,liunw,isvice))
877# endif
878 ENDIF
879!
880 RETURN
881 END SUBROUTINE ice_elastic_tile
882#endif
883 END MODULE ice_elastic_mod
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_u2d_tile
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:440
exchange_2d_mod::exchange_v2d_tile
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:627
ice_elastic_mod
Definition ice_elastic.F:2
ice_uibc_mod
Definition ice_uibc.F:2
ice_vibc_mod
Definition ice_vibc.F:2
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
mod_coupling
Definition mod_coupling.F:2
mod_coupling::coupling
type(t_coupling), dimension(:), allocatable coupling
Definition mod_coupling.F:102
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
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::ievp
integer, dimension(:), allocatable ievp
Definition ice_mod.h:212
mod_ice::dtevp
integer, dimension(:), allocatable dtevp
Definition ice_mod.h:218
mod_ice::ichsse
integer, parameter ichsse
Definition ice_mod.h:174
mod_ice::icuavg
integer, parameter icuavg
Definition ice_mod.h:187
mod_ice::isisxx
integer, parameter isisxx
Definition ice_mod.h:142
mod_ice::icerho
real(r8), dimension(:), allocatable icerho
Definition ice_mod.h:223
mod_ice::ice
type(t_ice), dimension(:), allocatable ice
Definition ice_mod.h:283
mod_ice::nevp
integer, dimension(:), allocatable nevp
Definition ice_mod.h:213
mod_ice::icvavg
integer, parameter icvavg
Definition ice_mod.h:188
mod_ice::iciomt
integer, parameter iciomt
Definition ice_mod.h:177
mod_ice::isisxy
integer, parameter isisxy
Definition ice_mod.h:143
mod_ice::isuevp
integer, parameter isuevp
Definition ice_mod.h:150
mod_ice::icaius
integer, parameter icaius
Definition ice_mod.h:171
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::ishice
integer, parameter ishice
Definition ice_mod.h:138
mod_ice::icaivs
integer, parameter icaivs
Definition ice_mod.h:172
mod_ice::isisyy
integer, parameter isisyy
Definition ice_mod.h:144
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
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::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
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