ice_evp_sig.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE ice_evp_sig_mod
 
#if defined ICE_MOMENTUM && defined ICE_EVP
!
!git $Id$
!=======================================================================
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license           W. Paul Budgell    !
!    See License_ROMS.md                            Katherine Hedstrom !
!================================================== Hernan G. Arango ===
!                                                                      !
!  This routine timesteps the EVP stresses term.                       !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_boundary
      USE mod_grid
      USE mod_ice
      USE mod_scalars
!
      USE exchange_2d_mod, ONLY : exchange_r2d_tile
      USE ice_bc2d_mod,    ONLY : ice_bc2d_tile
# ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
!
      implicit none
!
      PRIVATE
      PUBLIC  :: ice_evp_sig
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE ice_evp_sig (ng, tile, model)
!***********************************************************************
!
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
!
!  Local variable declarations.
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, model, 42, __line__, myfile)
# endif
      CALL ice_evp_sig_tile (ng, tile, model,                           &
     &                       lbi, ubi, lbj, ubj,                        &
     &                       imins, imaxs, jmins, jmaxs,                &
     &                       liold(ng), lieol(ng), lienw(ng),           &
     &                       grid(ng) % pm,                             &
     &                       grid(ng) % pn,                             &
     &                       ice(ng) % Fi,                              &
     &                       ice(ng) % Si)
# ifdef PROFILE
      CALL wclock_off (ng, model, 42, __line__, myfile)
# endif
!
      RETURN
      END SUBROUTINE ice_evp_sig
!
!***********************************************************************
      SUBROUTINE ice_evp_sig_tile (ng, tile, model,                     &
     &                             LBi, UBi, LBj, UBj,                  &
     &                             IminS, ImaxS, JminS, JmaxS,          &
     &                             liold, lieol, lienw,                 &
     &                             pm, pn,                              &
     &                             Fi, Si)
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: liold, lieol, lienw
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: pm(LBi:,LBj:)
      real(r8), intent(in) :: pn(LBi:,LBj:)
      real(r8), intent(in) :: Fi(LBi:,LBj:,:)
      real(r8), intent(inout) :: Si(LBi:,LBj:,:,:)
# else
      real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Fi(LBi:UBi,LBj:UBj,nIceF)
      real(r8), intent(inout) :: Si(LBi:UBi,LBj:UBj,2,nIceS)
# endif
!
!  Local variable definitions.
!
      integer :: i, j
!
      real(r8) :: alfa, beta, cff, f1, f2, f3, f4, gamma
      real(r8) :: E, E0, ee, ees, ep, epx, epy
!
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: eps_xx
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: eps_xy
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: eps_yy
 
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Timestep internal ice stress tensor term.
!-----------------------------------------------------------------------
!
!  Initial value for Young modulus, E.  It has values between 0 and 1,
!  and makes the elastic term small compared to other terms.
!
      e0=0.25_r8
!
!  Comput stress tensor.
!
      DO j=jstr,jend
        DO i=istr,iend
          eps_xx(i,j)=(si(i+1,j,lieol,isuevp)-                          &
     &                 si(i  ,j,lieol,isuevp))*pm(i,j)
          eps_yy(i,j)=(si(i,j+1,lieol,isvevp)-                          &
     &                 si(i,j  ,lieol,isvevp))*pn(i,j)
          epx=0.25_r8*(si(i+1,j+1,lieol,isvevp)+                        &
     &                 si(i+1,j  ,lieol,isvevp)-                        &
     &                 si(i-1,j+1,lieol,isvevp)-                        &
     &                 si(i-1,j  ,lieol,isvevp))*pm(i,j)
          epy=0.25_r8*(si(i+1,j+1,lieol,isuevp)+                        &
     &                 si(i  ,j+1,lieol,isuevp)-                        &
     &                 si(i+1,j-1,lieol,isuevp)-                        &
     &                 si(i  ,j-1,lieol,isuevp))*pn(i,j)
          eps_xy(i,j)=0.5_r8*(epx+epy)
        END DO
      END DO
