LCOV - code coverage report
Current view: top level - cicecore/cicedyn/dynamics - ice_dyn_evp_1d.F90 (source / functions) Hit Total Coverage
Test: 231018-211459:8916b9ff2c:1:quick Lines: 0 833 0.00 %
Date: 2023-10-18 15:30:36 Functions: 0 23 0.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : !
       3             : ! Elastic-viscous-plastic sea ice dynamics model (1D implementations)
       4             : ! Computes ice velocity and deformation
       5             : !
       6             : ! authors: Jacob Weismann Poulsen, DMI
       7             : !          Mads Hvid Ribergaard, DMI
       8             : 
       9             : module ice_dyn_evp_1d
      10             : 
      11             :    use ice_kinds_mod
      12             :    use ice_fileunits, only : nu_diag
      13             :    use ice_exit, only : abort_ice
      14             :    use icepack_intfc, only : icepack_query_parameters
      15             :    use icepack_intfc, only : icepack_warnings_flush, &
      16             :        icepack_warnings_aborted
      17             : 
      18             :    implicit none
      19             :    private
      20             :    public :: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_copyout, &
      21             :       ice_dyn_evp_1d_kernel
      22             : 
      23             :    integer(kind=int_kind) :: NA_len, NAVEL_len, domp_iam, domp_nt
      24             : #if defined (_OPENMP)
      25             :    real(kind=dbl_kind) :: rdomp_iam, rdomp_nt
      26             :    !$OMP THREADPRIVATE(domp_iam, domp_nt, rdomp_iam, rdomp_nt)
      27             : #endif
      28             :    logical(kind=log_kind), dimension(:), allocatable :: skiptcell, skipucell
      29             :    integer(kind=int_kind), dimension(:), allocatable :: ee, ne, se, &
      30             :       nw, sw, sse, indi, indj, indij, halo_parent
      31             :    real(kind=dbl_kind), dimension(:), allocatable :: cdn_ocn, aiu, &
      32             :       uocn, vocn, forcex, forcey, Tbu, tarear, umassdti, fm, uarear, &   ! LCOV_EXCL_LINE
      33             :       strintx, strinty, uvel_init, vvel_init, strength, uvel, vvel, &   ! LCOV_EXCL_LINE
      34             :       dxT, dyT, stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, &   ! LCOV_EXCL_LINE
      35             :       stressm_2, stressm_3, stressm_4, stress12_1, stress12_2, &   ! LCOV_EXCL_LINE
      36             :       stress12_3, stress12_4, divu, rdg_conv, rdg_shear, shear, taubx, &   ! LCOV_EXCL_LINE
      37             :       tauby, str1, str2, str3, str4, str5, str6, str7, str8, HTE, HTN, &   ! LCOV_EXCL_LINE
      38             :       HTEm1, HTNm1
      39             :    integer, parameter :: JPIM = selected_int_kind(9)
      40             : 
      41             :    interface evp1d_stress
      42             :       module procedure stress_iter
      43             :       module procedure stress_last
      44             :    end interface
      45             : 
      46             :    interface evp1d_stepu
      47             :       module procedure stepu_iter
      48             :       module procedure stepu_last
      49             :    end interface
      50             : 
      51             : !=======================================================================
      52             : 
      53             : contains
      54             : 
      55             : !=======================================================================
      56             : 
      57           0 :    subroutine domp_init
      58             : #if defined (_OPENMP)
      59             : 
      60             :       use omp_lib, only : omp_get_thread_num, omp_get_num_threads
      61             : #endif
      62             : 
      63             :       implicit none
      64             : 
      65             :       character(len=*), parameter :: subname = '(domp_init)'
      66             : 
      67             :       !$OMP PARALLEL DEFAULT(none)
      68             : #if defined (_OPENMP)
      69             :       domp_iam = omp_get_thread_num()
      70             :       rdomp_iam = real(domp_iam, dbl_kind)
      71             :       domp_nt = omp_get_num_threads()
      72             :       rdomp_nt = real(domp_nt, dbl_kind)
      73             : #else
      74           0 :       domp_iam = 0
      75           0 :       domp_nt = 1
      76             : #endif
      77             :       !$OMP END PARALLEL
      78             : 
      79           0 :    end subroutine domp_init
      80             : 
      81             : !=======================================================================
      82             : 
      83           0 :    subroutine domp_get_domain(lower, upper, d_lower, d_upper)
      84             : #if defined (_OPENMP)
      85             : 
      86             :       use omp_lib, only : omp_in_parallel
      87             :       use ice_constants, only : p5
      88             : #endif
      89             : 
      90             :       implicit none
      91             : 
      92             :       integer(kind=JPIM), intent(in) :: lower, upper
      93             :       integer(kind=JPIM), intent(out) :: d_lower, d_upper
      94             : 
      95             :       ! local variables
      96             : #if defined (_OPENMP)
      97             : 
      98           0 :       real(kind=dbl_kind) :: dlen
      99             : #endif
     100             : 
     101             :       character(len=*), parameter :: subname = '(domp_get_domain)'
     102             : 
     103             :       ! proper action in "null" case
     104           0 :       if (upper <= 0 .or. upper < lower) then
     105           0 :          d_lower = 0
     106           0 :          d_upper = -1
     107           0 :          return
     108             :       end if
     109             : 
     110             :       ! proper action in serial case
     111           0 :       d_lower = lower
     112           0 :       d_upper = upper
     113             : #if defined (_OPENMP)
     114             : 
     115           0 :       if (omp_in_parallel()) then
     116           0 :          dlen = real((upper - lower + 1), dbl_kind)
     117           0 :          d_lower = lower + floor(((rdomp_iam * dlen + p5) / rdomp_nt), JPIM)
     118           0 :          d_upper = lower - 1 + floor(((rdomp_iam * dlen + dlen + p5) / rdomp_nt), JPIM)
     119             :       end if
     120             : #endif
     121             : 
     122             :    end subroutine domp_get_domain
     123             : 
     124             : !=======================================================================
     125             : 
     126           0 :    subroutine stress_iter(NA_len, ee, ne, se, lb, ub, uvel, vvel, dxT, &
     127             :       dyT, hte, htn, htem1, htnm1, strength, stressp_1, stressp_2, &   ! LCOV_EXCL_LINE
     128             :       stressp_3, stressp_4, stressm_1, stressm_2, stressm_3, &   ! LCOV_EXCL_LINE
     129             :       stressm_4, stress12_1, stress12_2, stress12_3, stress12_4, str1, &   ! LCOV_EXCL_LINE
     130           0 :       str2, str3, str4, str5, str6, str7, str8, skiptcell)
     131             : 
     132             :       use ice_kinds_mod
     133             :       use ice_constants, only : p027, p055, p111, p166, p222, p25, &
     134             :           p333, p5, c1p5, c1
     135             :       use ice_dyn_shared, only : ecci, denom1, arlx1i, Ktens, revp, &
     136             :           deltaminEVP
     137             : 
     138             :       implicit none
     139             : 
     140             :       integer(kind=int_kind), intent(in) :: NA_len, lb, ub
     141             :       integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
     142             :          ee, ne, se
     143             :       real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
     144             :          strength, uvel, vvel, dxT, dyT, hte, htn, htem1, htnm1
     145             :       logical(kind=log_kind), intent(in), dimension(:) :: skiptcell
     146             :       real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
     147             :          stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, &   ! LCOV_EXCL_LINE
     148             :          stressm_2, stressm_3, stressm_4, stress12_1, stress12_2, &   ! LCOV_EXCL_LINE
     149             :          stress12_3, stress12_4
     150             :       real(kind=dbl_kind), dimension(:), intent(out), contiguous :: &
     151             :          str1, str2, str3, str4, str5, str6, str7, str8
     152             : 
     153             :       ! local variables
     154             : 
     155             :       integer(kind=int_kind) :: iw, il, iu
     156           0 :       real(kind=dbl_kind) :: divune, divunw, divuse, divusw, &
     157             :          tensionne, tensionnw, tensionse, tensionsw, shearne, shearnw, &   ! LCOV_EXCL_LINE
     158             :          shearse, shearsw, Deltane, Deltanw, Deltase, Deltasw, c0ne, &   ! LCOV_EXCL_LINE
     159             :          c0nw, c0se, c0sw, c1ne, c1nw, c1se, c1sw, ssigpn, ssigps, &   ! LCOV_EXCL_LINE
     160             :          ssigpe, ssigpw, ssigmn, ssigms, ssigme, ssigmw, ssig12n, &   ! LCOV_EXCL_LINE
     161             :          ssig12s, ssig12e, ssig12w, ssigp1, ssigp2, ssigm1, ssigm2, &   ! LCOV_EXCL_LINE
     162             :          ssig121, ssig122, csigpne, csigpnw, csigpse, csigpsw, &   ! LCOV_EXCL_LINE
     163             :          csigmne, csigmnw, csigmse, csigmsw, csig12ne, csig12nw, &   ! LCOV_EXCL_LINE
     164             :          csig12se, csig12sw, str12ew, str12we, str12ns, str12sn, &   ! LCOV_EXCL_LINE
     165             :          strp_tmp, strm_tmp, tmp_uvel_ee, tmp_vvel_se, tmp_vvel_ee, &   ! LCOV_EXCL_LINE
     166             :          tmp_vvel_ne, tmp_uvel_ne, tmp_uvel_se, dxhy, dyhx, cxp, cyp, &   ! LCOV_EXCL_LINE
     167           0 :          cxm, cym, tmparea, DminTarea
     168             : 
     169             :       character(len=*), parameter :: subname = '(stress_iter)'
     170             : 
     171             : #ifdef _OPENACC
     172             :       !$acc parallel &
     173             :       !$acc present(ee, ne, se, strength, uvel, vvel, dxT, dyT, hte, &   ! LCOV_EXCL_LINE
     174             :       !$acc    htn, htem1, htnm1, str1, str2, str3, str4, str5, str6, &   ! LCOV_EXCL_LINE
     175             :       !$acc    str7, str8, stressp_1, stressp_2, stressp_3, stressp_4, &   ! LCOV_EXCL_LINE
     176             :       !$acc    stressm_1, stressm_2, stressm_3, stressm_4, stress12_1, &   ! LCOV_EXCL_LINE
     177             :       !$acc    stress12_2, stress12_3, stress12_4, skiptcell)
     178             :       !$acc loop
     179             :       do iw = 1, NA_len
     180             : #else
     181           0 :       call domp_get_domain(lb, ub, il, iu)
     182           0 :       do iw = il, iu
     183             : #endif
     184             : 
     185           0 :          if (skiptcell(iw)) cycle
     186             : 
     187           0 :          tmparea = dxT(iw) * dyT(iw) ! necessary to split calc of DminTarea. Otherwize not binary identical
     188           0 :          DminTarea =  deltaminEVP * tmparea
     189           0 :          dxhy      = p5 * (hte(iw) - htem1(iw))
     190           0 :          dyhx      = p5 * (htn(iw) - htnm1(iw))
     191           0 :          cxp       = c1p5 * htn(iw) - p5 * htnm1(iw)
     192           0 :          cyp       = c1p5 * hte(iw) - p5 * htem1(iw)
     193           0 :          cxm       = -(c1p5 * htnm1(iw) - p5 * htn(iw))
     194           0 :          cym       = -(c1p5 * htem1(iw) - p5 * hte(iw))
     195             : 
     196             :          !--------------------------------------------------------------
     197             :          ! strain rates
     198             :          ! NOTE: these are actually strain rates * area (m^2/s)
     199             :          !--------------------------------------------------------------
     200             : 
     201           0 :          tmp_uvel_ne = uvel(ne(iw))
     202           0 :          tmp_uvel_se = uvel(se(iw))
     203           0 :          tmp_uvel_ee = uvel(ee(iw))
     204             : 
     205           0 :          tmp_vvel_ee = vvel(ee(iw))
     206           0 :          tmp_vvel_se = vvel(se(iw))
     207           0 :          tmp_vvel_ne = vvel(ne(iw))
     208             :          ! divergence = e_11 + e_22
     209           0 :          divune = cyp * uvel(iw)    - dyT(iw) * tmp_uvel_ee &
     210           0 :                 + cxp * vvel(iw)    - dxT(iw) * tmp_vvel_se
     211           0 :          divunw = cym * tmp_uvel_ee + dyT(iw) * uvel(iw) &
     212           0 :                 + cxp * tmp_vvel_ee - dxT(iw) * tmp_vvel_ne
     213           0 :          divusw = cym * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
     214           0 :                 + cxm * tmp_vvel_ne + dxT(iw) * tmp_vvel_ee
     215           0 :          divuse = cyp * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
     216           0 :                 + cxm * tmp_vvel_se + dxT(iw) * vvel(iw)
     217             : 
     218             :          ! tension strain rate = e_11 - e_22
     219           0 :          tensionne = -cym * uvel(iw)    - dyT(iw) * tmp_uvel_ee &
     220           0 :                    +  cxm * vvel(iw)    + dxT(iw) * tmp_vvel_se
     221           0 :          tensionnw = -cyp * tmp_uvel_ee + dyT(iw) * uvel(iw) &
     222           0 :                    +  cxm * tmp_vvel_ee + dxT(iw) * tmp_vvel_ne
     223           0 :          tensionsw = -cyp * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
     224           0 :                    +  cxp * tmp_vvel_ne - dxT(iw) * tmp_vvel_ee
     225           0 :          tensionse = -cym * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
     226           0 :                    +  cxp * tmp_vvel_se - dxT(iw) * vvel(iw)
     227             : 
     228             :          ! shearing strain rate = 2 * e_12
     229           0 :          shearne = -cym * vvel(iw)    - dyT(iw) * tmp_vvel_ee &
     230           0 :                  -  cxm * uvel(iw)    - dxT(iw) * tmp_uvel_se
     231           0 :          shearnw = -cyp * tmp_vvel_ee + dyT(iw) * vvel(iw) &
     232           0 :                  -  cxm * tmp_uvel_ee - dxT(iw) * tmp_uvel_ne
     233           0 :          shearsw = -cyp * tmp_vvel_ne + dyT(iw) * tmp_vvel_se &
     234           0 :                  -  cxp * tmp_uvel_ne + dxT(iw) * tmp_uvel_ee
     235           0 :          shearse = -cym * tmp_vvel_se - dyT(iw) * tmp_vvel_ne &
     236           0 :                  -  cxp * tmp_uvel_se + dxT(iw) * uvel(iw)
     237             : 
     238             :          ! Delta (in the denominator of zeta and eta)
     239           0 :          Deltane = sqrt(divune**2 + ecci * (tensionne**2 + shearne**2))
     240           0 :          Deltanw = sqrt(divunw**2 + ecci * (tensionnw**2 + shearnw**2))
     241           0 :          Deltasw = sqrt(divusw**2 + ecci * (tensionsw**2 + shearsw**2))
     242           0 :          Deltase = sqrt(divuse**2 + ecci * (tensionse**2 + shearse**2))
     243             : 
     244             :          !--------------------------------------------------------------
     245             :          ! replacement pressure/Delta (kg/s)
     246             :          ! save replacement pressure for principal stress calculation
     247             :          !--------------------------------------------------------------
     248             : 
     249           0 :          c0ne = strength(iw) / max(Deltane, DminTarea)
     250           0 :          c0nw = strength(iw) / max(Deltanw, DminTarea)
     251           0 :          c0sw = strength(iw) / max(Deltasw, DminTarea)
     252           0 :          c0se = strength(iw) / max(Deltase, DminTarea)
     253             : 
     254           0 :          c1ne = c0ne * arlx1i
     255           0 :          c1nw = c0nw * arlx1i
     256           0 :          c1sw = c0sw * arlx1i
     257           0 :          c1se = c0se * arlx1i
     258             : 
     259           0 :          c0ne = c1ne * ecci
     260           0 :          c0nw = c1nw * ecci
     261           0 :          c0sw = c1sw * ecci
     262           0 :          c0se = c1se * ecci
     263             : 
     264             :          !--------------------------------------------------------------
     265             :          ! the stresses (kg/s^2)
     266             :          ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
     267             :          !--------------------------------------------------------------
     268             : 
     269           0 :          stressp_1(iw) = (stressp_1(iw) * (c1 - arlx1i * revp) &
     270           0 :                        + c1ne * (divune * (c1 + Ktens) - Deltane * (c1 - Ktens))) * denom1
     271           0 :          stressp_2(iw) = (stressp_2(iw) * (c1 - arlx1i * revp) &
     272           0 :                        + c1nw * (divunw * (c1 + Ktens) - Deltanw * (c1 - Ktens))) * denom1
     273           0 :          stressp_3(iw) = (stressp_3(iw) * (c1 - arlx1i * revp) &
     274           0 :                        + c1sw * (divusw * (c1 + Ktens) - Deltasw * (c1 - Ktens))) * denom1
     275           0 :          stressp_4(iw) = (stressp_4(iw) * (c1 - arlx1i * revp) &
     276           0 :                        + c1se * (divuse * (c1 + Ktens) - Deltase * (c1 - Ktens))) * denom1
     277             : 
     278           0 :          stressm_1(iw) = (stressm_1(iw) * (c1 - arlx1i * revp) + c0ne * tensionne * (c1 + Ktens)) * denom1
     279           0 :          stressm_2(iw) = (stressm_2(iw) * (c1 - arlx1i * revp) + c0nw * tensionnw * (c1 + Ktens)) * denom1
     280           0 :          stressm_3(iw) = (stressm_3(iw) * (c1 - arlx1i * revp) + c0sw * tensionsw * (c1 + Ktens)) * denom1
     281           0 :          stressm_4(iw) = (stressm_4(iw) * (c1 - arlx1i * revp) + c0se * tensionse * (c1 + Ktens)) * denom1
     282             : 
     283           0 :          stress12_1(iw) = (stress12_1(iw) * (c1 - arlx1i * revp) + c0ne * shearne * p5 * (c1 + Ktens)) * denom1
     284           0 :          stress12_2(iw) = (stress12_2(iw) * (c1 - arlx1i * revp) + c0nw * shearnw * p5 * (c1 + Ktens)) * denom1
     285           0 :          stress12_3(iw) = (stress12_3(iw) * (c1 - arlx1i * revp) + c0sw * shearsw * p5 * (c1 + Ktens)) * denom1
     286           0 :          stress12_4(iw) = (stress12_4(iw) * (c1 - arlx1i * revp) + c0se * shearse * p5 * (c1 + Ktens)) * denom1
     287             : 
     288             :          !--------------------------------------------------------------
     289             :          ! combinations of the stresses for the momentum equation
     290             :          ! (kg/s^2)
     291             :          !--------------------------------------------------------------
     292             : 
     293           0 :          ssigpn =  stressp_1(iw) + stressp_2(iw)
     294           0 :          ssigps =  stressp_3(iw) + stressp_4(iw)
     295           0 :          ssigpe =  stressp_1(iw) + stressp_4(iw)
     296           0 :          ssigpw =  stressp_2(iw) + stressp_3(iw)
     297           0 :          ssigp1 = (stressp_1(iw) + stressp_3(iw)) * p055
     298           0 :          ssigp2 = (stressp_2(iw) + stressp_4(iw)) * p055
     299             : 
     300           0 :          ssigmn =  stressm_1(iw) + stressm_2(iw)
     301           0 :          ssigms =  stressm_3(iw) + stressm_4(iw)
     302           0 :          ssigme =  stressm_1(iw) + stressm_4(iw)
     303           0 :          ssigmw =  stressm_2(iw) + stressm_3(iw)
     304           0 :          ssigm1 = (stressm_1(iw) + stressm_3(iw)) * p055
     305           0 :          ssigm2 = (stressm_2(iw) + stressm_4(iw)) * p055
     306             : 
     307           0 :          ssig12n =  stress12_1(iw) + stress12_2(iw)
     308           0 :          ssig12s =  stress12_3(iw) + stress12_4(iw)
     309           0 :          ssig12e =  stress12_1(iw) + stress12_4(iw)
     310           0 :          ssig12w =  stress12_2(iw) + stress12_3(iw)
     311           0 :          ssig121 = (stress12_1(iw) + stress12_3(iw)) * p111
     312           0 :          ssig122 = (stress12_2(iw) + stress12_4(iw)) * p111
     313             : 
     314           0 :          csigpne = p111 * stressp_1(iw) + ssigp2 + p027 * stressp_3(iw)
     315           0 :          csigpnw = p111 * stressp_2(iw) + ssigp1 + p027 * stressp_4(iw)
     316           0 :          csigpsw = p111 * stressp_3(iw) + ssigp2 + p027 * stressp_1(iw)
     317           0 :          csigpse = p111 * stressp_4(iw) + ssigp1 + p027 * stressp_2(iw)
     318             : 
     319           0 :          csigmne = p111 * stressm_1(iw) + ssigm2 + p027 * stressm_3(iw)
     320           0 :          csigmnw = p111 * stressm_2(iw) + ssigm1 + p027 * stressm_4(iw)
     321           0 :          csigmsw = p111 * stressm_3(iw) + ssigm2 + p027 * stressm_1(iw)
     322           0 :          csigmse = p111 * stressm_4(iw) + ssigm1 + p027 * stressm_2(iw)
     323             : 
     324           0 :          csig12ne = p222 * stress12_1(iw) + ssig122 + p055 * stress12_3(iw)
     325           0 :          csig12nw = p222 * stress12_2(iw) + ssig121 + p055 * stress12_4(iw)
     326           0 :          csig12sw = p222 * stress12_3(iw) + ssig122 + p055 * stress12_1(iw)
     327           0 :          csig12se = p222 * stress12_4(iw) + ssig121 + p055 * stress12_2(iw)
     328             : 
     329           0 :          str12ew = p5 * dxT(iw) * (p333 * ssig12e + p166 * ssig12w)
     330           0 :          str12we = p5 * dxT(iw) * (p333 * ssig12w + p166 * ssig12e)
     331           0 :          str12ns = p5 * dyT(iw) * (p333 * ssig12n + p166 * ssig12s)
     332           0 :          str12sn = p5 * dyT(iw) * (p333 * ssig12s + p166 * ssig12n)
     333             : 
     334             :          !--------------------------------------------------------------
     335             :          ! for dF/dx (u momentum)
     336             :          !--------------------------------------------------------------
     337             : 
     338           0 :          strp_tmp = p25 * dyT(iw) * (p333 * ssigpn + p166 * ssigps)
     339           0 :          strm_tmp = p25 * dyT(iw) * (p333 * ssigmn + p166 * ssigms)
     340             : 
     341             :          ! northeast (i,j)
     342           0 :          str1(iw) = -strp_tmp - strm_tmp - str12ew &
     343           0 :                   + dxhy * (-csigpne + csigmne) + dyhx * csig12ne
     344             : 
     345             :          ! northwest (i+1,j)
     346           0 :          str2(iw) = strp_tmp + strm_tmp - str12we &
     347           0 :                   + dxhy * (-csigpnw + csigmnw) + dyhx * csig12nw
     348             : 
     349           0 :          strp_tmp = p25 * dyT(iw) * (p333 * ssigps + p166 * ssigpn)
     350           0 :          strm_tmp = p25 * dyT(iw) * (p333 * ssigms + p166 * ssigmn)
     351             : 
     352             :          ! southeast (i,j+1)
     353           0 :          str3(iw) = -strp_tmp - strm_tmp + str12ew &
     354           0 :                   + dxhy * (-csigpse + csigmse) + dyhx * csig12se
     355             : 
     356             :          ! southwest (i+1,j+1)
     357           0 :          str4(iw) = strp_tmp + strm_tmp + str12we &
     358           0 :                   + dxhy * (-csigpsw + csigmsw) + dyhx * csig12sw
     359             : 
     360             :          !--------------------------------------------------------------
     361             :          ! for dF/dy (v momentum)
     362             :          !--------------------------------------------------------------
     363             : 
     364           0 :          strp_tmp = p25 * dxT(iw) * (p333 * ssigpe + p166 * ssigpw)
     365           0 :          strm_tmp = p25 * dxT(iw) * (p333 * ssigme + p166 * ssigmw)
     366             : 
     367             :          ! northeast (i,j)
     368           0 :          str5(iw) = -strp_tmp + strm_tmp - str12ns &
     369           0 :                   - dyhx * (csigpne + csigmne) + dxhy * csig12ne
     370             : 
     371             :          ! southeast (i,j+1)
     372           0 :          str6(iw) = strp_tmp - strm_tmp - str12sn &
     373           0 :                   - dyhx * (csigpse + csigmse) + dxhy * csig12se
     374             : 
     375           0 :          strp_tmp = p25 * dxT(iw) * (p333 * ssigpw + p166 * ssigpe)
     376           0 :          strm_tmp = p25 * dxT(iw) * (p333 * ssigmw + p166 * ssigme)
     377             : 
     378             :          ! northwest (i+1,j)
     379           0 :          str7(iw) = -strp_tmp + strm_tmp + str12ns &
     380           0 :                   - dyhx * (csigpnw + csigmnw) + dxhy * csig12nw
     381             : 
     382             :          ! southwest (i+1,j+1)
     383           0 :          str8(iw) = strp_tmp - strm_tmp + str12sn &
     384           0 :                   - dyhx * (csigpsw + csigmsw) + dxhy * csig12sw
     385             : 
     386             :       end do
     387             : #ifdef _OPENACC
     388             :       !$acc end parallel
     389             : #endif
     390             : 
     391           0 :    end subroutine stress_iter
     392             : 
     393             : !=======================================================================
     394             : 
     395           0 :    subroutine stress_last(NA_len, ee, ne, se, lb, ub, uvel, vvel, dxT, &
     396             :       dyT, hte, htn, htem1, htnm1, strength, stressp_1, stressp_2, &   ! LCOV_EXCL_LINE
     397             :       stressp_3, stressp_4, stressm_1, stressm_2, stressm_3, &   ! LCOV_EXCL_LINE
     398             :       stressm_4, stress12_1, stress12_2, stress12_3, stress12_4, str1, &   ! LCOV_EXCL_LINE
     399             :       str2, str3, str4, str5, str6, str7, str8, skiptcell, tarear, divu, &   ! LCOV_EXCL_LINE
     400           0 :       rdg_conv, rdg_shear, shear)
     401             : 
     402             :       use ice_kinds_mod
     403             :       use ice_constants, only : p027, p055, p111, p166, p222, p25, &
     404             :           p333, p5, c1p5, c1, c0
     405             :       use ice_dyn_shared, only : ecci, denom1, arlx1i, Ktens, revp,&
     406             :           deltaminEVP
     407             : 
     408             :       implicit none
     409             : 
     410             :       integer(kind=int_kind), intent(in) :: NA_len, lb, ub
     411             :       integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
     412             :          ee, ne, se
     413             :       real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
     414             :          strength, uvel, vvel, dxT, dyT, hte, htn, htem1, htnm1, tarear
     415             :       logical(kind=log_kind), intent(in), dimension(:) :: skiptcell
     416             :       real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
     417             :          stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, &   ! LCOV_EXCL_LINE
     418             :          stressm_2, stressm_3, stressm_4, stress12_1, stress12_2, &   ! LCOV_EXCL_LINE
     419             :          stress12_3, stress12_4
     420             :       real(kind=dbl_kind), dimension(:), intent(out), contiguous :: &
     421             :          str1, str2, str3, str4, str5, str6, str7, str8, divu, &   ! LCOV_EXCL_LINE
     422             :          rdg_conv, rdg_shear, shear
     423             : 
     424             :       ! local variables
     425             : 
     426             :       integer(kind=int_kind) :: iw, il, iu
     427           0 :       real(kind=dbl_kind) :: divune, divunw, divuse, divusw, &
     428             :          tensionne, tensionnw, tensionse, tensionsw, shearne, shearnw, &   ! LCOV_EXCL_LINE
     429             :          shearse, shearsw, Deltane, Deltanw, Deltase, Deltasw, c0ne, &   ! LCOV_EXCL_LINE
     430             :          c0nw, c0se, c0sw, c1ne, c1nw, c1se, c1sw, ssigpn, ssigps, &   ! LCOV_EXCL_LINE
     431             :          ssigpe, ssigpw, ssigmn, ssigms, ssigme, ssigmw, ssig12n, &   ! LCOV_EXCL_LINE
     432             :          ssig12s, ssig12e, ssig12w, ssigp1, ssigp2, ssigm1, ssigm2, &   ! LCOV_EXCL_LINE
     433             :          ssig121, ssig122, csigpne, csigpnw, csigpse, csigpsw, &   ! LCOV_EXCL_LINE
     434             :          csigmne, csigmnw, csigmse, csigmsw, csig12ne, csig12nw, &   ! LCOV_EXCL_LINE
     435             :          csig12se, csig12sw, str12ew, str12we, str12ns, str12sn, &   ! LCOV_EXCL_LINE
     436             :          strp_tmp, strm_tmp, tmp_uvel_ee, tmp_vvel_se, tmp_vvel_ee, &   ! LCOV_EXCL_LINE
     437             :          tmp_vvel_ne, tmp_uvel_ne, tmp_uvel_se, dxhy, dyhx, cxp, cyp, &   ! LCOV_EXCL_LINE
     438           0 :          cxm, cym, tmparea, DminTarea
     439             : 
     440             :       character(len=*), parameter :: subname = '(stress_last)'
     441             : 
     442             : #ifdef _OPENACC
     443             :       !$acc parallel &
     444             :       !$acc present(ee, ne, se, strength, uvel, vvel, dxT, dyT, hte, &   ! LCOV_EXCL_LINE
     445             :       !$acc    htn, htem1, htnm1, str1, str2, str3, str4, str5, str6, &   ! LCOV_EXCL_LINE
     446             :       !$acc    str7, str8, stressp_1, stressp_2, stressp_3, stressp_4, &   ! LCOV_EXCL_LINE
     447             :       !$acc    stressm_1, stressm_2, stressm_3, stressm_4, stress12_1, &   ! LCOV_EXCL_LINE
     448             :       !$acc    stress12_2, stress12_3, stress12_4, tarear, divu, &   ! LCOV_EXCL_LINE
     449             :       !$acc    rdg_conv, rdg_shear, shear, skiptcell)
     450             :       !$acc loop
     451             :       do iw = 1, NA_len
     452             : #else
     453           0 :       call domp_get_domain(lb, ub, il, iu)
     454           0 :       do iw = il, iu
     455             : #endif
     456             : 
     457           0 :          if (skiptcell(iw)) cycle
     458             : 
     459           0 :          tmparea = dxT(iw) * dyT(iw) ! necessary to split calc of DminTarea. Otherwize not binary identical
     460           0 :          DminTarea = deltaminEVP * tmparea
     461           0 :          dxhy      = p5 * (hte(iw) - htem1(iw))
     462           0 :          dyhx      = p5 * (htn(iw) - htnm1(iw))
     463           0 :          cxp       = c1p5 * htn(iw) - p5 * htnm1(iw)
     464           0 :          cyp       = c1p5 * hte(iw) - p5 * htem1(iw)
     465           0 :          cxm       = -(c1p5 * htnm1(iw) - p5 * htn(iw))
     466           0 :          cym       = -(c1p5 * htem1(iw) - p5 * hte(iw))
     467             : 
     468             :          !--------------------------------------------------------------
     469             :          ! strain rates
     470             :          ! NOTE: these are actually strain rates * area (m^2/s)
     471             :          !--------------------------------------------------------------
     472             : 
     473           0 :          tmp_uvel_ne = uvel(ne(iw))
     474           0 :          tmp_uvel_se = uvel(se(iw))
     475           0 :          tmp_uvel_ee = uvel(ee(iw))
     476             : 
     477           0 :          tmp_vvel_ee = vvel(ee(iw))
     478           0 :          tmp_vvel_se = vvel(se(iw))
     479           0 :          tmp_vvel_ne = vvel(ne(iw))
     480             : 
     481             :          ! divergence = e_11 + e_22
     482           0 :          divune = cyp * uvel(iw)    - dyT(iw) * tmp_uvel_ee &
     483           0 :                 + cxp * vvel(iw)    - dxT(iw) * tmp_vvel_se
     484           0 :          divunw = cym * tmp_uvel_ee + dyT(iw) * uvel(iw) &
     485           0 :                 + cxp * tmp_vvel_ee - dxT(iw) * tmp_vvel_ne
     486           0 :          divusw = cym * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
     487           0 :                 + cxm * tmp_vvel_ne + dxT(iw) * tmp_vvel_ee
     488           0 :          divuse = cyp * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
     489           0 :                 + cxm * tmp_vvel_se + dxT(iw) * vvel(iw)
     490             : 
     491             :          ! tension strain rate = e_11 - e_22
     492           0 :          tensionne = -cym * uvel(iw)    - dyT(iw) * tmp_uvel_ee &
     493           0 :                    +  cxm * vvel(iw)    + dxT(iw) * tmp_vvel_se
     494           0 :          tensionnw = -cyp * tmp_uvel_ee + dyT(iw) * uvel(iw) &
     495           0 :                    +  cxm * tmp_vvel_ee + dxT(iw) * tmp_vvel_ne
     496           0 :          tensionsw = -cyp * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
     497           0 :                    +  cxp * tmp_vvel_ne - dxT(iw) * tmp_vvel_ee
     498           0 :          tensionse = -cym * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
     499           0 :                    +  cxp * tmp_vvel_se - dxT(iw) * vvel(iw)
     500             : 
     501             :          ! shearing strain rate = 2 * e_12
     502           0 :          shearne = -cym * vvel(iw)    - dyT(iw) * tmp_vvel_ee &
     503           0 :                  -  cxm * uvel(iw)    - dxT(iw) * tmp_uvel_se
     504           0 :          shearnw = -cyp * tmp_vvel_ee + dyT(iw) * vvel(iw) &
     505           0 :                  -  cxm * tmp_uvel_ee - dxT(iw) * tmp_uvel_ne
     506           0 :          shearsw = -cyp * tmp_vvel_ne + dyT(iw) * tmp_vvel_se &
     507           0 :                  -  cxp * tmp_uvel_ne + dxT(iw) * tmp_uvel_ee
     508           0 :          shearse = -cym * tmp_vvel_se - dyT(iw) * tmp_vvel_ne &
     509           0 :                  -  cxp * tmp_uvel_se + dxT(iw) * uvel(iw)
     510             : 
     511             :          ! Delta (in the denominator of zeta and eta)
     512           0 :          Deltane = sqrt(divune**2 + ecci * (tensionne**2 + shearne**2))
     513           0 :          Deltanw = sqrt(divunw**2 + ecci * (tensionnw**2 + shearnw**2))
     514           0 :          Deltasw = sqrt(divusw**2 + ecci * (tensionsw**2 + shearsw**2))
     515           0 :          Deltase = sqrt(divuse**2 + ecci * (tensionse**2 + shearse**2))
     516             : 
     517             :          !--------------------------------------------------------------
     518             :          ! on last subcycle, save quantities for mechanical
     519             :          ! redistribution
     520             :          !--------------------------------------------------------------
     521             : 
     522           0 :          divu(iw) = p25 * (divune + divunw + divuse + divusw) * tarear(iw)
     523           0 :          rdg_conv(iw) = -min(divu(iw), c0)  ! TODO: Could move outside the entire kernel
     524           0 :          rdg_shear(iw) = p5 * (p25 * (Deltane + Deltanw + Deltase + Deltasw) * tarear(iw) - abs(divu(iw)))
     525             : 
     526             :          ! diagnostic only
     527             :          ! shear = sqrt(tension**2 + shearing**2)
     528           0 :          shear(iw) = p25 * tarear(iw) * sqrt((tensionne + tensionnw + tensionse + tensionsw)**2 &
     529           0 :                                              + (shearne + shearnw + shearse + shearsw)**2)
     530             : 
     531             :          !--------------------------------------------------------------
     532             :          ! replacement pressure/Delta (kg/s)
     533             :          ! save replacement pressure for principal stress calculation
     534             :          !--------------------------------------------------------------
     535             : 
     536           0 :          c0ne = strength(iw) / max(Deltane, DminTarea)
     537           0 :          c0nw = strength(iw) / max(Deltanw, DminTarea)
     538           0 :          c0sw = strength(iw) / max(Deltasw, DminTarea)
     539           0 :          c0se = strength(iw) / max(Deltase, DminTarea)
     540             : 
     541           0 :          c1ne = c0ne * arlx1i
     542           0 :          c1nw = c0nw * arlx1i
     543           0 :          c1sw = c0sw * arlx1i
     544           0 :          c1se = c0se * arlx1i
     545             : 
     546           0 :          c0ne = c1ne * ecci
     547           0 :          c0nw = c1nw * ecci
     548           0 :          c0sw = c1sw * ecci
     549           0 :          c0se = c1se * ecci
     550             : 
     551             :          !--------------------------------------------------------------
     552             :          ! the stresses (kg/s^2)
     553             :          ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
     554             :          !--------------------------------------------------------------
     555             : 
     556           0 :          stressp_1(iw) = (stressp_1(iw) * (c1 - arlx1i * revp) &
     557           0 :                        + c1ne * (divune * (c1 + Ktens) - Deltane * (c1 - Ktens))) * denom1
     558           0 :          stressp_2(iw) = (stressp_2(iw) * (c1 - arlx1i * revp) &
     559           0 :                        + c1nw * (divunw * (c1 + Ktens) - Deltanw * (c1 - Ktens))) * denom1
     560           0 :          stressp_3(iw) = (stressp_3(iw) * (c1 - arlx1i * revp) &
     561           0 :                        + c1sw * (divusw * (c1 + Ktens) - Deltasw * (c1 - Ktens))) * denom1
     562           0 :          stressp_4(iw) = (stressp_4(iw) * (c1 - arlx1i * revp) &
     563           0 :                        + c1se * (divuse * (c1 + Ktens) - Deltase * (c1 - Ktens))) * denom1
     564             : 
     565           0 :          stressm_1(iw) = (stressm_1(iw) * (c1 - arlx1i * revp) + c0ne * tensionne * (c1 + Ktens)) * denom1
     566           0 :          stressm_2(iw) = (stressm_2(iw) * (c1 - arlx1i * revp) + c0nw * tensionnw * (c1 + Ktens)) * denom1
     567           0 :          stressm_3(iw) = (stressm_3(iw) * (c1 - arlx1i * revp) + c0sw * tensionsw * (c1 + Ktens)) * denom1
     568           0 :          stressm_4(iw) = (stressm_4(iw) * (c1 - arlx1i * revp) + c0se * tensionse * (c1 + Ktens)) * denom1
     569             : 
     570           0 :          stress12_1(iw) = (stress12_1(iw) * (c1 - arlx1i * revp) + c0ne * shearne * p5 * (c1 + Ktens)) * denom1
     571           0 :          stress12_2(iw) = (stress12_2(iw) * (c1 - arlx1i * revp) + c0nw * shearnw * p5 * (c1 + Ktens)) * denom1
     572           0 :          stress12_3(iw) = (stress12_3(iw) * (c1 - arlx1i * revp) + c0sw * shearsw * p5 * (c1 + Ktens)) * denom1
     573           0 :          stress12_4(iw) = (stress12_4(iw) * (c1 - arlx1i * revp) + c0se * shearse * p5 * (c1 + Ktens)) * denom1
     574             : 
     575             :          !--------------------------------------------------------------
     576             :          ! combinations of the stresses for the momentum equation
     577             :          ! (kg/s^2)
     578             :          !--------------------------------------------------------------
     579             : 
     580           0 :          ssigpn =  stressp_1(iw) + stressp_2(iw)
     581           0 :          ssigps =  stressp_3(iw) + stressp_4(iw)
     582           0 :          ssigpe =  stressp_1(iw) + stressp_4(iw)
     583           0 :          ssigpw =  stressp_2(iw) + stressp_3(iw)
     584           0 :          ssigp1 = (stressp_1(iw) + stressp_3(iw)) * p055
     585           0 :          ssigp2 = (stressp_2(iw) + stressp_4(iw)) * p055
     586             : 
     587           0 :          ssigmn =  stressm_1(iw) + stressm_2(iw)
     588           0 :          ssigms =  stressm_3(iw) + stressm_4(iw)
     589           0 :          ssigme =  stressm_1(iw) + stressm_4(iw)
     590           0 :          ssigmw =  stressm_2(iw) + stressm_3(iw)
     591           0 :          ssigm1 = (stressm_1(iw) + stressm_3(iw)) * p055
     592           0 :          ssigm2 = (stressm_2(iw) + stressm_4(iw)) * p055
     593             : 
     594           0 :          ssig12n =  stress12_1(iw) + stress12_2(iw)
     595           0 :          ssig12s =  stress12_3(iw) + stress12_4(iw)
     596           0 :          ssig12e =  stress12_1(iw) + stress12_4(iw)
     597           0 :          ssig12w =  stress12_2(iw) + stress12_3(iw)
     598           0 :          ssig121 = (stress12_1(iw) + stress12_3(iw)) * p111
     599           0 :          ssig122 = (stress12_2(iw) + stress12_4(iw)) * p111
     600             : 
     601           0 :          csigpne = p111 * stressp_1(iw) + ssigp2 + p027 * stressp_3(iw)
     602           0 :          csigpnw = p111 * stressp_2(iw) + ssigp1 + p027 * stressp_4(iw)
     603           0 :          csigpsw = p111 * stressp_3(iw) + ssigp2 + p027 * stressp_1(iw)
     604           0 :          csigpse = p111 * stressp_4(iw) + ssigp1 + p027 * stressp_2(iw)
     605             : 
     606           0 :          csigmne = p111 * stressm_1(iw) + ssigm2 + p027 * stressm_3(iw)
     607           0 :          csigmnw = p111 * stressm_2(iw) + ssigm1 + p027 * stressm_4(iw)
     608           0 :          csigmsw = p111 * stressm_3(iw) + ssigm2 + p027 * stressm_1(iw)
     609           0 :          csigmse = p111 * stressm_4(iw) + ssigm1 + p027 * stressm_2(iw)
     610             : 
     611           0 :          csig12ne = p222 * stress12_1(iw) + ssig122 + p055 * stress12_3(iw)
     612           0 :          csig12nw = p222 * stress12_2(iw) + ssig121 + p055 * stress12_4(iw)
     613           0 :          csig12sw = p222 * stress12_3(iw) + ssig122 + p055 * stress12_1(iw)
     614           0 :          csig12se = p222 * stress12_4(iw) + ssig121 + p055 * stress12_2(iw)
     615             : 
     616           0 :          str12ew = p5 * dxT(iw) * (p333 * ssig12e + p166 * ssig12w)
     617           0 :          str12we = p5 * dxT(iw) * (p333 * ssig12w + p166 * ssig12e)
     618           0 :          str12ns = p5 * dyT(iw) * (p333 * ssig12n + p166 * ssig12s)
     619           0 :          str12sn = p5 * dyT(iw) * (p333 * ssig12s + p166 * ssig12n)
     620             : 
     621             :          !--------------------------------------------------------------
     622             :          ! for dF/dx (u momentum)
     623             :          !--------------------------------------------------------------
     624             : 
     625           0 :          strp_tmp = p25 * dyT(iw) * (p333 * ssigpn + p166 * ssigps)
     626           0 :          strm_tmp = p25 * dyT(iw) * (p333 * ssigmn + p166 * ssigms)
     627             : 
     628             :          ! northeast (i,j)
     629           0 :          str1(iw) = -strp_tmp - strm_tmp - str12ew &
     630           0 :                   + dxhy * (-csigpne + csigmne) + dyhx * csig12ne
     631             : 
     632             :          ! northwest (i+1,j)
     633           0 :          str2(iw) = strp_tmp + strm_tmp - str12we &
     634           0 :                   + dxhy * (-csigpnw + csigmnw) + dyhx * csig12nw
     635             : 
     636           0 :          strp_tmp = p25 * dyT(iw) * (p333 * ssigps + p166 * ssigpn)
     637           0 :          strm_tmp = p25 * dyT(iw) * (p333 * ssigms + p166 * ssigmn)
     638             : 
     639             :          ! southeast (i,j+1)
     640           0 :          str3(iw) = -strp_tmp - strm_tmp + str12ew &
     641           0 :                   + dxhy * (-csigpse + csigmse) + dyhx * csig12se
     642             : 
     643             :          ! southwest (i+1,j+1)
     644           0 :          str4(iw) = strp_tmp + strm_tmp + str12we &
     645           0 :                   + dxhy * (-csigpsw + csigmsw) + dyhx * csig12sw
     646             : 
     647             :          !--------------------------------------------------------------
     648             :          ! for dF/dy (v momentum)
     649             :          !--------------------------------------------------------------
     650             : 
     651           0 :          strp_tmp = p25 * dxT(iw) * (p333 * ssigpe + p166 * ssigpw)
     652           0 :          strm_tmp = p25 * dxT(iw) * (p333 * ssigme + p166 * ssigmw)
     653             : 
     654             :          ! northeast (i,j)
     655           0 :          str5(iw) = -strp_tmp + strm_tmp - str12ns &
     656           0 :                   - dyhx * (csigpne + csigmne) + dxhy * csig12ne
     657             : 
     658             :          ! southeast (i,j+1)
     659           0 :          str6(iw) = strp_tmp - strm_tmp - str12sn &
     660           0 :                   - dyhx * (csigpse + csigmse) + dxhy * csig12se
     661             : 
     662           0 :          strp_tmp = p25 * dxT(iw) * (p333 * ssigpw + p166 * ssigpe)
     663           0 :          strm_tmp = p25 * dxT(iw) * (p333 * ssigmw + p166 * ssigme)
     664             : 
     665             :          ! northwest (i+1,j)
     666           0 :          str7(iw) = -strp_tmp + strm_tmp + str12ns &
     667           0 :                   - dyhx * (csigpnw + csigmnw) + dxhy * csig12nw
     668             : 
     669             :          ! southwest (i+1,j+1)
     670           0 :          str8(iw) = strp_tmp - strm_tmp + str12sn &
     671           0 :                   - dyhx * (csigpsw + csigmsw) + dxhy * csig12sw
     672             : 
     673             :       end do
     674             : #ifdef _OPENACC
     675             :       !$acc end parallel
     676             : #endif
     677             : 
     678           0 :    end subroutine stress_last
     679             : 
     680             : !=======================================================================
     681             : 
     682           0 :    subroutine stepu_iter(NA_len, rhow, lb, ub, Cw, aiu, uocn, vocn, &
     683             :       forcex, forcey, umassdti, fm, uarear, Tbu, uvel_init, vvel_init, &   ! LCOV_EXCL_LINE
     684             :       uvel, vvel, str1, str2, str3, str4, str5, str6, str7, str8, nw, &   ! LCOV_EXCL_LINE
     685           0 :       sw, sse, skipucell)
     686             : 
     687             :       use ice_kinds_mod
     688             :       use ice_constants, only : c0, c1
     689             :       use ice_dyn_shared, only : brlx, revp, u0, cosw, sinw
     690             : 
     691             :       implicit none
     692             : 
     693             :       integer(kind=int_kind), intent(in) :: NA_len, lb, ub
     694             :       real(kind=dbl_kind), intent(in) :: rhow
     695             :       logical(kind=log_kind), intent(in), dimension(:) :: skipucell
     696             :       integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
     697             :          nw, sw, sse
     698             :       real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
     699             :          uvel_init, vvel_init, aiu, forcex, forcey, umassdti, Tbu, &   ! LCOV_EXCL_LINE
     700             :          uocn, vocn, fm, uarear, Cw, str1, str2, str3, str4, str5, &   ! LCOV_EXCL_LINE
     701             :          str6, str7, str8
     702             :       real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
     703             :          uvel, vvel
     704             : 
     705             :       ! local variables
     706             : 
     707             :       integer(kind=int_kind) :: iw, il, iu
     708           0 :       real(kind=dbl_kind) :: uold, vold, vrel, cca, ccb, ab2, cc1, &
     709             :          cc2, taux, tauy, Cb, tmp_str2_nw, tmp_str3_sse, tmp_str4_sw, &   ! LCOV_EXCL_LINE
     710             :          tmp_str6_sse, tmp_str7_nw, tmp_str8_sw, waterx, watery, &   ! LCOV_EXCL_LINE
     711           0 :          tmp_strintx, tmp_strinty
     712             : 
     713             :       character(len=*), parameter :: subname = '(stepu_iter)'
     714             : 
     715             : #ifdef _OPENACC
     716             :       !$acc parallel &
     717             :       !$acc present(Cw, aiu, uocn, vocn, forcex, forcey, umassdti, fm, &   ! LCOV_EXCL_LINE
     718             :       !$acc    uarear, Tbu, uvel_init, vvel_init, nw, sw, sse, skipucell, &   ! LCOV_EXCL_LINE
     719             :       !$acc    str1, str2, str3, str4, str5, str6, str7, str8, uvel, &   ! LCOV_EXCL_LINE
     720             :       !$acc    vvel)
     721             :       !$acc loop
     722             :       do iw = 1, NA_len
     723             : #else
     724           0 :       call domp_get_domain(lb, ub, il, iu)
     725           0 :       do iw = il, iu
     726             : #endif
     727             : 
     728           0 :          if (skipucell(iw)) cycle
     729             : 
     730           0 :          uold = uvel(iw)
     731           0 :          vold = vvel(iw)
     732             : 
     733           0 :          vrel = aiu(iw) * rhow * Cw(iw) * sqrt((uocn(iw) - uold)**2 + (vocn(iw) - vold)**2)
     734             : 
     735           0 :          waterx = uocn(iw) * cosw - vocn(iw) * sinw * sign(c1, fm(iw))
     736           0 :          watery = vocn(iw) * cosw + uocn(iw) * sinw * sign(c1, fm(iw))
     737             : 
     738           0 :          taux = vrel * waterx
     739           0 :          tauy = vrel * watery
     740             : 
     741           0 :          Cb = Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
     742             : 
     743           0 :          cca = (brlx + revp) * umassdti(iw) + vrel * cosw + Cb
     744           0 :          ccb = fm(iw) + sign(c1, fm(iw)) * vrel * sinw
     745             : 
     746           0 :          ab2 = cca**2 + ccb**2
     747             : 
     748           0 :          tmp_str2_nw = str2(nw(iw))
     749           0 :          tmp_str3_sse = str3(sse(iw))
     750           0 :          tmp_str4_sw = str4(sw(iw))
     751           0 :          tmp_str6_sse = str6(sse(iw))
     752           0 :          tmp_str7_nw = str7(nw(iw))
     753           0 :          tmp_str8_sw = str8(sw(iw))
     754             : 
     755           0 :          tmp_strintx = uarear(iw) * (str1(iw) + tmp_str2_nw + tmp_str3_sse + tmp_str4_sw)
     756           0 :          tmp_strinty = uarear(iw) * (str5(iw) + tmp_str6_sse + tmp_str7_nw + tmp_str8_sw)
     757             : 
     758           0 :          cc1 = tmp_strintx + forcex(iw) + taux &
     759           0 :              + umassdti(iw) * (brlx * uold + revp * uvel_init(iw))
     760           0 :          cc2 = tmp_strinty + forcey(iw) + tauy &
     761           0 :              + umassdti(iw) * (brlx * vold + revp * vvel_init(iw))
     762             : 
     763           0 :          uvel(iw) = (cca * cc1 + ccb * cc2) / ab2
     764           0 :          vvel(iw) = (cca * cc2 - ccb * cc1) / ab2
     765             : 
     766             :       end do
     767             : #ifdef _OPENACC
     768             :       !$acc end parallel
     769             : #endif
     770             : 
     771           0 :    end subroutine stepu_iter
     772             : 
     773             : !=======================================================================
     774             : 
     775           0 :    subroutine stepu_last(NA_len, rhow, lb, ub, Cw, aiu, uocn, vocn, &
     776             :       forcex, forcey, umassdti, fm, uarear, Tbu, uvel_init, vvel_init, &   ! LCOV_EXCL_LINE
     777             :       uvel, vvel, str1, str2, str3, str4, str5, str6, str7, str8, nw, &   ! LCOV_EXCL_LINE
     778           0 :       sw, sse, skipucell, strintx, strinty, taubx, tauby)
     779             : 
     780             :       use ice_kinds_mod
     781             :       use ice_constants, only : c0, c1
     782             :       use ice_dyn_shared, only : brlx, revp, u0, cosw, sinw
     783             : 
     784             :       implicit none
     785             : 
     786             :       integer(kind=int_kind), intent(in) :: NA_len, lb, ub
     787             :       real(kind=dbl_kind), intent(in) :: rhow
     788             :       logical(kind=log_kind), intent(in), dimension(:) :: skipucell
     789             :       integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
     790             :          nw, sw, sse
     791             :       real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
     792             :          uvel_init, vvel_init, aiu, forcex, forcey, umassdti, Tbu, &   ! LCOV_EXCL_LINE
     793             :          uocn, vocn, fm, uarear, Cw, str1, str2, str3, str4, str5, &   ! LCOV_EXCL_LINE
     794             :          str6, str7, str8
     795             :       real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
     796             :          uvel, vvel, strintx, strinty, taubx, tauby
     797             : 
     798             :       ! local variables
     799             : 
     800             :       integer(kind=int_kind) :: iw, il, iu
     801           0 :       real(kind=dbl_kind) :: uold, vold, vrel, cca, ccb, ab2, cc1, &
     802             :          cc2, taux, tauy, Cb, tmp_str2_nw, tmp_str3_sse, tmp_str4_sw, &   ! LCOV_EXCL_LINE
     803           0 :          tmp_str6_sse, tmp_str7_nw, tmp_str8_sw, waterx, watery
     804             : 
     805             :       character(len=*), parameter :: subname = '(stepu_last)'
     806             : 
     807             : #ifdef _OPENACC
     808             :       !$acc parallel &
     809             :       !$acc present(Cw, aiu, uocn, vocn, forcex, forcey, umassdti, fm, &   ! LCOV_EXCL_LINE
     810             :       !$acc    uarear, Tbu, uvel_init, vvel_init, nw, sw, sse, skipucell, &   ! LCOV_EXCL_LINE
     811             :       !$acc    str1, str2, str3, str4, str5, str6, str7, str8, uvel, &   ! LCOV_EXCL_LINE
     812             :       !$acc    vvel, strintx, strinty, taubx, tauby)
     813             :       !$acc loop
     814             :       do iw = 1, NA_len
     815             : #else
     816           0 :       call domp_get_domain(lb, ub, il, iu)
     817           0 :       do iw = il, iu
     818             : #endif
     819             : 
     820           0 :          if (skipucell(iw)) cycle
     821             : 
     822           0 :          uold = uvel(iw)
     823           0 :          vold = vvel(iw)
     824             : 
     825           0 :          vrel = aiu(iw) * rhow * Cw(iw) * sqrt((uocn(iw) - uold)**2 + (vocn(iw) - vold)**2)
     826             : 
     827           0 :          waterx = uocn(iw) * cosw - vocn(iw) * sinw * sign(c1, fm(iw))
     828           0 :          watery = vocn(iw) * cosw + uocn(iw) * sinw * sign(c1, fm(iw))
     829             : 
     830           0 :          taux = vrel * waterx
     831           0 :          tauy = vrel * watery
     832             : 
     833           0 :          Cb = Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
     834             : 
     835           0 :          cca = (brlx + revp) * umassdti(iw) + vrel * cosw + Cb
     836           0 :          ccb = fm(iw) + sign(c1, fm(iw)) * vrel * sinw
     837             : 
     838           0 :          ab2 = cca**2 + ccb**2
     839             : 
     840           0 :          tmp_str2_nw = str2(nw(iw))
     841           0 :          tmp_str3_sse = str3(sse(iw))
     842           0 :          tmp_str4_sw = str4(sw(iw))
     843           0 :          tmp_str6_sse = str6(sse(iw))
     844           0 :          tmp_str7_nw = str7(nw(iw))
     845           0 :          tmp_str8_sw = str8(sw(iw))
     846             : 
     847           0 :          strintx(iw) = uarear(iw) * (str1(iw) + tmp_str2_nw + tmp_str3_sse + tmp_str4_sw)
     848           0 :          strinty(iw) = uarear(iw) * (str5(iw) + tmp_str6_sse + tmp_str7_nw + tmp_str8_sw)
     849             : 
     850           0 :          cc1 = strintx(iw) + forcex(iw) + taux &
     851           0 :              + umassdti(iw) * (brlx * uold + revp * uvel_init(iw))
     852           0 :          cc2 = strinty(iw) + forcey(iw) + tauy &
     853           0 :              + umassdti(iw) * (brlx * vold + revp * vvel_init(iw))
     854             : 
     855           0 :          uvel(iw) = (cca * cc1 + ccb * cc2) / ab2
     856           0 :          vvel(iw) = (cca * cc2 - ccb * cc1) / ab2
     857             : 
     858             :          ! calculate seabed stress component for outputs
     859           0 :          taubx(iw) = -uvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
     860           0 :          tauby(iw) = -vvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
     861             : 
     862             :       end do
     863             : #ifdef _OPENACC
     864             :       !$acc end parallel
     865             : #endif
     866             : 
     867           0 :    end subroutine stepu_last
     868             : 
     869             : !=======================================================================
     870             : 
     871           0 :    subroutine evp1d_halo_update(NAVEL_len, lb, ub, uvel, vvel, &
     872           0 :       halo_parent)
     873             : 
     874             :       use ice_kinds_mod
     875             : 
     876             :       implicit none
     877             : 
     878             :       integer(kind=int_kind), intent(in) :: NAVEL_len, lb, ub
     879             :       integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
     880             :          halo_parent
     881             :       real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
     882             :          uvel, vvel
     883             : 
     884             :       ! local variables
     885             : 
     886             :       integer (kind=int_kind) :: iw, il, iu
     887             : 
     888             :       character(len=*), parameter :: subname = '(evp1d_halo_update)'
     889             : 
     890             : #ifdef _OPENACC
     891             :       !$acc parallel &
     892             :       !$acc present(uvel, vvel)
     893             :       !$acc loop
     894             :       do iw = 1, NAVEL_len
     895             :          if (halo_parent(iw) == 0) cycle
     896             :          uvel(iw) = uvel(halo_parent(iw))
     897             :          vvel(iw) = vvel(halo_parent(iw))
     898             :       end do
     899             :       !$acc end parallel
     900             : #else
     901           0 :       call domp_get_domain(lb, ub, il, iu)
     902           0 :       do iw = il, iu
     903           0 :          if (halo_parent(iw) == 0) cycle
     904           0 :          uvel(iw) = uvel(halo_parent(iw))
     905           0 :          vvel(iw) = vvel(halo_parent(iw))
     906             :       end do
     907           0 :       call domp_get_domain(ub + 1, NAVEL_len, il, iu)
     908           0 :       do iw = il, iu
     909           0 :          if (halo_parent(iw) == 0) cycle
     910           0 :          uvel(iw) = uvel(halo_parent(iw))
     911           0 :          vvel(iw) = vvel(halo_parent(iw))
     912             :       end do
     913             : #endif
     914             : 
     915           0 :    end subroutine evp1d_halo_update
     916             : 
     917             : !=======================================================================
     918             : 
     919           0 :    subroutine alloc1d(na)
     920             : 
     921             :       implicit none
     922             : 
     923             :       integer(kind=int_kind), intent(in) :: na
     924             : 
     925             :       ! local variables
     926             : 
     927             :       integer(kind=int_kind) :: ierr
     928             : 
     929             :       character(len=*), parameter :: subname = '(alloc1d)'
     930             : 
     931             :       allocate( &
     932             :          ! helper indices for neighbours
     933             :          indj(1:na), indi(1:na), ee(1:na), ne(1:na), se(1:na), &
     934             :          nw(1:na), sw(1:na), sse(1:na), skipucell(1:na), &   ! LCOV_EXCL_LINE
     935             :          skiptcell(1:na), &   ! LCOV_EXCL_LINE
     936             :          ! grid distances and their "-1 neighbours"
     937             :          HTE(1:na), HTN(1:na), HTEm1(1:na), HTNm1(1:na), &
     938             :          ! T cells
     939             :          strength(1:na), dxT(1:na), dyT(1:na), tarear(1:na), &
     940             :          stressp_1(1:na), stressp_2(1:na), stressp_3(1:na), &   ! LCOV_EXCL_LINE
     941             :          stressp_4(1:na), stressm_1(1:na), stressm_2(1:na), &   ! LCOV_EXCL_LINE
     942             :          stressm_3(1:na), stressm_4(1:na), stress12_1(1:na), &   ! LCOV_EXCL_LINE
     943             :          stress12_2(1:na), stress12_3(1:na), stress12_4(1:na), &   ! LCOV_EXCL_LINE
     944             :          divu(1:na), rdg_conv(1:na), rdg_shear(1:na), shear(1:na), &   ! LCOV_EXCL_LINE
     945             :          ! U cells
     946             :          cdn_ocn(1:na), aiu(1:na), uocn(1:na), vocn(1:na), &
     947             :          forcex(1:na), forcey(1:na), Tbu(1:na), umassdti(1:na), &   ! LCOV_EXCL_LINE
     948             :          fm(1:na), uarear(1:na), strintx(1:na), strinty(1:na), &   ! LCOV_EXCL_LINE
     949             :          uvel_init(1:na), vvel_init(1:na), taubx(1:na), tauby(1:na), &   ! LCOV_EXCL_LINE
     950             :          ! error handling
     951           0 :          stat=ierr)
     952             : 
     953           0 :       if (ierr /= 0) call abort_ice(subname &
     954           0 :          // ' ERROR: could not allocate 1D arrays')
     955             : 
     956           0 :    end subroutine alloc1d
     957             : 
     958             : !=======================================================================
     959             : 
     960           0 :    subroutine alloc1d_navel(navel)
     961             : 
     962             :       implicit none
     963             : 
     964             :       integer(kind=int_kind), intent(in) :: navel
     965             : 
     966             :       ! local variables
     967             : 
     968             :       integer(kind=int_kind) :: ierr
     969             : 
     970             :       character(len=*), parameter :: subname = '(alloc1d_navel)'
     971             : 
     972             :       allocate(uvel(1:navel), vvel(1:navel), indij(1:navel), &
     973             :          halo_parent(1:navel), str1(1:navel), str2(1:navel), &   ! LCOV_EXCL_LINE
     974             :          str3(1:navel), str4(1:navel), str5(1:navel), str6(1:navel), &   ! LCOV_EXCL_LINE
     975           0 :          str7(1:navel), str8(1:navel), stat=ierr)
     976             : 
     977           0 :       if (ierr /= 0) call abort_ice(subname &
     978           0 :          // ' ERROR: could not allocate 1D arrays')
     979             : 
     980           0 :    end subroutine alloc1d_navel
     981             : 
     982             : !=======================================================================
     983             : 
     984           0 :    subroutine dealloc1d
     985             : 
     986             :       implicit none
     987             : 
     988             :       ! local variables
     989             : 
     990             :       integer(kind=int_kind) :: ierr
     991             : 
     992             :       character(len=*), parameter :: subname = '(dealloc1d)'
     993             : 
     994             :       deallocate( &
     995             :          ! helper indices for neighbours
     996             :          indj, indi, ee, ne, se, nw, sw, sse, skipucell, skiptcell, &
     997             :          ! grid distances and their "-1 neighbours"
     998             :          HTE, HTN, HTEm1, HTNm1, &
     999             :          ! T cells
    1000             :          strength, dxT, dyT, tarear, stressp_1, stressp_2, stressp_3, &
    1001             :          stressp_4, stressm_1, stressm_2, stressm_3, stressm_4, &   ! LCOV_EXCL_LINE
    1002             :          stress12_1, stress12_2, stress12_3, stress12_4, str1, str2, &   ! LCOV_EXCL_LINE
    1003             :          str3, str4, str5, str6, str7, str8, divu, rdg_conv, &   ! LCOV_EXCL_LINE
    1004             :          rdg_shear, shear, &   ! LCOV_EXCL_LINE
    1005             :          ! U cells
    1006             :          cdn_ocn, aiu, uocn, vocn, forcex, forcey, Tbu, umassdti, fm, &
    1007             :          uarear, strintx, strinty, uvel_init, vvel_init, taubx, tauby, &   ! LCOV_EXCL_LINE
    1008             :          uvel, vvel, indij, halo_parent, &   ! LCOV_EXCL_LINE
    1009             :          ! error handling
    1010           0 :          stat=ierr)
    1011             : 
    1012           0 :       if (ierr /= 0) call abort_ice(subname &
    1013           0 :          // ' ERROR: could not deallocate 1D arrays')
    1014             : 
    1015           0 :    end subroutine dealloc1d
    1016             : 
    1017             : !=======================================================================
    1018             : 
    1019           0 :    subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, &
    1020             :       I_iceTmask, I_iceUmask, I_cdn_ocn, I_aiu, I_uocn, I_vocn, &   ! LCOV_EXCL_LINE
    1021             :       I_forcex, I_forcey, I_Tbu, I_umassdti, I_fm, I_uarear, I_tarear, &   ! LCOV_EXCL_LINE
    1022             :       I_strintx, I_strinty, I_uvel_init, I_vvel_init, I_strength, &   ! LCOV_EXCL_LINE
    1023             :       I_uvel, I_vvel, I_dxT, I_dyT, I_stressp_1, I_stressp_2, &   ! LCOV_EXCL_LINE
    1024             :       I_stressp_3, I_stressp_4, I_stressm_1, I_stressm_2, I_stressm_3, &   ! LCOV_EXCL_LINE
    1025             :       I_stressm_4, I_stress12_1, I_stress12_2, I_stress12_3, &   ! LCOV_EXCL_LINE
    1026           0 :       I_stress12_4)
    1027             : 
    1028             :       use ice_gather_scatter, only : gather_global_ext
    1029             :       use ice_domain, only : distrb_info
    1030             :       use ice_communicate, only : my_task, master_task
    1031             :       use ice_grid, only : G_HTE, G_HTN
    1032             :       use ice_constants, only : c0
    1033             : 
    1034             :       implicit none
    1035             : 
    1036             :       integer(int_kind), intent(in) :: nx, ny, nblk, nx_glob, ny_glob
    1037             :       logical(kind=log_kind), dimension(nx, ny, nblk), intent(in) :: &
    1038             :          I_iceTmask, I_iceUmask
    1039             :       real(kind=dbl_kind), dimension(nx, ny, nblk), intent(in) :: &
    1040             :          I_cdn_ocn, I_aiu, I_uocn, I_vocn, I_forcex, I_forcey, I_Tbu, &   ! LCOV_EXCL_LINE
    1041             :          I_umassdti, I_fm, I_uarear, I_tarear, I_strintx, I_strinty, &   ! LCOV_EXCL_LINE
    1042             :          I_uvel_init, I_vvel_init, I_strength, I_uvel, I_vvel, I_dxT, &   ! LCOV_EXCL_LINE
    1043             :          I_dyT, I_stressp_1, I_stressp_2, I_stressp_3, I_stressp_4, &   ! LCOV_EXCL_LINE
    1044             :          I_stressm_1, I_stressm_2, I_stressm_3, I_stressm_4, &   ! LCOV_EXCL_LINE
    1045             :          I_stress12_1, I_stress12_2, I_stress12_3, I_stress12_4
    1046             : 
    1047             :       ! local variables
    1048             : 
    1049             :       logical(kind=log_kind), dimension(nx_glob, ny_glob) :: &
    1050           0 :          G_iceTmask, G_iceUmask
    1051             :       real(kind=dbl_kind), dimension(nx_glob, ny_glob) :: &
    1052             :          G_cdn_ocn, G_aiu, G_uocn, G_vocn, G_forcex, G_forcey, G_Tbu, &   ! LCOV_EXCL_LINE
    1053             :          G_umassdti, G_fm, G_uarear, G_tarear, G_strintx, G_strinty, &   ! LCOV_EXCL_LINE
    1054             :          G_uvel_init, G_vvel_init, G_strength, G_uvel, G_vvel, G_dxT, &   ! LCOV_EXCL_LINE
    1055             :          G_dyT, G_stressp_1, G_stressp_2, G_stressp_3, G_stressp_4, &   ! LCOV_EXCL_LINE
    1056             :          G_stressm_1, G_stressm_2, G_stressm_3, G_stressm_4, &   ! LCOV_EXCL_LINE
    1057           0 :          G_stress12_1, G_stress12_2, G_stress12_3, G_stress12_4
    1058             : 
    1059             :       character(len=*), parameter :: &
    1060             :          subname = '(ice_dyn_evp_1d_copyin)'
    1061             : 
    1062           0 :       call gather_global_ext(G_iceTmask,   I_iceTmask,   master_task, distrb_info    )
    1063           0 :       call gather_global_ext(G_iceUmask,   I_iceUmask,   master_task, distrb_info    )
    1064           0 :       call gather_global_ext(G_cdn_ocn,    I_cdn_ocn,    master_task, distrb_info    )
    1065           0 :       call gather_global_ext(G_aiu,        I_aiu,        master_task, distrb_info    )
    1066           0 :       call gather_global_ext(G_uocn,       I_uocn,       master_task, distrb_info    )
    1067           0 :       call gather_global_ext(G_vocn,       I_vocn,       master_task, distrb_info    )
    1068           0 :       call gather_global_ext(G_forcex,     I_forcex,     master_task, distrb_info    )
    1069           0 :       call gather_global_ext(G_forcey,     I_forcey,     master_task, distrb_info    )
    1070           0 :       call gather_global_ext(G_Tbu,        I_Tbu,        master_task, distrb_info    )
    1071           0 :       call gather_global_ext(G_umassdti,   I_umassdti,   master_task, distrb_info    )
    1072           0 :       call gather_global_ext(G_fm,         I_fm,         master_task, distrb_info    )
    1073           0 :       call gather_global_ext(G_uarear,     I_uarear,     master_task, distrb_info    )
    1074           0 :       call gather_global_ext(G_tarear,     I_tarear,     master_task, distrb_info    )
    1075           0 :       call gather_global_ext(G_strintx,    I_strintx,    master_task, distrb_info    )
    1076           0 :       call gather_global_ext(G_strinty,    I_strinty,    master_task, distrb_info    )
    1077           0 :       call gather_global_ext(G_uvel_init,  I_uvel_init,  master_task, distrb_info    )
    1078           0 :       call gather_global_ext(G_vvel_init,  I_vvel_init,  master_task, distrb_info    )
    1079           0 :       call gather_global_ext(G_strength,   I_strength,   master_task, distrb_info    )
    1080           0 :       call gather_global_ext(G_uvel,       I_uvel,       master_task, distrb_info, c0)
    1081           0 :       call gather_global_ext(G_vvel,       I_vvel,       master_task, distrb_info, c0)
    1082           0 :       call gather_global_ext(G_dxT,        I_dxT,        master_task, distrb_info    )
    1083           0 :       call gather_global_ext(G_dyT,        I_dyT,        master_task, distrb_info    )
    1084           0 :       call gather_global_ext(G_stressp_1,  I_stressp_1,  master_task, distrb_info    )
    1085           0 :       call gather_global_ext(G_stressp_2,  I_stressp_2,  master_task, distrb_info    )
    1086           0 :       call gather_global_ext(G_stressp_3,  I_stressp_3,  master_task, distrb_info    )
    1087           0 :       call gather_global_ext(G_stressp_4,  I_stressp_4,  master_task, distrb_info    )
    1088           0 :       call gather_global_ext(G_stressm_1,  I_stressm_1,  master_task, distrb_info    )
    1089           0 :       call gather_global_ext(G_stressm_2,  I_stressm_2,  master_task, distrb_info    )
    1090           0 :       call gather_global_ext(G_stressm_3,  I_stressm_3,  master_task, distrb_info    )
    1091           0 :       call gather_global_ext(G_stressm_4,  I_stressm_4,  master_task, distrb_info    )
    1092           0 :       call gather_global_ext(G_stress12_1, I_stress12_1, master_task, distrb_info    )
    1093           0 :       call gather_global_ext(G_stress12_2, I_stress12_2, master_task, distrb_info    )
    1094           0 :       call gather_global_ext(G_stress12_3, I_stress12_3, master_task, distrb_info    )
    1095           0 :       call gather_global_ext(G_stress12_4, I_stress12_4, master_task, distrb_info    )
    1096             : 
    1097             :       ! all calculations id done on master task
    1098           0 :       if (my_task == master_task) then
    1099             :          ! find number of active points and allocate 1D vectors
    1100           0 :          call calc_na(nx_glob, ny_glob, NA_len, G_iceTmask, G_iceUmask)
    1101           0 :          call alloc1d(NA_len)
    1102           0 :          call calc_2d_indices(nx_glob, ny_glob, NA_len, G_iceTmask, G_iceUmask)
    1103           0 :          call calc_navel(nx_glob, ny_glob, NA_len, NAVEL_len)
    1104           0 :          call alloc1d_navel(NAVEL_len)
    1105             :          ! initialize OpenMP. FIXME: ought to be called from main
    1106           0 :          call domp_init()
    1107             :          !$OMP PARALLEL DEFAULT(shared)
    1108           0 :          call numainit(1, NA_len, NAVEL_len)
    1109             :          !$OMP END PARALLEL
    1110             :          ! map 2D arrays to 1D arrays
    1111             :          call convert_2d_1d(nx_glob, ny_glob, NA_len, NAVEL_len, &
    1112             :             G_HTE, G_HTN, G_cdn_ocn, G_aiu, G_uocn, G_vocn, G_forcex, &   ! LCOV_EXCL_LINE
    1113             :             G_forcey, G_Tbu, G_umassdti, G_fm, G_uarear, G_tarear, &   ! LCOV_EXCL_LINE
    1114             :             G_strintx, G_strinty, G_uvel_init, G_vvel_init, &   ! LCOV_EXCL_LINE
    1115             :             G_strength, G_uvel, G_vvel, G_dxT, G_dyT, G_stressp_1, &   ! LCOV_EXCL_LINE
    1116             :             G_stressp_2, G_stressp_3, G_stressp_4, G_stressm_1, &   ! LCOV_EXCL_LINE
    1117             :             G_stressm_2, G_stressm_3, G_stressm_4, G_stress12_1, &   ! LCOV_EXCL_LINE
    1118           0 :             G_stress12_2, G_stress12_3, G_stress12_4)
    1119           0 :          call calc_halo_parent(nx_glob, ny_glob, NA_len, NAVEL_len, G_iceTmask)
    1120             :       end if
    1121             : 
    1122           0 :    end subroutine ice_dyn_evp_1d_copyin
    1123             : 
    1124             : !=======================================================================
    1125             : 
    1126           0 :    subroutine ice_dyn_evp_1d_copyout(nx, ny, nblk, nx_glob, ny_glob, &
    1127             :       I_uvel, I_vvel, I_strintx, I_strinty, I_stressp_1, I_stressp_2, &   ! LCOV_EXCL_LINE
    1128             :       I_stressp_3, I_stressp_4, I_stressm_1, I_stressm_2, I_stressm_3, &   ! LCOV_EXCL_LINE
    1129             :       I_stressm_4, I_stress12_1, I_stress12_2, I_stress12_3, &   ! LCOV_EXCL_LINE
    1130             :       I_stress12_4, I_divu, I_rdg_conv, I_rdg_shear, I_shear, I_taubx, &   ! LCOV_EXCL_LINE
    1131           0 :       I_tauby)
    1132             : 
    1133             :       use ice_constants, only : c0
    1134             :       use ice_gather_scatter, only : scatter_global_ext
    1135             :       use ice_domain, only : distrb_info
    1136             :       use ice_communicate, only : my_task, master_task
    1137             : 
    1138             :       implicit none
    1139             : 
    1140             :       integer(int_kind), intent(in) :: nx, ny, nblk, nx_glob, ny_glob
    1141             :       real(dbl_kind), dimension(nx, ny, nblk), intent(out) :: I_uvel, &
    1142             :          I_vvel, I_strintx, I_strinty, I_stressp_1, I_stressp_2, &   ! LCOV_EXCL_LINE
    1143             :          I_stressp_3, I_stressp_4, I_stressm_1, I_stressm_2, &   ! LCOV_EXCL_LINE
    1144             :          I_stressm_3, I_stressm_4, I_stress12_1, I_stress12_2, &   ! LCOV_EXCL_LINE
    1145             :          I_stress12_3, I_stress12_4, I_divu, I_rdg_conv, I_rdg_shear, &   ! LCOV_EXCL_LINE
    1146             :          I_shear, I_taubx, I_tauby
    1147             : 
    1148             :       ! local variables
    1149             : 
    1150             :       integer(int_kind) :: iw, lo, up, j, i
    1151           0 :       real(dbl_kind), dimension(nx_glob, ny_glob) :: G_uvel, G_vvel, &
    1152             :          G_strintx, G_strinty, G_stressp_1, G_stressp_2, G_stressp_3, &   ! LCOV_EXCL_LINE
    1153             :          G_stressp_4, G_stressm_1, G_stressm_2, G_stressm_3, &   ! LCOV_EXCL_LINE
    1154             :          G_stressm_4, G_stress12_1, G_stress12_2, G_stress12_3, &   ! LCOV_EXCL_LINE
    1155             :          G_stress12_4, G_divu, G_rdg_conv, G_rdg_shear, G_shear, &   ! LCOV_EXCL_LINE
    1156           0 :          G_taubx, G_tauby
    1157             : 
    1158             :       character(len=*), parameter :: &
    1159             :          subname = '(ice_dyn_evp_1d_copyout)'
    1160             : 
    1161             :       ! remap 1D arrays into 2D arrays
    1162           0 :       if (my_task == master_task) then
    1163             : 
    1164           0 :          G_uvel       = c0
    1165           0 :          G_vvel       = c0
    1166           0 :          G_strintx    = c0
    1167           0 :          G_strinty    = c0
    1168           0 :          G_stressp_1  = c0
    1169           0 :          G_stressp_2  = c0
    1170           0 :          G_stressp_3  = c0
    1171           0 :          G_stressp_4  = c0
    1172           0 :          G_stressm_1  = c0
    1173           0 :          G_stressm_2  = c0
    1174           0 :          G_stressm_3  = c0
    1175           0 :          G_stressm_4  = c0
    1176           0 :          G_stress12_1 = c0
    1177           0 :          G_stress12_2 = c0
    1178           0 :          G_stress12_3 = c0
    1179           0 :          G_stress12_4 = c0
    1180           0 :          G_divu       = c0
    1181           0 :          G_rdg_conv   = c0
    1182           0 :          G_rdg_shear  = c0
    1183           0 :          G_shear      = c0
    1184           0 :          G_taubx      = c0
    1185           0 :          G_tauby      = c0
    1186             : 
    1187           0 :          !$OMP PARALLEL PRIVATE(iw, lo, up, j, i)
    1188           0 :          call domp_get_domain(1, NA_len, lo, up)
    1189           0 :          do iw = lo, up
    1190             :             ! get 2D indices
    1191           0 :             i = indi(iw)
    1192           0 :             j = indj(iw)
    1193             :             ! remap
    1194           0 :             G_strintx(i, j)    = strintx(iw)
    1195           0 :             G_strinty(i, j)    = strinty(iw)
    1196           0 :             G_stressp_1(i, j)  = stressp_1(iw)
    1197           0 :             G_stressp_2(i, j)  = stressp_2(iw)
    1198           0 :             G_stressp_3(i, j)  = stressp_3(iw)
    1199           0 :             G_stressp_4(i, j)  = stressp_4(iw)
    1200           0 :             G_stressm_1(i, j)  = stressm_1(iw)
    1201           0 :             G_stressm_2(i, j)  = stressm_2(iw)
    1202           0 :             G_stressm_3(i, j)  = stressm_3(iw)
    1203           0 :             G_stressm_4(i, j)  = stressm_4(iw)
    1204           0 :             G_stress12_1(i, j) = stress12_1(iw)
    1205           0 :             G_stress12_2(i, j) = stress12_2(iw)
    1206           0 :             G_stress12_3(i, j) = stress12_3(iw)
    1207           0 :             G_stress12_4(i, j) = stress12_4(iw)
    1208           0 :             G_divu(i, j)       = divu(iw)
    1209           0 :             G_rdg_conv(i, j)   = rdg_conv(iw)
    1210           0 :             G_rdg_shear(i, j)  = rdg_shear(iw)
    1211           0 :             G_shear(i, j)      = shear(iw)
    1212           0 :             G_taubx(i, j)      = taubx(iw)
    1213           0 :             G_tauby(i, j)      = tauby(iw)
    1214           0 :             G_uvel(i, j) = uvel(iw)
    1215           0 :             G_vvel(i, j) = vvel(iw)
    1216             :          end do
    1217           0 :          call domp_get_domain(NA_len + 1, NAVEL_len, lo, up)
    1218           0 :          do iw = lo, up
    1219             :             ! get 2D indices
    1220           0 :             j = int((indij(iw) - 1) / (nx_glob)) + 1
    1221           0 :             i = indij(iw) - (j - 1) * nx_glob
    1222             :             ! remap
    1223           0 :             G_uvel(i, j) = uvel(iw)
    1224           0 :             G_vvel(i, j) = vvel(iw)
    1225             :          end do
    1226             :          !$OMP END PARALLEL
    1227             : 
    1228           0 :          call dealloc1d()
    1229             : 
    1230             :       end if
    1231             : 
    1232             :       ! scatter data on all tasks
    1233           0 :       call scatter_global_ext(I_uvel,       G_uvel,       master_task, distrb_info)
    1234           0 :       call scatter_global_ext(I_vvel,       G_vvel,       master_task, distrb_info)
    1235           0 :       call scatter_global_ext(I_strintx,    G_strintx,    master_task, distrb_info)
    1236           0 :       call scatter_global_ext(I_strinty,    G_strinty,    master_task, distrb_info)
    1237           0 :       call scatter_global_ext(I_stressp_1,  G_stressp_1,  master_task, distrb_info)
    1238           0 :       call scatter_global_ext(I_stressp_2,  G_stressp_2,  master_task, distrb_info)
    1239           0 :       call scatter_global_ext(I_stressp_3,  G_stressp_3,  master_task, distrb_info)
    1240           0 :       call scatter_global_ext(I_stressp_4,  G_stressp_4,  master_task, distrb_info)
    1241           0 :       call scatter_global_ext(I_stressm_1,  G_stressm_1,  master_task, distrb_info)
    1242           0 :       call scatter_global_ext(I_stressm_2,  G_stressm_2,  master_task, distrb_info)
    1243           0 :       call scatter_global_ext(I_stressm_3,  G_stressm_3,  master_task, distrb_info)
    1244           0 :       call scatter_global_ext(I_stressm_4,  G_stressm_4,  master_task, distrb_info)
    1245           0 :       call scatter_global_ext(I_stress12_1, G_stress12_1, master_task, distrb_info)
    1246           0 :       call scatter_global_ext(I_stress12_2, G_stress12_2, master_task, distrb_info)
    1247           0 :       call scatter_global_ext(I_stress12_3, G_stress12_3, master_task, distrb_info)
    1248           0 :       call scatter_global_ext(I_stress12_4, G_stress12_4, master_task, distrb_info)
    1249           0 :       call scatter_global_ext(I_divu,       G_divu,       master_task, distrb_info)
    1250           0 :       call scatter_global_ext(I_rdg_conv,   G_rdg_conv,   master_task, distrb_info)
    1251           0 :       call scatter_global_ext(I_rdg_shear,  G_rdg_shear,  master_task, distrb_info)
    1252           0 :       call scatter_global_ext(I_shear,      G_shear,      master_task, distrb_info)
    1253           0 :       call scatter_global_ext(I_taubx,      G_taubx,      master_task, distrb_info)
    1254           0 :       call scatter_global_ext(I_tauby,      G_tauby,      master_task, distrb_info)
    1255             : 
    1256           0 :    end subroutine ice_dyn_evp_1d_copyout
    1257             : 
    1258             : !=======================================================================
    1259             : 
    1260           0 :    subroutine ice_dyn_evp_1d_kernel
    1261             : 
    1262             :       use ice_constants, only : c0
    1263             :       use ice_dyn_shared, only : ndte
    1264             :       use ice_communicate, only : my_task, master_task
    1265             : 
    1266             :       implicit none
    1267             : 
    1268             :       ! local variables
    1269             : 
    1270           0 :       real(kind=dbl_kind) :: rhow
    1271             :       integer(kind=int_kind) :: ksub
    1272             : 
    1273             :       character(len=*), parameter :: &
    1274             :          subname = '(ice_dyn_evp_1d_kernel)'
    1275             : 
    1276             :       ! all calculations is done on master task
    1277           0 :       if (my_task == master_task) then
    1278             : 
    1279             :          ! read constants
    1280           0 :          call icepack_query_parameters(rhow_out = rhow)
    1281           0 :          call icepack_warnings_flush(nu_diag)
    1282           0 :          if (icepack_warnings_aborted()) then
    1283             :             call abort_ice(error_message=subname, file=__FILE__, &
    1284           0 :                line=__LINE__)
    1285             :          end if
    1286             : 
    1287           0 :          if (ndte < 2) call abort_ice(subname &
    1288           0 :             // ' ERROR: ndte must be 2 or higher for this kernel')
    1289             : 
    1290             :          ! tcraig, turn off the OMP directives here, Jan, 2022
    1291             :          ! This produces non bit-for-bit results with different thread counts.
    1292             :          ! Seems like there isn't an opportunity for safe threading here ???
    1293             :          !$XXXOMP PARALLEL PRIVATE(ksub)
    1294           0 :          do ksub = 1, ndte - 1
    1295             :             call evp1d_stress(NA_len, ee, ne, se, 1, NA_len, uvel, &
    1296             :                vvel, dxT, dyT, hte, htn, htem1, htnm1, strength, &   ! LCOV_EXCL_LINE
    1297             :                stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, &   ! LCOV_EXCL_LINE
    1298             :                stressm_2, stressm_3, stressm_4, stress12_1, &   ! LCOV_EXCL_LINE
    1299             :                stress12_2, stress12_3, stress12_4, str1, str2, str3, &   ! LCOV_EXCL_LINE
    1300           0 :                str4, str5, str6, str7, str8, skiptcell)
    1301             :             !$XXXOMP BARRIER
    1302             :             call evp1d_stepu(NA_len, rhow, 1, NA_len, cdn_ocn, aiu, &
    1303             :                uocn, vocn, forcex, forcey, umassdti, fm, uarear, Tbu, &   ! LCOV_EXCL_LINE
    1304             :                uvel_init, vvel_init, uvel, vvel, str1, str2, str3, &   ! LCOV_EXCL_LINE
    1305           0 :                str4, str5, str6, str7, str8, nw, sw, sse, skipucell)
    1306             :             !$XXXOMP BARRIER
    1307             :             call evp1d_halo_update(NAVEL_len, 1, NA_len, uvel, vvel, &
    1308           0 :                halo_parent)
    1309             :             !$XXXOMP BARRIER
    1310             :          end do
    1311             : 
    1312             :          call evp1d_stress(NA_len, ee, ne, se, 1, NA_len, uvel, vvel, &
    1313             :             dxT, dyT, hte, htn, htem1, htnm1, strength, stressp_1, &   ! LCOV_EXCL_LINE
    1314             :             stressp_2, stressp_3, stressp_4, stressm_1, stressm_2, &   ! LCOV_EXCL_LINE
    1315             :             stressm_3, stressm_4, stress12_1, stress12_2, stress12_3, &   ! LCOV_EXCL_LINE
    1316             :             stress12_4, str1, str2, str3, str4, str5, str6, str7, &   ! LCOV_EXCL_LINE
    1317           0 :             str8, skiptcell, tarear, divu, rdg_conv, rdg_shear, shear)
    1318             :          !$XXXOMP BARRIER
    1319             :          call evp1d_stepu(NA_len, rhow, 1, NA_len, cdn_ocn, aiu, uocn, &
    1320             :             vocn, forcex, forcey, umassdti, fm, uarear, Tbu, &   ! LCOV_EXCL_LINE
    1321             :             uvel_init, vvel_init, uvel, vvel, str1, str2, str3, str4, &   ! LCOV_EXCL_LINE
    1322             :             str5, str6, str7, str8, nw, sw, sse, skipucell, strintx, &   ! LCOV_EXCL_LINE
    1323           0 :             strinty, taubx, tauby)
    1324             :          !$XXXOMP BARRIER
    1325             :          call evp1d_halo_update(NAVEL_len, 1, NA_len, uvel, vvel, &
    1326           0 :             halo_parent)
    1327             :          !$XXXOMP END PARALLEL
    1328             : 
    1329             :       end if  ! master task
    1330             : 
    1331           0 :    end subroutine ice_dyn_evp_1d_kernel
    1332             : 
    1333             : !=======================================================================
    1334             : 
    1335           0 :    subroutine calc_na(nx, ny, na, iceTmask, iceUmask)
    1336             :       ! Calculate number of active points
    1337             : 
    1338             :       use ice_blocks, only : nghost
    1339             : 
    1340             :       implicit none
    1341             : 
    1342             :       integer(kind=int_kind), intent(in) :: nx, ny
    1343             :       logical(kind=log_kind), dimension(nx, ny), intent(in) :: &
    1344             :          iceTmask, iceUmask
    1345             :       integer(kind=int_kind), intent(out) :: na
    1346             : 
    1347             :       ! local variables
    1348             : 
    1349             :       integer(kind=int_kind) :: i, j
    1350             : 
    1351             :       character(len=*), parameter :: subname = '(calc_na)'
    1352             : 
    1353           0 :       na = 0
    1354             :       ! NOTE: T mask includes northern and eastern ghost cells
    1355           0 :       do j = 1 + nghost, ny
    1356           0 :          do i = 1 + nghost, nx
    1357           0 :             if (iceTmask(i,j) .or. iceUmask(i,j)) na = na + 1
    1358             :          end do
    1359             :       end do
    1360             : 
    1361           0 :    end subroutine calc_na
    1362             : 
    1363             : !=======================================================================
    1364             : 
    1365           0 :    subroutine calc_2d_indices(nx, ny, na, iceTmask, iceUmask)
    1366             : 
    1367             :       use ice_blocks, only : nghost
    1368             : 
    1369             :       implicit none
    1370             : 
    1371             :       integer(kind=int_kind), intent(in) :: nx, ny, na
    1372             :       logical(kind=log_kind), dimension(nx, ny), intent(in) :: &
    1373             :          iceTmask, iceUmask
    1374             : 
    1375             :       ! local variables
    1376             : 
    1377             :       integer(kind=int_kind) :: i, j, Nmaskt
    1378             : 
    1379             :       character(len=*), parameter :: subname = '(calc_2d_indices)'
    1380             : 
    1381           0 :       skipucell(:) = .false.
    1382           0 :       skiptcell(:) = .false.
    1383           0 :       indi = 0
    1384           0 :       indj = 0
    1385           0 :       Nmaskt = 0
    1386             :       ! NOTE: T mask includes northern and eastern ghost cells
    1387           0 :       do j = 1 + nghost, ny
    1388           0 :          do i = 1 + nghost, nx
    1389           0 :             if (iceTmask(i,j) .or. iceUmask(i,j)) then
    1390           0 :                Nmaskt = Nmaskt + 1
    1391           0 :                indi(Nmaskt) = i
    1392           0 :                indj(Nmaskt) = j
    1393           0 :                if (.not. iceTmask(i,j)) skiptcell(Nmaskt) = .true.
    1394           0 :                if (.not. iceUmask(i,j)) skipucell(Nmaskt) = .true.
    1395             :                ! NOTE: U mask does not include northern and eastern
    1396             :                ! ghost cells. Skip northern and eastern ghost cells
    1397           0 :                if (i == nx) skipucell(Nmaskt) = .true.
    1398           0 :                if (j == ny) skipucell(Nmaskt) = .true.
    1399             :             end if
    1400             :          end do
    1401             :       end do
    1402             : 
    1403           0 :    end subroutine calc_2d_indices
    1404             : 
    1405             : !=======================================================================
    1406             : 
    1407           0 :    subroutine calc_navel(nx_block, ny_block, na, navel)
    1408             :       ! Calculate number of active points, including halo points
    1409             : 
    1410             :       implicit none
    1411             : 
    1412             :       integer(kind=int_kind), intent(in) :: nx_block, ny_block, na
    1413             :       integer(kind=int_kind), intent(out) :: navel
    1414             : 
    1415             :       ! local variables
    1416             : 
    1417             :       integer(kind=int_kind) :: iw, i, j
    1418           0 :       integer(kind=int_kind), dimension(1:na) :: Iin, Iee, Ine, Ise, &
    1419           0 :          Inw, Isw, Isse
    1420           0 :       integer(kind=int_kind), dimension(1:7 * na) :: util1, util2
    1421             : 
    1422             :       character(len=*), parameter :: subname = '(calc_navel)'
    1423             : 
    1424             :       ! calculate additional 1D indices used for finite differences
    1425           0 :       do iw = 1, na
    1426             :          ! get 2D indices
    1427           0 :          i = indi(iw)
    1428           0 :          j = indj(iw)
    1429             :          ! calculate 1D indices
    1430           0 :          Iin(iw)  = i     + (j - 1) * nx_block  ! ( 0,  0) target point
    1431           0 :          Iee(iw)  = i - 1 + (j - 1) * nx_block  ! (-1,  0)
    1432           0 :          Ine(iw)  = i - 1 + (j - 2) * nx_block  ! (-1, -1)
    1433           0 :          Ise(iw)  = i     + (j - 2) * nx_block  ! ( 0, -1)
    1434           0 :          Inw(iw)  = i + 1 + (j - 1) * nx_block  ! (+1,  0)
    1435           0 :          Isw(iw)  = i + 1 + (j - 0) * nx_block  ! (+1, +1)
    1436           0 :          Isse(iw) = i +     (j - 0) * nx_block  ! ( 0, +1)
    1437             :       end do
    1438             : 
    1439             :       ! find number of points needed for finite difference calculations
    1440           0 :       call union(Iin,   Iee,  na, na, util1, i    )
    1441           0 :       call union(util1, Ine,  i,  na, util2, j    )
    1442           0 :       call union(util2, Ise,  j,  na, util1, i    )
    1443           0 :       call union(util1, Inw,  i,  na, util2, j    )
    1444           0 :       call union(util2, Isw,  j,  na, util1, i    )
    1445           0 :       call union(util1, Isse, i,  na, util2, navel)
    1446             : 
    1447           0 :    end subroutine calc_navel
    1448             : 
    1449             : !=======================================================================
    1450             : 
    1451           0 :    subroutine convert_2d_1d(nx, ny, na, navel, I_HTE, I_HTN, &
    1452             :       I_cdn_ocn, I_aiu, I_uocn, I_vocn, I_forcex, I_forcey, I_Tbu, &   ! LCOV_EXCL_LINE
    1453             :       I_umassdti, I_fm, I_uarear, I_tarear, I_strintx, I_strinty, &   ! LCOV_EXCL_LINE
    1454             :       I_uvel_init, I_vvel_init, I_strength, I_uvel, I_vvel, I_dxT, &   ! LCOV_EXCL_LINE
    1455             :       I_dyT, I_stressp_1, I_stressp_2, I_stressp_3, I_stressp_4, &   ! LCOV_EXCL_LINE
    1456             :       I_stressm_1, I_stressm_2, I_stressm_3, I_stressm_4, &   ! LCOV_EXCL_LINE
    1457           0 :       I_stress12_1, I_stress12_2, I_stress12_3, I_stress12_4)
    1458             : 
    1459             :       implicit none
    1460             : 
    1461             :       integer(kind=int_kind), intent(in) :: nx, ny, na, navel
    1462             :       real (kind=dbl_kind), dimension(nx, ny), intent(in) :: I_HTE, &
    1463             :          I_HTN, I_cdn_ocn, I_aiu, I_uocn, I_vocn, I_forcex, I_forcey, &   ! LCOV_EXCL_LINE
    1464             :          I_Tbu, I_umassdti, I_fm, I_uarear, I_tarear, I_strintx, &   ! LCOV_EXCL_LINE
    1465             :          I_strinty, I_uvel_init, I_vvel_init, I_strength, I_uvel, &   ! LCOV_EXCL_LINE
    1466             :          I_vvel, I_dxT, I_dyT, I_stressp_1, I_stressp_2, I_stressp_3, &   ! LCOV_EXCL_LINE
    1467             :          I_stressp_4, I_stressm_1, I_stressm_2, I_stressm_3, &   ! LCOV_EXCL_LINE
    1468             :          I_stressm_4, I_stress12_1, I_stress12_2, I_stress12_3, &   ! LCOV_EXCL_LINE
    1469             :          I_stress12_4
    1470             : 
    1471             :       ! local variables
    1472             : 
    1473             :       integer(kind=int_kind) :: iw, lo, up, j, i, nachk
    1474           0 :       integer(kind=int_kind), dimension(1:na) :: Iin, Iee, Ine, Ise, &
    1475           0 :          Inw, Isw, Isse
    1476           0 :       integer(kind=int_kind), dimension(1:7 * na) :: util1, util2
    1477             : 
    1478             :       character(len=*), parameter :: subname = '(convert_2d_1d)'
    1479             : 
    1480             :       ! calculate additional 1D indices used for finite differences
    1481           0 :       do iw = 1, na
    1482             :          ! get 2D indices
    1483           0 :          i = indi(iw)
    1484           0 :          j = indj(iw)
    1485             :          ! calculate 1D indices
    1486           0 :          Iin(iw)  = i     + (j - 1) * nx  ! ( 0, 0) target point
    1487           0 :          Iee(iw)  = i - 1 + (j - 1) * nx  ! (-1, 0)
    1488           0 :          Ine(iw)  = i - 1 + (j - 2) * nx  ! (-1,-1)
    1489           0 :          Ise(iw)  = i     + (j - 2) * nx  ! ( 0,-1)
    1490           0 :          Inw(iw)  = i + 1 + (j - 1) * nx  ! (+1, 0)
    1491           0 :          Isw(iw)  = i + 1 + (j - 0) * nx  ! (+1,+1)
    1492           0 :          Isse(iw) = i     + (j - 0) * nx  ! ( 0,+1)
    1493             :       end do
    1494             : 
    1495             :       ! find number of points needed for finite difference calculations
    1496           0 :       call union(Iin,   Iee,  na, na, util1, i    )
    1497           0 :       call union(util1, Ine,  i,  na, util2, j    )
    1498           0 :       call union(util2, Ise,  j,  na, util1, i    )
    1499           0 :       call union(util1, Inw,  i,  na, util2, j    )
    1500           0 :       call union(util2, Isw,  j,  na, util1, i    )
    1501           0 :       call union(util1, Isse, i,  na, util2, nachk)
    1502             : 
    1503             :       ! index vector with sorted target points
    1504           0 :       do iw = 1, na
    1505           0 :          indij(iw) = Iin(iw)
    1506             :       end do
    1507             : 
    1508             :       ! sorted additional points
    1509           0 :       call setdiff(util2, Iin, navel, na, util1, j)
    1510           0 :       do iw = na + 1, navel
    1511           0 :          indij(iw) = util1(iw - na)
    1512             :       end do
    1513             : 
    1514             :       ! indices for additional points needed for uvel and vvel
    1515           0 :       call findXinY(Iee,  indij, na, navel, ee)
    1516           0 :       call findXinY(Ine,  indij, na, navel, ne)
    1517           0 :       call findXinY(Ise,  indij, na, navel, se)
    1518           0 :       call findXinY(Inw,  indij, na, navel, nw)
    1519           0 :       call findXinY(Isw,  indij, na, navel, sw)
    1520           0 :       call findXinY(Isse, indij, na, navel, sse)
    1521             : 
    1522           0 :       !$OMP PARALLEL PRIVATE(iw, lo, up, j, i)
    1523             :       ! write 1D arrays from 2D arrays (target points)
    1524           0 :       call domp_get_domain(1, na, lo, up)
    1525           0 :       do iw = lo, up
    1526             :          ! get 2D indices
    1527           0 :          i = indi(iw)
    1528           0 :          j = indj(iw)
    1529             :          ! map
    1530           0 :          uvel(iw)       = I_uvel(i, j)
    1531           0 :          vvel(iw)       = I_vvel(i, j)
    1532           0 :          cdn_ocn(iw)    = I_cdn_ocn(i, j)
    1533           0 :          aiu(iw)        = I_aiu(i, j)
    1534           0 :          uocn(iw)       = I_uocn(i, j)
    1535           0 :          vocn(iw)       = I_vocn(i, j)
    1536           0 :          forcex(iw)     = I_forcex(i, j)
    1537           0 :          forcey(iw)     = I_forcey(i, j)
    1538           0 :          Tbu(iw)        = I_Tbu(i, j)
    1539           0 :          umassdti(iw)   = I_umassdti(i, j)
    1540           0 :          fm(iw)         = I_fm(i, j)
    1541           0 :          tarear(iw)     = I_tarear(i, j)
    1542           0 :          uarear(iw)     = I_uarear(i, j)
    1543           0 :          strintx(iw)    = I_strintx(i, j)
    1544           0 :          strinty(iw)    = I_strinty(i, j)
    1545           0 :          uvel_init(iw)  = I_uvel_init(i, j)
    1546           0 :          vvel_init(iw)  = I_vvel_init(i, j)
    1547           0 :          strength(iw)   = I_strength(i, j)
    1548           0 :          dxT(iw)        = I_dxT(i, j)
    1549           0 :          dyT(iw)        = I_dyT(i, j)
    1550           0 :          stressp_1(iw)  = I_stressp_1(i, j)
    1551           0 :          stressp_2(iw)  = I_stressp_2(i, j)
    1552           0 :          stressp_3(iw)  = I_stressp_3(i, j)
    1553           0 :          stressp_4(iw)  = I_stressp_4(i, j)
    1554           0 :          stressm_1(iw)  = I_stressm_1(i, j)
    1555           0 :          stressm_2(iw)  = I_stressm_2(i, j)
    1556           0 :          stressm_3(iw)  = I_stressm_3(i, j)
    1557           0 :          stressm_4(iw)  = I_stressm_4(i, j)
    1558           0 :          stress12_1(iw) = I_stress12_1(i, j)
    1559           0 :          stress12_2(iw) = I_stress12_2(i, j)
    1560           0 :          stress12_3(iw) = I_stress12_3(i, j)
    1561           0 :          stress12_4(iw) = I_stress12_4(i, j)
    1562           0 :          HTE(iw)        = I_HTE(i, j)
    1563           0 :          HTN(iw)        = I_HTN(i, j)
    1564           0 :          HTEm1(iw)      = I_HTE(i - 1, j)
    1565           0 :          HTNm1(iw)      = I_HTN(i, j - 1)
    1566             :       end do
    1567             :       ! write 1D arrays from 2D arrays (additional points)
    1568           0 :       call domp_get_domain(na + 1, navel, lo, up)
    1569           0 :       do iw = lo, up
    1570             :          ! get 2D indices
    1571           0 :          j = int((indij(iw) - 1) / (nx)) + 1
    1572           0 :          i = indij(iw) - (j - 1) * nx
    1573             :          ! map
    1574           0 :          uvel(iw) = I_uvel(i, j)
    1575           0 :          vvel(iw) = I_vvel(i, j)
    1576             :       end do
    1577             :       !$OMP END PARALLEL
    1578             : 
    1579           0 :    end subroutine convert_2d_1d
    1580             : 
    1581             : !=======================================================================
    1582             : 
    1583           0 :    subroutine calc_halo_parent(nx, ny, na, navel, I_iceTmask)
    1584             : 
    1585             :       implicit none
    1586             : 
    1587             :       integer(kind=int_kind), intent(in) :: nx, ny, na, navel
    1588             :       logical(kind=log_kind), dimension(nx, ny), intent(in) :: &
    1589             :          I_iceTmask
    1590             : 
    1591             :       ! local variables
    1592             : 
    1593             :       integer(kind=int_kind) :: iw, i, j
    1594           0 :       integer(kind=int_kind), dimension(1:navel) :: Ihalo
    1595             : 
    1596             :       character(len=*), parameter :: subname = '(calc_halo_parent)'
    1597             : 
    1598             :       !-----------------------------------------------------------------
    1599             :       ! Indices for halo update:
    1600             :       !     0: no halo point
    1601             :       !    >0: index for halo point parent, related to indij vector
    1602             :       !
    1603             :       ! TODO: Implement for nghost > 1
    1604             :       ! TODO: Implement for tripole grids
    1605             :       !-----------------------------------------------------------------
    1606             : 
    1607           0 :       Ihalo(:) = 0
    1608           0 :       halo_parent(:) = 0
    1609             : 
    1610           0 :       do iw = 1, navel
    1611           0 :          j = int((indij(iw) - 1) / (nx)) + 1
    1612           0 :          i = indij(iw) - (j - 1) * nx
    1613             :          ! if within ghost zone
    1614           0 :          if (i == nx .and. I_iceTmask(2, j)      ) Ihalo(iw) = 2 + (j - 1) * nx
    1615           0 :          if (i == 1  .and. I_iceTmask(nx - 1, j) ) Ihalo(iw) = (nx - 1) + (j - 1) * nx
    1616           0 :          if (j == ny .and. I_iceTmask(i, 2)      ) Ihalo(iw) = i + nx
    1617           0 :          if (j == 1  .and. I_iceTmask(i, ny - 1) ) Ihalo(iw) = i + (ny - 2) * nx
    1618             :       end do
    1619             : 
    1620             :       ! relate halo indices to indij vector
    1621           0 :       call findXinY_halo(Ihalo, indij, navel, navel, halo_parent)
    1622             : 
    1623           0 :    end subroutine calc_halo_parent
    1624             : 
    1625             : !=======================================================================
    1626             : 
    1627           0 :    subroutine union(x, y, nx, ny, xy, nxy)
    1628             :       ! Find union (xy) of two sorted integer vectors (x and y), i.e.
    1629             :       ! combined values of the two vectors with no repetitions
    1630             : 
    1631             :       implicit none
    1632             : 
    1633             :       integer(int_kind), intent(in) :: nx, ny
    1634             :       integer(int_kind), intent(in) :: x(1:nx), y(1:ny)
    1635             :       integer(int_kind), intent(out) :: xy(1:nx + ny)
    1636             :       integer(int_kind), intent(out) :: nxy
    1637             : 
    1638             :       ! local variables
    1639             : 
    1640             :       integer(int_kind) :: i, j, k
    1641             : 
    1642             :       character(len=*), parameter :: subname = '(union)'
    1643             : 
    1644           0 :       i = 1
    1645           0 :       j = 1
    1646           0 :       k = 1
    1647           0 :       do while (i <= nx .and. j <= ny)
    1648           0 :          if (x(i) < y(j)) then
    1649           0 :             xy(k) = x(i)
    1650           0 :             i = i + 1
    1651           0 :          else if (x(i) > y(j)) then
    1652           0 :             xy(k) = y(j)
    1653           0 :             j = j + 1
    1654             :          else
    1655           0 :             xy(k) = x(i)
    1656           0 :             i = i + 1
    1657           0 :             j = j + 1
    1658             :          end if
    1659           0 :          k = k + 1
    1660             :       end do
    1661             : 
    1662             :       ! the rest
    1663           0 :       do while (i <= nx)
    1664           0 :          xy(k) = x(i)
    1665           0 :          i = i + 1
    1666           0 :          k = k + 1
    1667             :       end do
    1668           0 :       do while (j <= ny)
    1669           0 :          xy(k) = y(j)
    1670           0 :          j = j + 1
    1671           0 :          k = k + 1
    1672             :       end do
    1673           0 :       nxy = k - 1
    1674             : 
    1675           0 :    end subroutine union
    1676             : 
    1677             : !=======================================================================
    1678             : 
    1679           0 :    subroutine setdiff(x, y, nx, ny, xy, nxy)
    1680             :       ! Find element (xy) of two sorted integer vectors (x and y) that
    1681             :       ! are in x, but not in y, or in y, but not in x
    1682             : 
    1683             :       implicit none
    1684             : 
    1685             :       integer(kind=int_kind), intent(in) :: nx, ny
    1686             :       integer(kind=int_kind), intent(in) :: x(1:nx), y(1:ny)
    1687             :       integer(kind=int_kind), intent(out) :: xy(1:nx + ny)
    1688             :       integer(kind=int_kind), intent(out) :: nxy
    1689             : 
    1690             :       ! local variables
    1691             : 
    1692             :       integer(kind=int_kind) :: i, j, k
    1693             : 
    1694             :       character(len=*), parameter :: subname = '(setdiff)'
    1695             : 
    1696           0 :       i = 1
    1697           0 :       j = 1
    1698           0 :       k = 1
    1699           0 :       do while (i <= nx .and. j <= ny)
    1700           0 :          if (x(i) < y(j)) then
    1701           0 :             xy(k) = x(i)
    1702           0 :             i = i + 1
    1703           0 :             k = k + 1
    1704           0 :          else if (x(i) > y(j)) then
    1705           0 :             xy(k) = y(j)
    1706           0 :             j = j + 1
    1707           0 :             k = k + 1
    1708             :          else
    1709           0 :             i = i + 1
    1710           0 :             j = j + 1
    1711             :          end if
    1712             :       end do
    1713             : 
    1714             :       ! the rest
    1715           0 :       do while (i <= nx)
    1716           0 :          xy(k) = x(i)
    1717           0 :          i = i + 1
    1718           0 :          k = k + 1
    1719             :       end do
    1720           0 :       do while (j <= ny)
    1721           0 :          xy(k) = y(j)
    1722           0 :          j = j + 1
    1723           0 :          k = k + 1
    1724             :       end do
    1725           0 :       nxy = k - 1
    1726             : 
    1727           0 :    end subroutine setdiff
    1728             : 
    1729             : !========================================================================
    1730             : 
    1731           0 :    subroutine findXinY(x, y, nx, ny, indx)
    1732             :       ! Find indx vector so that x(1:na) = y(indx(1:na))
    1733             :       !
    1734             :       !  Conditions:
    1735             :       !   * EVERY item in x is found in y
    1736             :       !   * x(1:nx) is a sorted integer vector
    1737             :       !   * y(1:ny) consists of two sorted integer vectors:
    1738             :       !        [y(1:nx); y(nx + 1:ny)]
    1739             :       !   * ny >= nx
    1740             :       !
    1741             :       !  Return: indx(1:na)
    1742             : 
    1743             :       implicit none
    1744             : 
    1745             :       integer (kind=int_kind), intent(in) :: nx, ny
    1746             :       integer (kind=int_kind), intent(in) :: x(1:nx), y(1:ny)
    1747             :       integer (kind=int_kind), intent(out) :: indx(1:nx)
    1748             : 
    1749             :       ! local variables
    1750             : 
    1751             :       integer (kind=int_kind) :: i, j1, j2
    1752             : 
    1753             :       character(len=*), parameter :: subname = '(findXinY)'
    1754             : 
    1755           0 :       i = 1
    1756           0 :       j1 = 1
    1757           0 :       j2 = nx + 1
    1758           0 :       do while (i <= nx)
    1759           0 :          if (x(i) == y(j1)) then
    1760           0 :             indx(i) = j1
    1761           0 :             i = i + 1
    1762           0 :             j1 = j1 + 1
    1763           0 :          else if (x(i) == y(j2)) then
    1764           0 :             indx(i) = j2
    1765           0 :             i = i + 1
    1766           0 :             j2 = j2 + 1
    1767           0 :          else if (x(i) > y(j1)) then
    1768           0 :             j1 = j1 + 1
    1769           0 :          else if (x(i) > y(j2)) then
    1770           0 :             j2 = j2 + 1
    1771             :          else
    1772             :             call abort_ice(subname &
    1773           0 :                // ': ERROR: conditions not met')
    1774             :          end if
    1775             :       end do
    1776             : 
    1777           0 :    end subroutine findXinY
    1778             : 
    1779             : !=======================================================================
    1780             : 
    1781           0 :    subroutine findXinY_halo(x, y, nx, ny, indx)
    1782             :       ! Find indx vector so that x(1:na) = y(indx(1:na))
    1783             :       !
    1784             :       !  Conditions:
    1785             :       !   * EVERY item in x is found in y,
    1786             :       !        except for x == 0, where indx = 0 is returned
    1787             :       !   * x(1:nx) is a non-sorted integer vector
    1788             :       !   * y(1:ny) is a sorted integer vector
    1789             :       !   * ny >= nx
    1790             :       !
    1791             :       !  Return: indx(1:na)
    1792             : 
    1793             :       implicit none
    1794             : 
    1795             :       integer (kind=int_kind), intent(in) :: nx, ny
    1796             :       integer (kind=int_kind), intent(in) :: x(1:nx), y(1:ny)
    1797             :       integer (kind=int_kind), intent(out) :: indx(1:nx)
    1798             : 
    1799             :       ! local variables
    1800             : 
    1801             :       integer (kind=int_kind) :: i, j1, nloop
    1802             : 
    1803             :       character(len=*), parameter :: subname = '(findXinY_halo)'
    1804             : 
    1805           0 :       nloop = 1
    1806           0 :       i = 1
    1807           0 :       j1 = int((ny + 1) / 2)  ! initial guess in the middle
    1808           0 :       do while (i <= nx)
    1809           0 :          if (x(i) == 0) then
    1810           0 :             indx(i) = 0
    1811           0 :             i = i + 1
    1812           0 :             nloop = 1
    1813           0 :          else if (x(i) == y(j1)) then
    1814           0 :             indx(i) = j1
    1815           0 :             i = i + 1
    1816           0 :             j1 = j1 + 1
    1817             :             ! initial guess in the middle
    1818           0 :             if (j1 > ny) j1 = int((ny + 1) / 2)
    1819           0 :             nloop = 1
    1820           0 :          else if (x(i) < y(j1)) then
    1821           0 :             j1 = 1
    1822           0 :          else if (x(i) > y(j1)) then
    1823           0 :             j1 = j1 + 1
    1824           0 :             if (j1 > ny) then
    1825           0 :                j1 = 1
    1826           0 :                nloop = nloop + 1
    1827           0 :                if (nloop > 2) then
    1828             :                   ! stop for infinite loop. This check should not be
    1829             :                   ! necessary for halo
    1830           0 :                   call abort_ice(subname // ' ERROR: too many loops')
    1831             :                end if
    1832             :             end if
    1833             :          end if
    1834             :       end do
    1835             : 
    1836           0 :    end subroutine findXinY_halo
    1837             : 
    1838             : !=======================================================================
    1839             : 
    1840           0 :    subroutine numainit(l, u, uu)
    1841             : 
    1842             :       use ice_constants, only : c0
    1843             : 
    1844             :       implicit none
    1845             : 
    1846             :       integer(kind=int_kind), intent(in) :: l, u, uu
    1847             : 
    1848             :       ! local variables
    1849             : 
    1850             :       integer(kind=int_kind) :: lo, up
    1851             : 
    1852             :       character(len=*), parameter :: subname = '(numainit)'
    1853             : 
    1854           0 :       call domp_get_domain(l, u, lo, up)
    1855           0 :       ee(lo:up)          = 0
    1856           0 :       ne(lo:up)          = 0
    1857           0 :       se(lo:up)          = 0
    1858           0 :       sse(lo:up)         = 0
    1859           0 :       nw(lo:up)          = 0
    1860           0 :       sw(lo:up)          = 0
    1861           0 :       halo_parent(lo:up) = 0
    1862           0 :       strength(lo:up)    = c0
    1863           0 :       uvel(lo:up)        = c0
    1864           0 :       vvel(lo:up)        = c0
    1865           0 :       uvel_init(lo:up)   = c0
    1866           0 :       vvel_init(lo:up)   = c0
    1867           0 :       uocn(lo:up)        = c0
    1868           0 :       vocn(lo:up)        = c0
    1869           0 :       dxT(lo:up)         = c0
    1870           0 :       dyT(lo:up)         = c0
    1871           0 :       HTE(lo:up)         = c0
    1872           0 :       HTN(lo:up)         = c0
    1873           0 :       HTEm1(lo:up)       = c0
    1874           0 :       HTNm1(lo:up)       = c0
    1875           0 :       stressp_1(lo:up)   = c0
    1876           0 :       stressp_2(lo:up)   = c0
    1877           0 :       stressp_3(lo:up)   = c0
    1878           0 :       stressp_4(lo:up)   = c0
    1879           0 :       stressm_1(lo:up)   = c0
    1880           0 :       stressm_2(lo:up)   = c0
    1881           0 :       stressm_3(lo:up)   = c0
    1882           0 :       stressm_4(lo:up)   = c0
    1883           0 :       stress12_1(lo:up)  = c0
    1884           0 :       stress12_2(lo:up)  = c0
    1885           0 :       stress12_3(lo:up)  = c0
    1886           0 :       stress12_4(lo:up)  = c0
    1887           0 :       tarear(lo:up)      = c0
    1888           0 :       Tbu(lo:up)         = c0
    1889           0 :       taubx(lo:up)       = c0
    1890           0 :       tauby(lo:up)       = c0
    1891           0 :       divu(lo:up)        = c0
    1892           0 :       rdg_conv(lo:up)    = c0
    1893           0 :       rdg_shear(lo:up)   = c0
    1894           0 :       shear(lo:up)       = c0
    1895           0 :       str1(lo:up)        = c0
    1896           0 :       str2(lo:up)        = c0
    1897           0 :       str3(lo:up)        = c0
    1898           0 :       str4(lo:up)        = c0
    1899           0 :       str5(lo:up)        = c0
    1900           0 :       str6(lo:up)        = c0
    1901           0 :       str7(lo:up)        = c0
    1902           0 :       str8(lo:up)        = c0
    1903             : 
    1904           0 :       call domp_get_domain(u + 1, uu, lo, up)
    1905           0 :       halo_parent(lo:up) = 0
    1906           0 :       uvel(lo:up)        = c0
    1907           0 :       vvel(lo:up)        = c0
    1908           0 :       str1(lo:up)        = c0
    1909           0 :       str2(lo:up)        = c0
    1910           0 :       str3(lo:up)        = c0
    1911           0 :       str4(lo:up)        = c0
    1912           0 :       str5(lo:up)        = c0
    1913           0 :       str6(lo:up)        = c0
    1914           0 :       str7(lo:up)        = c0
    1915           0 :       str8(lo:up)        = c0
    1916             : 
    1917           0 :    end subroutine numainit
    1918             : 
    1919             : !=======================================================================
    1920             : 
    1921             : end module ice_dyn_evp_1d

Generated by: LCOV version 1.14-6-g40580cd