!
!  Timestep internal ice stress tensor term. A split explicit scheme is
!  used with shorter steps (dts) to resolve the elastic wave velocity.
!
      DO j=jstr,jend
        DO i=istr,iend
          IF (si(i,j,liold,ishice).gt.0.01_r8) THEN
            e  =2.0_r8*e0*icerho(ng)*si(i,j,liold,ishice)/              &
     &          (pm(i,j)*dtevp(ng))**2
            ep =e*fi(i,j,icpice)/(4.0_r8*fi(i,j,icbvis)+1.0e-8_r8)
            ee =e/(2.0_r8*fi(i,j,icsvis)+1.0e-8_r8)
            ees=e*(fi(i,j,icsvis)-fi(i,j,icbvis))/                      &
     &          (4.0_r8*fi(i,j,icsvis)*fi(i,j,icbvis)+1.0e-8_r8)
!
            cff  =1.0_r8/dtevp(ng)
            alfa =cff+ee+ees
            beta =ees
            gamma=1.0_r8/dtevp(ng)+ee
            f1=e*eps_xx(i,j)+cff*si(i,j,lieol,isisxx)-ep
            f2=e*eps_yy(i,j)+cff*si(i,j,lieol,isisyy)-ep
            f3=e*eps_xy(i,j)+cff*si(i,j,lieol,isisxy)
            f4=1.0_r8/(alfa**2-beta**2)
            si(i,j,lienw,isisxx)=f4*(alfa*f1-beta*f2)
            si(i,j,lienw,isisxy)=f3/gamma
            si(i,j,lienw,isisyy)=f4*(alfa*f2-beta*f1)
          ELSE
            si(i,j,lienw,isisxx)=2.0_r8*fi(i,j,icsvis)*eps_xx(i,j)+     &
     &                           (fi(i,j,icbvis)-fi(i,j,icsvis))*       &
     &                           (eps_xx(i,j)+eps_yy(i,j))-             &
     &                           0.5_r8*fi(i,j,icpice)
            si(i,j,lienw,isisyy)=2.0_r8*fi(i,j,icsvis)*eps_yy(i,j)+     &
     &                           (fi(i,j,icbvis)-fi(i,j,icsvis))*       &
     &                           (eps_xx(i,j)+eps_yy(i,j))-             &
     &                           0.5_r8*fi(i,j,icpice)
            si(i,j,lienw,isisxy)=2.0_r8*fi(i,j,icsvis)*eps_xy(i,j)
          END IF
        END DO
      END DO
!
!  Apply lateral boundary conditions.
!
      CALL ice_bc2d_tile (ng, tile, model, isisxx,                      &
     &                    lbi, ubi, lbj, ubj,                           &
     &                    imins, imaxs, jmins, jmaxs,                   &
     &                    lieol, lienw,                                 &
     &                    si(:,:,:,isuevp),                             &
     &                    si(:,:,:,isvevp),                             &
     &                    si(:,:,:,isisxx),                             &
     &                    lbc(:,ibice(isisxx),ng))
 
      CALL ice_bc2d_tile (ng, tile, model, isisxy,                      &
     &                    lbi, ubi, lbj, ubj,                           &
     &                    imins, imaxs, jmins, jmaxs,                   &
     &                    lieol, lienw,                                 &
     &                    si(:,:,:,isuevp),                             &
     &                    si(:,:,:,isvevp),                             &
     &                    si(:,:,:,isisxy),                             &
     &                    lbc(:,ibice(isisxy),ng))
 
      CALL ice_bc2d_tile (ng, tile, model, isisyy,                      &
     &                    lbi, ubi, lbj, ubj,                           &
     &                    imins, imaxs, jmins, jmaxs,                   &
     &                    lieol, lienw,                                 &
     &                    si(:,:,:,isuevp),                             &
     &                    si(:,:,:,isvevp),                             &
     &                    si(:,:,:,isisyy),                             &
     &                    lbc(:,ibice(isisyy),ng))
!
     IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
        CALL exchange_r2d_tile (ng, tile,                               &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          si(:,:,lienw,isisxx))
 
        CALL exchange_r2d_tile (ng, tile,                               &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          si(:,:,lienw,isisxy))
 
        CALL exchange_r2d_tile (ng, tile,                               &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          si(:,:,lienw,isisyy))
      END IF
 
# ifdef DISTRIBUTE
!
      CALL mp_exchange2d (ng, tile, model, 3,                           &
     &                    lbi, ubi, lbj, ubj,                           &
     &                    nghostpoints, ewperiodic(ng), nsperiodic(ng), &
     &                    si(:,:,lienw,isisxx),                         &
     &                    si(:,:,lienw,isisxy),                         &
     &                    si(:,:,lienw,isisyy))
# endif
!
      RETURN
      END SUBROUTINE ice_evp_sig_tile
#endif
      END MODULE ice_evp_sig_mod