New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14117 for NEMO – NEMO

Changeset 14117 for NEMO


Ignore:
Timestamp:
2020-12-07T11:21:42+01:00 (4 years ago)
Author:
acc
Message:

Re-working of icedyn_rhg_eap.F90 for better adhesion to coding standards and a reduction in the number of functions. AGRIF is no longer confused by this routine so the bypass introduced at revision 14020 has been removed. May need to force a clean remake of AGRIF-based SETTE tests to pick up changes. No changes in actual SETTE results.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_rhg_eap.F90

    r14072 r14117  
    1616   !!                                           CICE code (Tsamados, Heorton) 
    1717   !!---------------------------------------------------------------------- 
    18 #if defined key_si3 && ! defined key_agrif 
     18#if defined key_si3  
    1919   !!---------------------------------------------------------------------- 
    2020   !!   'key_si3'                                       SI3 sea-ice model 
     
    6666   INTEGER ::   ncvgid   ! netcdf file id 
    6767   INTEGER ::   nvarid   ! netcdf variable id 
    68    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
     68   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   aimsk00 
     69   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   eap_res  , aimsk15 
    6970   !!---------------------------------------------------------------------- 
    7071   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    202203      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 
    203204      ! 
    204       ! for diagnostics and convergence tests 
    205       ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 
     205      IF( kt == nit000 )  THEN  
     206         ! 
     207         ! for diagnostics  
     208         ALLOCATE( aimsk00(jpi,jpj) ) 
     209         ! for convergence tests 
     210         IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) ) 
     211      ENDIF 
     212      ! 
    206213      DO_2D( 1, 1, 1, 1 ) 
    207          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
    208          zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     214         aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
    209215      END_2D 
     216      IF( nn_rhg_chkcvg > 0 ) THEN 
     217         DO_2D( 1, 1, 1, 1 ) 
     218            aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     219         END_2D 
     220      ENDIF 
    210221      ! 
    211222!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
     
    349360            zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    350361            ! ice-bottom stress at U points 
    351             zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 
     362            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 
    352363            ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    353364            ! ice-bottom stress at V points 
    354             zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 
     365            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 
    355366            ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    356367            ! ice_bottom stress at T points 
    357             zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 
     368            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) )    ! if grounded icebergs are read: ocean depth = 0 
    358369            tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    359370         END_2D 
     
    819830            &                                  ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 
    820831         ! 
    821          CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
    822          CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 
    823          CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 
    824          CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 
    825          CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 
    826          CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 
     832         CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 ) 
     833         CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 ) 
     834         CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 ) 
     835         CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 ) 
     836         CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 ) 
     837         CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 ) 
    827838      ENDIF 
    828839 
    829840      ! --- divergence, shear and strength --- ! 
    830       IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    831       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
    832       IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    833       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
     841      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * aimsk00 )   ! divergence 
     842      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * aimsk00 )   ! shear 
     843      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * aimsk00 )   ! delta 
     844      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * aimsk00 )   ! strength 
    834845 
    835846      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     
    856867         ! 
    857868         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 
    858          IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 
    859          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
     869         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * aimsk00(:,:) ) ! Normal stress 
     870         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress 
    860871 
    861872         DEALLOCATE ( zsig_I, zsig_II ) 
     
    903914         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp ) 
    904915 
    905          CALL iom_put( 'yield11', zyield11 * zmsk00 ) 
    906          CALL iom_put( 'yield22', zyield22 * zmsk00 ) 
    907          CALL iom_put( 'yield12', zyield12 * zmsk00 ) 
     916         CALL iom_put( 'yield11', zyield11 * aimsk00 ) 
     917         CALL iom_put( 'yield22', zyield22 * aimsk00 ) 
     918         CALL iom_put( 'yield12', zyield12 * aimsk00 ) 
    908919      ENDIF 
    909920 
     
    911922      IF( iom_use('aniso') ) THEN 
    912923         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp ) 
    913          CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) 
     924         CALL iom_put( 'aniso' , paniso_11 * aimsk00 ) 
    914925      ENDIF 
    915926 
     
    922933            &                                    zfU, 'U', -1.0_wp,   zfV, 'V', -1.0_wp ) 
    923934 
    924          CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x) 
    925          CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y) 
    926          CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x) 
    927          CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y) 
    928          CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x) 
    929          CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y) 
     935         CALL iom_put( 'dssh_dx' , zspgU * aimsk00 )   ! Sea-surface tilt term in force balance (x) 
     936         CALL iom_put( 'dssh_dy' , zspgV * aimsk00 )   ! Sea-surface tilt term in force balance (y) 
     937         CALL iom_put( 'corstrx' , zCorU * aimsk00 )   ! Coriolis force term in force balance (x) 
     938         CALL iom_put( 'corstry' , zCorV * aimsk00 )   ! Coriolis force term in force balance (y) 
     939         CALL iom_put( 'intstrx' , zfU   * aimsk00 )   ! Internal force term in force balance (x) 
     940         CALL iom_put( 'intstry' , zfV   * aimsk00 )   ! Internal force term in force balance (y) 
    930941      ENDIF 
    931942 
     
    938949         DO_2D( 0, 0, 0, 0 ) 
    939950            ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    940             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
    941             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 
     951            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj) 
     952            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj) 
    942953 
    943954            zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
     
    973984            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
    974985               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 
    975                   &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     986                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) ) 
    976987            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
    977988               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 
    978                   &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     989                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) ) 
    979990            ENDIF 
    980991         ENDIF 
    981992      ENDIF 
    982       ! 
    983       DEALLOCATE( zmsk00, zmsk15 ) 
    984993      ! 
    985994   END SUBROUTINE ice_dyn_rhg_eap 
     
    10061015      REAL(wp)          ::   zresm           ! local real 
    10071016      CHARACTER(len=20) ::   clname 
    1008       REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence 
    10091017      !!---------------------------------------------------------------------- 
    10101018 
     
    10371045      ELSE 
    10381046         DO_2D( 1, 1, 1, 1 ) 
    1039             zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
    1040                &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 
     1047            eap_res(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1048               &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * aimsk15(ji,jj) 
    10411049         END_2D 
    1042          zresm = MAXVAL( zres ) 
     1050 
     1051         zresm = MAXVAL( eap_res ) 
    10431052         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
    10441053      ENDIF 
     
    10801089      REAL(wp) ::   zsig11, zsig12, zsig22 
    10811090      REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22 
    1082       REAL(wp) ::   zinvstressconviso 
    10831091      REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha 
    10841092      REAL(wp) ::   zTany_1, zTany_2 
    1085       REAL(wp) ::   zx, zy, zdx, zdy, zda, zkxw, kyw, kaw 
     1093      REAL(wp) ::   zx, zy, zkxw, kyw, kaw 
    10861094      REAL(wp) ::   zinvdx, zinvdy, zinvda 
    1087       REAL(wp) ::   zdtemp1, zdtemp2, zatempprime, zinvsin 
    1088  
    1089       REAL(wp), PARAMETER ::   kfriction = 0.45_wp 
    1090       !!--------------------------------------------------------------------- 
     1095      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime 
     1096 
     1097      REAL(wp), PARAMETER ::   ppkfriction = 0.45_wp 
    10911098      ! Factor to maintain the same stress as in EVP (see Section 3) 
    10921099      ! Can be set to 1 otherwise 
    1093 !      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
    1094       zinvstressconviso = 1._wp 
    1095  
    1096       zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso 
    1097       !now uses phi as set in higher code 
     1100!     REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp/(1._wp+ppkfriction*ppkfriction) 
     1101      REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp 
     1102 
     1103      ! next statement uses pphi set in main module (icedyn_rhg_eap) 
     1104      REAL(wp), PARAMETER ::   ppinvsin = 1._wp/sin(2._wp*pphi) * ppinvstressconviso 
    10981105 
    10991106      ! compute eigenvalues, eigenvectors and angles for structure tensor, strain 
     
    11751182 
    11761183      ! 3) update anisotropic stress tensor 
    1177       zdx   = rpi/real(nx_yield-1,kind=wp) 
    1178       zdy   = rpi/real(ny_yield-1,kind=wp) 
    1179       zda   = 0.5_wp/real(na_yield-1,kind=wp) 
    1180       zinvdx = 1._wp/zdx 
    1181       zinvdy = 1._wp/zdy 
    1182       zinvda = 1._wp/zda 
     1184      zinvdx = real(nx_yield-1,kind=wp)/rpi 
     1185      zinvdy = real(ny_yield-1,kind=wp)/rpi 
     1186      zinvda = 2._wp*real(na_yield-1,kind=wp) 
    11831187 
    11841188      ! % need 8 coords and 8 weights 
     
    12581262      ! Tsamados 2013) 
    12591263 
    1260       zsig11  = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin 
    1261       zsig12  = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin 
    1262       zsig22  = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin 
     1264      zsig11  = pstrength*(zstemp11r + ppkfriction*zstemp11s) * ppinvsin 
     1265      zsig12  = pstrength*(zstemp12r + ppkfriction*zstemp12s) * ppinvsin 
     1266      zsig22  = pstrength*(zstemp22r + ppkfriction*zstemp22s) * ppinvsin 
    12631267 
    12641268      ! Back - rotation of the stress from principal axes into general coordinates 
     
    13191323      REAL (wp) ::   zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 
    13201324 
    1321 !!$      REAL (wp), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation 
    1322       REAL (wp), PARAMETER ::   kfrac = 1.e-3_wp   ! rate of fracture formation 
    1323       REAL (wp), PARAMETER ::   threshold = 0.3_wp  ! critical confinement ratio 
     1325!!$   REAL (wp), PARAMETER ::   ppkfrac = 0.0001_wp   ! rate of fracture formation 
     1326      REAL (wp), PARAMETER ::   ppkfrac = 1.e-3_wp    ! rate of fracture formation 
     1327      REAL (wp), PARAMETER ::   ppthreshold = 0.3_wp  ! critical confinement ratio 
    13241328      !!--------------------------------------------------------------- 
    13251329      ! 
     
    13631367      ! which leads to the loss of their shape, so we again model it through diffusion 
    13641368      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN 
    1365          pmresult11 = - kfrac * (pa11 - zQ12Q12) 
    1366          pmresult12 = - kfrac * (pa12 + zQ11Q12) 
     1369         pmresult11 = - ppkfrac * (pa11 - zQ12Q12) 
     1370         pmresult12 = - ppkfrac * (pa12 + zQ11Q12) 
    13671371 
    13681372      ! Shear faulting 
     
    13701374         pmresult11 = 0.0_wp 
    13711375         pmresult12 = 0.0_wp 
    1372       ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 
    1373          pmresult11 = - kfrac * (pa11 - zQ12Q12) 
    1374          pmresult12 = - kfrac * (pa12 + zQ11Q12) 
     1376      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= ppthreshold)) THEN 
     1377         pmresult11 = - ppkfrac * (pa11 - zQ12Q12) 
     1378         pmresult12 = - ppkfrac * (pa12 + zQ11Q12) 
    13751379 
    13761380      ! Horizontal spalling 
     
    14051409      !!clem 
    14061410      REAL(wp) ::   zw1, zw2, zfac, ztemp 
    1407       REAL(wp) ::   idx, idy, idz 
     1411      REAL(wp) ::   zidx, zidy, zidz 
     1412      REAL(wp) ::   zsaak(6)                  ! temporary array 
    14081413 
    14091414      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp 
     
    15221527            zw2 = w2(ainit+ia*da) 
    15231528            DO iz = 1, nz 
    1524                idz = zinit+iz*dz 
     1529               zidz = zinit+iz*dz 
    15251530               ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz)) 
    15261531               DO iy = 1, ny_yield 
    1527                   idy = yinit+iy*dy 
     1532                  zidy = yinit+iy*dy 
    15281533                  DO ix = 1, nx_yield 
    1529                      idx = xinit+ix*dx 
    1530                      s11r(ix,iy,ia) = s11r(ix,iy,ia) + ztemp * s11kr(idx,idy,idz)*zfac 
    1531                      s12r(ix,iy,ia) = s12r(ix,iy,ia) + ztemp * s12kr(idx,idy,idz)*zfac 
    1532                      s22r(ix,iy,ia) = s22r(ix,iy,ia) + ztemp * s22kr(idx,idy,idz)*zfac 
    1533                      s11s(ix,iy,ia) = s11s(ix,iy,ia) + ztemp * s11ks(idx,idy,idz)*zfac 
    1534                      s12s(ix,iy,ia) = s12s(ix,iy,ia) + ztemp * s12ks(idx,idy,idz)*zfac 
    1535                      s22s(ix,iy,ia) = s22s(ix,iy,ia) + ztemp * s22ks(idx,idy,idz)*zfac 
     1534                     zidx = xinit+ix*dx 
     1535                     call all_skr_sks(zidx,zidy,zidz,zsaak) 
     1536                     zsaak = ztemp*zsaak*zfac 
     1537                     s11r(ix,iy,ia) = s11r(ix,iy,ia) + zsaak(1) 
     1538                     s12r(ix,iy,ia) = s12r(ix,iy,ia) + zsaak(2) 
     1539                     s22r(ix,iy,ia) = s22r(ix,iy,ia) + zsaak(3) 
     1540                     s11s(ix,iy,ia) = s11s(ix,iy,ia) + zsaak(4) 
     1541                     s12s(ix,iy,ia) = s12s(ix,iy,ia) + zsaak(5) 
     1542                     s22s(ix,iy,ia) = s22s(ix,iy,ia) + zsaak(6) 
    15361543                  END DO 
    15371544               END DO 
    15381545            END DO 
    15391546         END DO 
    1540  
    15411547         zfac = 1._wp/sin(2._wp*pphi) 
    15421548         ia = na_yield 
    15431549         DO iy = 1, ny_yield 
    1544             idy = yinit+iy*dy 
     1550            zidy = yinit+iy*dy 
    15451551            DO ix = 1, nx_yield 
    1546                idx = xinit+ix*dx 
    1547                s11r(ix,iy,ia) = 0.5_wp*s11kr(idx,idy,0._wp)*zfac 
    1548                s12r(ix,iy,ia) = 0.5_wp*s12kr(idx,idy,0._wp)*zfac 
    1549                s22r(ix,iy,ia) = 0.5_wp*s22kr(idx,idy,0._wp)*zfac 
    1550                s11s(ix,iy,ia) = 0.5_wp*s11ks(idx,idy,0._wp)*zfac 
    1551                s12s(ix,iy,ia) = 0.5_wp*s12ks(idx,idy,0._wp)*zfac 
    1552                s22s(ix,iy,ia) = 0.5_wp*s22ks(idx,idy,0._wp)*zfac 
     1552               zidx = xinit+ix*dx 
     1553               call all_skr_sks(zidx,zidy,0._wp,zsaak) 
     1554               zsaak = 0.5_wp*zsaak*zfac 
     1555               s11r(ix,iy,ia) = zsaak(1) 
     1556               s12r(ix,iy,ia) = zsaak(2) 
     1557               s22r(ix,iy,ia) = zsaak(3) 
     1558               s11s(ix,iy,ia) = zsaak(4) 
     1559               s12s(ix,iy,ia) = zsaak(5) 
     1560               s22s(ix,iy,ia) = zsaak(6) 
    15531561            ENDDO 
    15541562         ENDDO 
     
    16161624   END FUNCTION w2 
    16171625 
    1618    FUNCTION s11kr(px,py,pz) 
    1619       !!------------------------------------------------------------------- 
    1620       !! Function : s11kr 
    1621       !!------------------------------------------------------------------- 
     1626   SUBROUTINE all_skr_sks( px, py, pz, allsk ) 
    16221627      REAL(wp), INTENT(in   ) ::   px,py,pz 
    1623  
    1624       REAL(wp) ::   s11kr, zpih 
    1625  
     1628      REAL(wp), INTENT(out  ) ::   allsk(6) 
     1629 
     1630      REAL(wp) ::   zs12r0, zs21r0 
     1631      REAL(wp) ::   zs12s0, zs21s0 
     1632 
     1633      REAL(wp) ::   zpih 
    16261634      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    16271635      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     
    16731681      ENDIF 
    16741682 
    1675       s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 
    1676  
    1677    END FUNCTION s11kr 
    1678  
    1679    FUNCTION s12kr(px,py,pz) 
     1683      !!------------------------------------------------------------------- 
     1684      !! Function : s11kr 
     1685      !!------------------------------------------------------------------- 
     1686      allsk(1) = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 
    16801687      !!------------------------------------------------------------------- 
    16811688      !! Function : s12kr 
    16821689      !!------------------------------------------------------------------- 
    1683       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1684  
    1685       REAL(wp) ::   s12kr, zs12r0, zs21r0, zpih 
    1686  
    1687       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1688       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1689       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1690       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1691       REAL(wp) ::   zd11, zd12, zd22 
    1692       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1693       REAL(wp) ::   zHen1t2, zHen2t1 
    1694       !!------------------------------------------------------------------- 
    1695       zpih = 0.5_wp*rpi 
    1696  
    1697       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1698       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1699       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1700       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1701       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1702       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1703       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1704       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1705       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1706       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1707       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1708       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1709       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1710       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1711       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1712       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1713       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1714       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1715       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1716       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1717       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1718       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1719  
    1720       IF (-zIIn1t2>=rsmall) THEN 
    1721       zHen1t2 = 1._wp 
    1722       ELSE 
    1723       zHen1t2 = 0._wp 
    1724       ENDIF 
    1725  
    1726       IF (-zIIn2t1>=rsmall) THEN 
    1727       zHen2t1 = 1._wp 
    1728       ELSE 
    1729       zHen2t1 = 0._wp 
    1730       ENDIF 
    1731  
    17321690      zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 
    17331691      zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21) 
    1734       s12kr=0.5_wp*(zs12r0+zs21r0) 
    1735  
    1736    END FUNCTION s12kr 
    1737  
    1738    FUNCTION s22kr(px,py,pz) 
     1692      allsk(2)=0.5_wp*(zs12r0+zs21r0) 
    17391693      !!------------------------------------------------------------------- 
    17401694      !! Function : s22kr 
    17411695      !!------------------------------------------------------------------- 
    1742       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1743  
    1744       REAL(wp) ::   s22kr, zpih 
    1745  
    1746       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1747       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1748       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1749       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1750       REAL(wp) ::   zd11, zd12, zd22 
    1751       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1752       REAL(wp) ::   zHen1t2, zHen2t1 
    1753       !!------------------------------------------------------------------- 
    1754       zpih = 0.5_wp*rpi 
    1755  
    1756       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1757       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1758       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1759       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1760       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1761       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1762       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1763       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1764       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1765       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1766       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1767       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1768       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1769       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1770       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1771       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1772       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1773       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1774       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1775       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1776       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1777       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1778  
    1779       IF (-zIIn1t2>=rsmall) THEN 
    1780       zHen1t2 = 1._wp 
    1781       ELSE 
    1782       zHen1t2 = 0._wp 
    1783       ENDIF 
    1784  
    1785       IF (-zIIn2t1>=rsmall) THEN 
    1786       zHen2t1 = 1._wp 
    1787       ELSE 
    1788       zHen2t1 = 0._wp 
    1789       ENDIF 
    1790  
    1791       s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 
    1792  
    1793    END FUNCTION s22kr 
    1794  
    1795    FUNCTION s11ks(px,py,pz) 
     1696      allsk(3) = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 
    17961697      !!------------------------------------------------------------------- 
    17971698      !! Function : s11ks 
    17981699      !!------------------------------------------------------------------- 
    1799       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1800  
    1801       REAL(wp) ::   s11ks, zpih 
    1802  
    1803       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1804       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1805       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1806       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1807       REAL(wp) ::   zd11, zd12, zd22 
    1808       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1809       REAL(wp) ::   zHen1t2, zHen2t1 
    1810       !!------------------------------------------------------------------- 
    1811       zpih = 0.5_wp*rpi 
    1812  
    1813       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1814       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1815       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1816       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1817       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1818       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1819       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1820       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1821       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1822       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1823       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1824       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1825       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1826       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1827       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1828       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1829       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1830       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1831       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1832       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1833       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1834       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1835  
    1836       IF (-zIIn1t2>=rsmall) THEN 
    1837       zHen1t2 = 1._wp 
    1838       ELSE 
    1839       zHen1t2 = 0._wp 
    1840       ENDIF 
    1841  
    1842       IF (-zIIn2t1>=rsmall) THEN 
    1843       zHen2t1 = 1._wp 
    1844       ELSE 
    1845       zHen2t1 = 0._wp 
    1846       ENDIF 
    1847  
    1848       s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 
    1849  
    1850    END FUNCTION s11ks 
    1851  
    1852    FUNCTION s12ks(px,py,pz) 
     1700 
     1701      allsk(4) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 
    18531702      !!------------------------------------------------------------------- 
    18541703      !! Function : s12ks 
    18551704      !!------------------------------------------------------------------- 
    1856       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1857  
    1858       REAL(wp) ::   s12ks, zs12s0, zs21s0, zpih 
    1859  
    1860       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1861       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1862       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1863       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1864       REAL(wp) ::   zd11, zd12, zd22 
    1865       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1866       REAL(wp) ::   zHen1t2, zHen2t1 
    1867       !!------------------------------------------------------------------- 
    1868       zpih = 0.5_wp*rpi 
    1869  
    1870       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1871       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1872       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1873       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1874       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1875       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1876       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1877       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1878       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1879       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1880       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1881       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1882       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1883       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1884       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1885       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1886       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1887       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1888       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1889       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1890       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1891       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1892  
    1893       IF (-zIIn1t2>=rsmall) THEN 
    1894       zHen1t2 = 1._wp 
    1895       ELSE 
    1896       zHen1t2 = 0._wp 
    1897       ENDIF 
    1898  
    1899       IF (-zIIn2t1>=rsmall) THEN 
    1900       zHen2t1 = 1._wp 
    1901       ELSE 
    1902       zHen2t1 = 0._wp 
    1903       ENDIF 
    1904  
    19051705      zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 
    19061706      zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21) 
    1907       s12ks=0.5_wp*(zs12s0+zs21s0) 
    1908  
    1909    END FUNCTION s12ks 
    1910  
    1911    FUNCTION s22ks(px,py,pz) 
     1707      allsk(5)=0.5_wp*(zs12s0+zs21s0) 
    19121708      !!------------------------------------------------------------------- 
    19131709      !! Function : s22ks 
    19141710      !!------------------------------------------------------------------- 
    1915       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1916  
    1917       REAL(wp) ::   s22ks, zpih 
    1918  
    1919       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1920       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1921       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1922       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1923       REAL(wp) ::   zd11, zd12, zd22 
    1924       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1925       REAL(wp) ::   zHen1t2, zHen2t1 
    1926       !!------------------------------------------------------------------- 
    1927       zpih = 0.5_wp*rpi 
    1928  
    1929       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1930       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1931       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1932       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1933       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1934       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1935       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1936       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1937       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1938       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1939       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1940       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1941       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1942       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1943       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1944       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1945       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1946       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1947       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1948       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1949       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1950       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1951  
    1952       IF (-zIIn1t2>=rsmall) THEN 
    1953       zHen1t2 = 1._wp 
    1954       ELSE 
    1955       zHen1t2 = 0._wp 
    1956       ENDIF 
    1957  
    1958       IF (-zIIn2t1>=rsmall) THEN 
    1959       zHen2t1 = 1._wp 
    1960       ELSE 
    1961       zHen2t1 = 0._wp 
    1962       ENDIF 
    1963  
    1964       s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 
    1965  
    1966    END FUNCTION s22ks 
     1711      allsk(6) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 
     1712   END SUBROUTINE all_skr_sks 
    19671713 
    19681714#else 
  • NEMO/trunk/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90

    r14021 r14117  
    66   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code 
    77   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model 
    8    !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy  
     8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
    99   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case 
    10    !!            3.4  !  2011-01  (A. Porter)   dynamical allocation  
     10   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation 
    1111   !!            3.5  !  2012-08  (R. Benshila) AGRIF 
    1212   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + mEVP (Bouillon 2013) 
     
    1414   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube] 
    1515   !!                 !  2019     (S. Rynders, Y. Aksenov, C. Rousset)  change into eap rheology from 
    16    !!                                           CICE code (Tsamados, Heorton)  
     16   !!                                           CICE code (Tsamados, Heorton) 
    1717   !!---------------------------------------------------------------------- 
    1818#if defined key_si3 
     
    3030   USE icevar         ! ice_var_sshdyn 
    3131   USE icedyn_rdgrft  ! sea-ice: ice strength 
    32    USE bdy_oce , ONLY : ln_bdy  
    33    USE bdyice  
     32   USE bdy_oce , ONLY : ln_bdy 
     33   USE bdyice 
    3434#if defined key_agrif 
    3535   USE agrif_ice_interp 
     
    6666   INTEGER ::   ncvgid   ! netcdf file id 
    6767   INTEGER ::   nvarid   ! netcdf variable id 
    68    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
     68   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   aimsk00 
     69   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   eap_res  , aimsk15 
    6970   !!---------------------------------------------------------------------- 
    7071   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    8081      !! 
    8182      !! ** purpose : determines sea ice drift from wind stress, ice-ocean 
    82       !!  stress and sea-surface slope. Ice-ice interaction is described by  
    83       !!  a non-linear elasto-anisotropic-plastic (EAP) law including shear  
    84       !!  strength and a bulk rheology .   
     83      !!  stress and sea-surface slope. Ice-ice interaction is described by 
     84      !!  a non-linear elasto-anisotropic-plastic (EAP) law including shear 
     85      !!  strength and a bulk rheology . 
    8586      !! 
    8687      !!  The points in the C-grid look like this, dear reader 
     
    9091      !!                                 | 
    9192      !!                      (ji-1,jj)  |  (ji,jj) 
    92       !!                             ---------    
     93      !!                             --------- 
    9394      !!                            |         | 
    9495      !!                            | (ji,jj) |------(ji,jj) 
    9596      !!                            |         | 
    96       !!                             ---------    
     97      !!                             --------- 
    9798      !!                     (ji-1,jj-1)     (ji,jj-1) 
    9899      !! 
     
    101102      !!                snow total volume (vt_s) per unit area 
    102103      !! 
    103       !! ** Action  : - compute u_ice, v_ice : the components of the  
     104      !! ** Action  : - compute u_ice, v_ice : the components of the 
    104105      !!                sea-ice velocity vector 
    105106      !!              - compute delta_i, shear_i, divu_i, which are inputs 
     
    107108      !! 
    108109      !! ** Steps   : 0) compute mask at F point 
    109       !!              1) Compute ice snow mass, ice strength  
     110      !!              1) Compute ice snow mass, ice strength 
    110111      !!              2) Compute wind, oceanic stresses, mass terms and 
    111112      !!                 coriolis terms of the momentum equation 
     
    166167      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    167168      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    168       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     169      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points 
    169170      ! 
    170171      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
     
    172173      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    173174      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
    174       !                                                                 !    ice bottom surface if ice is embedded    
     175      !                                                                 !    ice bottom surface if ice is embedded 
    175176      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses 
    176177      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points 
     
    192193      !! --- diags 
    193194      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength 
    194       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p          
     195      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p 
    195196      !! --- SIMIP diags 
    196197      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
     
    199200      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 
    200201      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s) 
    201       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)       
     202      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s) 
    202203      !!------------------------------------------------------------------- 
    203204 
    204205      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 
    205206      ! 
    206       ! for diagnostics and convergence tests 
    207       ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 
     207      IF( kt == nit000 )  THEN  
     208         ! 
     209         ! for diagnostics  
     210         ALLOCATE( aimsk00(jpi,jpj) ) 
     211         ! for convergence tests 
     212         IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) ) 
     213      ENDIF 
     214      ! 
    208215      DO_2D( 1, 1, 1, 1 ) 
    209          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
    210          zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     216         aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
    211217      END_2D 
     218      IF( nn_rhg_chkcvg > 0 ) THEN 
     219         DO_2D( 1, 1, 1, 1 ) 
     220            aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     221         END_2D 
     222      ENDIF 
    212223      ! 
    213224!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
     
    249260      ! 1) define some variables and initialize arrays 
    250261      !------------------------------------------------------------------------------! 
    251       zrhoco = rho0 * rn_cio  
     262      zrhoco = rho0 * rn_cio 
    252263!extra code for test case boundary conditions 
    253264      zinvw=1._wp/(zrhoco*0.5_wp) 
     
    270281      ENDIF 
    271282      z1_dtevp = 1._wp / zdtevp 
    272           
    273       ! Initialise stress tensor  
    274       zs1 (:,:) = pstress1_i (:,:)  
     283 
     284      ! Initialise stress tensor 
     285      zs1 (:,:) = pstress1_i (:,:) 
    275286      zs2 (:,:) = pstress2_i (:,:) 
    276287      zs12(:,:) = pstress12_i(:,:) 
     
    319330         ! dt/m at T points (for alpha and beta coefficients) 
    320331         zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin ) 
    321           
     332 
    322333         ! m/dt 
    323334         zmU_t(ji,jj)    = zmassU * z1_dtevp 
    324335         zmV_t(ji,jj)    = zmassV * z1_dtevp 
    325           
     336 
    326337         ! Drag ice-atm. 
    327338         ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     
    353364            zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    354365            ! ice-bottom stress at U points 
    355             zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 
     366            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 
    356367            ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    357368            ! ice-bottom stress at V points 
    358             zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 
     369            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 
    359370            ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    360371            ! ice_bottom stress at T points 
    361             zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 
     372            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) )    ! if grounded icebergs are read: ocean depth = 0 
    362373            tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    363374         END_2D 
     
    377388      !                                               ! ==================== ! 
    378389      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    379          !                                            ! ==================== !         
     390         !                                            ! ==================== ! 
    380391         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    381392         ! 
     
    404415               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    405416               &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
    406             
     417 
    407418            ! divergence at T points 
    408419            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     
    410421               &    ) * r1_e1e2t(ji,jj) 
    411422            zdiv2 = zdiv * zdiv 
    412              
     423 
    413424            ! tension at T points 
    414425            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     
    418429 
    419430            ! delta at T points 
    420             zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     431            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
    421432 
    422433         END_2D 
    423434         CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp ) 
    424                 
     435 
    425436         ! P/delta at T points 
    426437         DO_2D( 1, 1, 1, 1 ) 
     
    430441         DO_2D( 0, 1, 0, 1 )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 
    431442 
    432              ! shear at T points  
     443             ! shear at T points 
    433444            zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
    434445               &   + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    435446               &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
    436              
     447 
    437448           ! divergence at T points (duplication to avoid communications) 
    438449            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    439450               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    440451               &    ) * r1_e1e2t(ji,jj) 
    441              
     452 
    442453            ! tension at T points (duplication to avoid communications) 
    443454            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     
    459470               zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp) 
    460471               zyield12(ji,jj) = zstress12tmp(ji,jj) 
    461                prdg_conv(ji,jj) = -min( zalphar, 0._wp)     
     472               prdg_conv(ji,jj) = -min( zalphar, 0._wp) 
    462473            ENDIF 
    463474 
     
    491502 
    492503         DO_2D( 1, 0, 1, 0 ) 
    493             ! stress12tmp at F points  
     504            ! stress12tmp at F points 
    494505            zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
    495506               &            + zstress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + zstress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  & 
     
    504515               ! zalph2 = zalph2 - 1._wp 
    505516            ENDIF 
    506              
     517 
    507518            ! stress at F points (zkt/=0 if landfast) 
    508519            zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1 
     
    570581                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    571582                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    572                      &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &  
     583                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
    573584                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
    574585                     &                                    ) / ( zbetav + 1._wp )                                              & 
     
    630641                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
    631642                     &                                    ) / ( zbetau + 1._wp )                                              & 
    632                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     643                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    633644                     &           )   * zmsk00x(ji,jj) 
    634645               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
     
    637648                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
    638649                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
    639                      &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                  & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     650                     &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                  & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    640651                     &            )   * zmsk00x(ji,jj) 
    641652               ENDIF 
     
    689700                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
    690701                     &                                    ) / ( zbetau + 1._wp )                                              & 
    691                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     702                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    692703                     &           )   * zmsk00x(ji,jj) 
    693704               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
     
    744755                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
    745756                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
    746                      &                                    ) / ( zbetav + 1._wp )                                              &  
     757                     &                                    ) / ( zbetav + 1._wp )                                              & 
    747758                     &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    748759                     &           )   * zmsk00y(ji,jj) 
     
    781792      ! 
    782793      !------------------------------------------------------------------------------! 
    783       ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)  
     794      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
    784795      !------------------------------------------------------------------------------! 
    785796      DO_2D( 1, 0, 1, 0 ) 
     
    791802 
    792803      END_2D 
    793        
     804 
    794805      DO_2D( 0, 0, 0, 0 ) 
    795           
     806 
    796807         ! tension**2 at T points 
    797808         zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     
    799810            &   ) * r1_e1e2t(ji,jj) 
    800811         zdt2 = zdt * zdt 
    801           
     812 
    802813         zten_i(ji,jj) = zdt 
    803814 
     
    806817            &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    807818            &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
    808           
     819 
    809820         ! shear at T points 
    810821         pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     
    814825            &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    815826            &             ) * r1_e1e2t(ji,jj) 
    816           
     827 
    817828         ! delta at T points 
    818          zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta   
     829         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 
    819830         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 
    820831         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 
     
    824835         &                                    zten_i, 'T', 1.0_wp, zs1    , 'T', 1.0_wp, zs2     , 'T', 1.0_wp, & 
    825836         &                                      zs12, 'F', 1.0_wp ) 
    826        
     837 
    827838      ! --- Store the stress tensor for the next time step --- ! 
    828839      pstress1_i (:,:) = zs1 (:,:) 
     
    841852            &                                  ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 
    842853         ! 
    843          CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
    844          CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 
    845          CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 
    846          CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 
    847          CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 
    848          CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 
    849       ENDIF 
    850         
     854         CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 ) 
     855         CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 ) 
     856         CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 ) 
     857         CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 ) 
     858         CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 ) 
     859         CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 ) 
     860      ENDIF 
     861 
    851862      ! --- divergence, shear and strength --- ! 
    852       IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    853       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
    854       IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    855       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
     863      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * aimsk00 )   ! divergence 
     864      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * aimsk00 )   ! shear 
     865      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * aimsk00 )   ! delta 
     866      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * aimsk00 )   ! strength 
    856867 
    857868      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     
    859870         ! 
    860871         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    861          !          
     872         ! 
    862873         DO_2D( 1, 1, 1, 1 ) 
    863           
     874 
    864875            ! Ice stresses 
    865876            ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 
    866877            ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 
    867878            ! I know, this can be confusing... 
    868             zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )  
     879            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 
    869880            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
    870881            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
    871882            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
    872              
     883 
    873884            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 
    874885            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
    875886            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
    876              
     887 
    877888         END_2D 
    878889         ! 
    879890         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 
    880          IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 
    881          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
    882           
     891         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * aimsk00(:,:) ) ! Normal stress 
     892         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress 
     893 
    883894         DEALLOCATE ( zsig_I, zsig_II ) 
    884           
     895 
    885896      ENDIF 
    886897 
     
    891902      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
    892903         ! 
    893          ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )          
    894          !          
     904         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
     905         ! 
    895906         DO_2D( 1, 1, 1, 1 ) 
    896           
    897             ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     907 
     908            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
    898909            !                        and **deformations** at current iterates 
    899910            !                        following Lemieux & Dupont (2020) 
     
    902913            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
    903914            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
    904              
     915 
    905916            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
    906917            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
    907918            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
    908     
     919 
    909920            ! Normalized  principal stresses (used to display the ellipse) 
    910921            z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) ) 
     
    913924         END_2D 
    914925         ! 
    915          CALL iom_put( 'sig1_pnorm' , zsig1_p )  
    916          CALL iom_put( 'sig2_pnorm' , zsig2_p )  
    917        
     926         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
     927         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
     928 
    918929         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
    919           
     930 
    920931      ENDIF 
    921932 
     
    925936         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp ) 
    926937 
    927          CALL iom_put( 'yield11', zyield11 * zmsk00 ) 
    928          CALL iom_put( 'yield22', zyield22 * zmsk00 ) 
    929          CALL iom_put( 'yield12', zyield12 * zmsk00 ) 
     938         CALL iom_put( 'yield11', zyield11 * aimsk00 ) 
     939         CALL iom_put( 'yield22', zyield22 * aimsk00 ) 
     940         CALL iom_put( 'yield12', zyield12 * aimsk00 ) 
    930941      ENDIF 
    931942 
    932943      ! --- anisotropy tensor --- ! 
    933       IF( iom_use('aniso') ) THEN        
    934          CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp )  
    935          CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) 
    936       ENDIF 
    937        
     944      IF( iom_use('aniso') ) THEN 
     945         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp ) 
     946         CALL iom_put( 'aniso' , paniso_11 * aimsk00 ) 
     947      ENDIF 
     948 
    938949      ! --- SIMIP --- ! 
    939950      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
     
    944955            &                                    zfU, 'U', -1.0_wp,   zfV, 'V', -1.0_wp ) 
    945956 
    946          CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x) 
    947          CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y) 
    948          CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x) 
    949          CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y) 
    950          CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x) 
    951          CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y) 
     957         CALL iom_put( 'dssh_dx' , zspgU * aimsk00 )   ! Sea-surface tilt term in force balance (x) 
     958         CALL iom_put( 'dssh_dy' , zspgV * aimsk00 )   ! Sea-surface tilt term in force balance (y) 
     959         CALL iom_put( 'corstrx' , zCorU * aimsk00 )   ! Coriolis force term in force balance (x) 
     960         CALL iom_put( 'corstry' , zCorV * aimsk00 )   ! Coriolis force term in force balance (y) 
     961         CALL iom_put( 'intstrx' , zfU   * aimsk00 )   ! Internal force term in force balance (x) 
     962         CALL iom_put( 'intstry' , zfV   * aimsk00 )   ! Internal force term in force balance (y) 
    952963      ENDIF 
    953964 
     
    960971         DO_2D( 0, 0, 0, 0 ) 
    961972            ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    962             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
    963             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 
     973            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj) 
     974            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj) 
    964975 
    965976            zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
     
    979990 
    980991         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s) 
    981          CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport  
     992         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport 
    982993         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s) 
    983994         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport 
     
    9951006            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
    9961007               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 
    997                   &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     1008                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) ) 
    9981009            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
    9991010               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 
    1000                   &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     1011                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) ) 
    10011012            ENDIF 
    10021013         ENDIF 
    1003       ENDIF       
    1004       ! 
    1005       DEALLOCATE( zmsk00, zmsk15 ) 
     1014      ENDIF 
    10061015      ! 
    10071016   END SUBROUTINE ice_dyn_rhg_eap 
    10081017 
    1009     
     1018 
    10101019   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 
    10111020      !!---------------------------------------------------------------------- 
    10121021      !!                    ***  ROUTINE rhg_cvg  *** 
    1013       !!                      
     1022      !! 
    10141023      !! ** Purpose :   check convergence of oce rheology 
    10151024      !! 
     
    10191028      !!                This routine is called every sub-iteration, so it is cpu expensive 
    10201029      !! 
    1021       !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     1030      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 
    10221031      !!---------------------------------------------------------------------- 
    10231032      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index 
     
    10261035      INTEGER           ::   it, idtime, istatus 
    10271036      INTEGER           ::   ji, jj          ! dummy loop indices 
    1028       REAL(wp)          ::   zresm           ! local real  
     1037      REAL(wp)          ::   zresm           ! local real 
    10291038      CHARACTER(len=20) ::   clname 
    1030       REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence 
    10311039      !!---------------------------------------------------------------------- 
    10321040 
     
    10531061      ! time 
    10541062      it = ( kt - 1 ) * kitermax + kiter 
    1055        
     1063 
    10561064      ! convergence 
    10571065      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 
     
    10591067      ELSE 
    10601068         DO_2D( 1, 1, 1, 1 ) 
    1061             zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
    1062                &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 
     1069            eap_res(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1070               &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * aimsk15(ji,jj) 
     1071            ! cut of the boundary of the box (forced velocities) 
     1072            IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN 
     1073               eap_res(ji,jj) = 0._wp 
     1074            END IF 
    10631075         END_2D 
    10641076 
    1065          ! cut of the boundary of the box (forced velocities) 
    1066          IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN 
    1067             zres(ji,jj) = 0._wp 
    1068          END IF 
    1069          zresm = MAXVAL( zres ) 
     1077         zresm = MAXVAL( eap_res ) 
    10701078         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
    10711079      ENDIF 
     
    10771085         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid) 
    10781086      ENDIF 
    1079        
     1087 
    10801088   END SUBROUTINE rhg_cvg 
    10811089 
     
    10851093      !!--------------------------------------------------------------------- 
    10861094      !!                   ***  SUBROUTINE update_stress_rdg  *** 
    1087       !!                      
     1095      !! 
    10881096      !! ** Purpose :   Updates the stress depending on values of strain rate and structure 
    10891097      !!                tensor and for the last subcycle step it computes closing and sliding rate 
     
    10981106      INTEGER  ::   kx ,ky, ka 
    10991107 
    1100       REAL(wp) ::   zstemp11r, zstemp12r, zstemp22r    
    1101       REAL(wp) ::   zstemp11s, zstemp12s, zstemp22s  
    1102       REAL(wp) ::   za22, zQd11Qd11, zQd11Qd12, zQd12Qd12  
    1103       REAL(wp) ::   zQ11Q11, zQ11Q12, zQ12Q12  
    1104       REAL(wp) ::   zdtemp11, zdtemp12, zdtemp22  
    1105       REAL(wp) ::   zrotstemp11r, zrotstemp12r, zrotstemp22r    
     1108      REAL(wp) ::   zstemp11r, zstemp12r, zstemp22r 
     1109      REAL(wp) ::   zstemp11s, zstemp12s, zstemp22s 
     1110      REAL(wp) ::   za22, zQd11Qd11, zQd11Qd12, zQd12Qd12 
     1111      REAL(wp) ::   zQ11Q11, zQ11Q12, zQ12Q12 
     1112      REAL(wp) ::   zdtemp11, zdtemp12, zdtemp22 
     1113      REAL(wp) ::   zrotstemp11r, zrotstemp12r, zrotstemp22r 
    11061114      REAL(wp) ::   zrotstemp11s, zrotstemp12s, zrotstemp22s 
    1107       REAL(wp) ::   zsig11, zsig12, zsig22  
    1108       REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22  
    1109       REAL(wp) ::   zinvstressconviso  
    1110       REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha  
    1111       REAL(wp) ::   zTany_1, zTany_2  
    1112       REAL(wp) ::   zx, zy, zdx, zdy, zda, zkxw, kyw, kaw 
    1113       REAL(wp) ::   zinvdx, zinvdy, zinvda    
    1114       REAL(wp) ::   zdtemp1, zdtemp2, zatempprime, zinvsin 
    1115  
    1116       REAL(wp), PARAMETER ::   kfriction = 0.45_wp 
    1117       !!--------------------------------------------------------------------- 
     1115      REAL(wp) ::   zsig11, zsig12, zsig22 
     1116      REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22 
     1117      REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha 
     1118      REAL(wp) ::   zTany_1, zTany_2 
     1119      REAL(wp) ::   zx, zy, zkxw, kyw, kaw 
     1120      REAL(wp) ::   zinvdx, zinvdy, zinvda 
     1121      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime 
     1122 
     1123      REAL(wp), PARAMETER ::   ppkfriction = 0.45_wp 
    11181124      ! Factor to maintain the same stress as in EVP (see Section 3) 
    11191125      ! Can be set to 1 otherwise 
    1120 !      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
    1121       zinvstressconviso = 1._wp 
    1122   
    1123       zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso  
    1124       !now uses phi as set in higher code 
    1125         
     1126!     REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp/(1._wp+ppkfriction*ppkfriction) 
     1127      REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp 
     1128 
     1129      ! next statement uses pphi set in main module (icedyn_rhg_eap) 
     1130      REAL(wp), PARAMETER ::   ppinvsin = 1._wp/sin(2._wp*pphi) * ppinvstressconviso 
     1131 
    11261132      ! compute eigenvalues, eigenvectors and angles for structure tensor, strain 
    11271133      ! rates 
     
    11321138      zQ12Q12 = rsmall 
    11331139      zQ11Q12 = rsmall 
    1134   
    1135       ! gamma: angle between general coordiantes and principal axis of A  
    1136       ! here Tan2gamma = 2 a12 / (a11 - a22)  
    1137       IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN       
     1140 
     1141      ! gamma: angle between general coordiantes and principal axis of A 
     1142      ! here Tan2gamma = 2 a12 / (a11 - a22) 
     1143      IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN 
    11381144         zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + & 
    11391145                              4._wp*pa12*pa12 ) 
    1140     
    1141          zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2  
     1146 
     1147         zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2 
    11421148         zQ12Q12 = 0.5_wp - ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2 
    1143          zQ11Q12 = pa12*zAngle_denom_gamma                     !CosSin  
    1144       ENDIF 
    1145   
     1149         zQ11Q12 = pa12*zAngle_denom_gamma                     !CosSin 
     1150      ENDIF 
     1151 
    11461152      ! rotation Q*atemp*Q^T 
    1147       zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22  
    1148        
     1153      zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22 
     1154 
    11491155      ! make first principal value the largest 
    11501156      zatempprime = max(zatempprime, 1.0_wp - zatempprime) 
    1151   
     1157 
    11521158      ! 2) strain rate 
    11531159      zdtemp11 = 0.5_wp*(pdivu + ptension) 
     
    11551161      zdtemp22 = 0.5_wp*(pdivu - ptension) 
    11561162 
    1157       ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)  
     1163      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 
    11581164 
    11591165      zQd11Qd11 = 1.0_wp 
     
    11661172                             ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12) 
    11671173 
    1168          zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2  
     1174         zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2 
    11691175         zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2 
    11701176         zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin 
     
    11771183      IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN 
    11781184         zx = atan2(zdtemp2,zdtemp1) 
    1179       ENDIF       
    1180                 
    1181       ! to ensure the angle lies between pi/4 and 9 pi/4  
     1185      ENDIF 
     1186 
     1187      ! to ensure the angle lies between pi/4 and 9 pi/4 
    11821188      IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp 
    1183        
     1189 
    11841190      ! y: angle between major principal axis of strain rate and structure 
    11851191      ! tensor 
    1186       ! y = gamma - alpha  
     1192      ! y = gamma - alpha 
    11871193      ! Expressesed componently with 
    11881194      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 
    1189        
     1195 
    11901196      zTany_1 = zQ11Q12 - zQd11Qd12 
    11911197      zTany_2 = zQ11Q11 - zQd12Qd12 
    1192        
     1198 
    11931199      zy = 0._wp 
    1194        
     1200 
    11951201      IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN 
    11961202         zy = atan2(zTany_1,zTany_2) 
    11971203      ENDIF 
    1198        
     1204 
    11991205      ! to make sure y is between 0 and pi 
    12001206      IF (zy > rpi) zy = zy - rpi 
    12011207      IF (zy < 0)  zy = zy + rpi 
    12021208 
    1203       ! 3) update anisotropic stress tensor  
    1204       zdx   = rpi/real(nx_yield-1,kind=wp) 
    1205       zdy   = rpi/real(ny_yield-1,kind=wp) 
    1206       zda   = 0.5_wp/real(na_yield-1,kind=wp) 
    1207       zinvdx = 1._wp/zdx 
    1208       zinvdy = 1._wp/zdy 
    1209       zinvda = 1._wp/zda 
     1209      ! 3) update anisotropic stress tensor 
     1210      zinvdx = real(nx_yield-1,kind=wp)/rpi 
     1211      zinvdy = real(ny_yield-1,kind=wp)/rpi 
     1212      zinvda = 2._wp*real(na_yield-1,kind=wp) 
    12101213 
    12111214      ! % need 8 coords and 8 weights 
     
    12131216      kx  = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 
    12141217      !!clem kx  = MAX( 1, MIN( nx_yield-1, INT((zx-rpi*0.25_wp-rpi)*zinvdx) + 1  ) ) 
    1215       zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx   
    1216        
     1218      zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 
     1219 
    12171220      ky  = int(zy*zinvdy) + 1 
    12181221      !!clem ky  = MAX( 1, MIN( ny_yield-1, INT(zy*zinvdy) + 1 ) ) 
    1219       kyw = ky - zy*zinvdy  
    1220        
     1222      kyw = ky - zy*zinvdy 
     1223 
    12211224      ka  = int((zatempprime-0.5_wp)*zinvda) + 1 
    12221225      !!clem ka  = MAX( 1, MIN( na_yield-1, INT((zatempprime-0.5_wp)*zinvda) + 1 ) ) 
    12231226      kaw = ka - (zatempprime-0.5_wp)*zinvda 
    1224        
     1227 
    12251228      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 
    12261229!!$      zstemp11r =  zkxw  * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) & 
     
    12481251!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) & 
    12491252!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 
    1250 !!$       
     1253!!$ 
    12511254!!$      zstemp11s =  zkxw  * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) & 
    12521255!!$        & + (1._wp-zkxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) & 
     
    12801283      zstemp12s = s12s(kx,ky,ka) 
    12811284      zstemp22s = s22s(kx,ky,ka) 
    1282        
    1283        
     1285 
     1286 
    12841287      ! Calculate mean ice stress over a collection of floes (Equation 3 in 
    12851288      ! Tsamados 2013) 
    12861289 
    1287       zsig11  = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin 
    1288       zsig12  = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin 
    1289       zsig22  = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin 
     1290      zsig11  = pstrength*(zstemp11r + ppkfriction*zstemp11s) * ppinvsin 
     1291      zsig12  = pstrength*(zstemp12r + ppkfriction*zstemp12s) * ppinvsin 
     1292      zsig22  = pstrength*(zstemp22r + ppkfriction*zstemp22s) * ppinvsin 
    12901293 
    12911294      ! Back - rotation of the stress from principal axes into general coordinates 
     
    13001303      pstressm  = zsgprm11 - zsgprm22 
    13011304 
    1302       ! Compute ridging and sliding functions in general coordinates  
     1305      ! Compute ridging and sliding functions in general coordinates 
    13031306      ! (Equation 11 in Tsamados 2013) 
    13041307      IF (ksub == kndte) THEN 
     
    13071310         zrotstemp12r = zQ11Q11*zstemp12r +       zQ11Q12*(zstemp11r-zstemp22r) & 
    13081311                      - zQ12Q12*zstemp12r 
    1309          zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r &  
     1312         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 
    13101313                      + zQ11Q11*zstemp22r 
    13111314 
     
    13141317         zrotstemp12s = zQ11Q11*zstemp12s +       zQ11Q12*(zstemp11s-zstemp22s) & 
    13151318                      - zQ12Q12*zstemp12s 
    1316          zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s &  
     1319         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 
    13171320                      + zQ11Q11*zstemp22s 
    13181321 
     
    13221325                  + zrotstemp22s*zdtemp22 
    13231326      ENDIF 
    1324    END SUBROUTINE update_stress_rdg  
     1327   END SUBROUTINE update_stress_rdg 
    13251328 
    13261329!======================================================================= 
     
    13311334      !!--------------------------------------------------------------------- 
    13321335      !!                   ***  ROUTINE calc_ffrac  *** 
    1333       !!                      
     1336      !! 
    13341337      !! ** Purpose :   Computes term in evolution equation for structure tensor 
    13351338      !!                which determines the ice floe re-orientation due to fracture 
     
    13461349      REAL (wp) ::   zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 
    13471350 
    1348 !!$      REAL (wp), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation  
    1349       REAL (wp), PARAMETER ::   kfrac = 1.e-3_wp   ! rate of fracture formation  
    1350       REAL (wp), PARAMETER ::   threshold = 0.3_wp  ! critical confinement ratio  
     1351!!$   REAL (wp), PARAMETER ::   ppkfrac = 0.0001_wp   ! rate of fracture formation 
     1352      REAL (wp), PARAMETER ::   ppkfrac = 1.e-3_wp    ! rate of fracture formation 
     1353      REAL (wp), PARAMETER ::   ppthreshold = 0.3_wp  ! critical confinement ratio 
    13511354      !!--------------------------------------------------------------- 
    13521355      ! 
    1353       zsigma11 = 0.5_wp*(pstressp+pstressm)  
    1354       zsigma12 = pstress12  
    1355       zsigma22 = 0.5_wp*(pstressp-pstressm)  
     1356      zsigma11 = 0.5_wp*(pstressp+pstressm) 
     1357      zsigma12 = pstress12 
     1358      zsigma22 = 0.5_wp*(pstressp-pstressm) 
    13561359 
    13571360      ! Here's the change - no longer calculate gamma, 
    13581361      ! use trig formulation, angle signs are kept correct, don't worry 
    1359     
     1362 
    13601363      ! rotate tensor to get into sigma principal axis 
    1361     
    1362       ! here Tan2gamma = 2 sig12 / (sig11 - sig12)  
    1363       ! add rsmall to the denominator to stop 1/0 errors, makes very little  
     1364 
     1365      ! here Tan2gamma = 2 sig12 / (sig11 - sig12) 
     1366      ! add rsmall to the denominator to stop 1/0 errors, makes very little 
    13641367      ! error to the calculated angles < rsmall 
    13651368 
     
    13731376                       zsigma22 ) + 4.0_wp*zsigma12*zsigma12) 
    13741377 
    1375          zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2  
     1378         zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2 
    13761379         zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Sin^2 
    1377          zQ11Q12 = zsigma12*zAngle_denom                      !CosSin  
     1380         zQ11Q12 = zsigma12*zAngle_denom                      !CosSin 
    13781381      ENDIF 
    13791382 
     
    13901393      ! which leads to the loss of their shape, so we again model it through diffusion 
    13911394      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN 
    1392          pmresult11 = - kfrac * (pa11 - zQ12Q12) 
    1393          pmresult12 = - kfrac * (pa12 + zQ11Q12) 
     1395         pmresult11 = - ppkfrac * (pa11 - zQ12Q12) 
     1396         pmresult12 = - ppkfrac * (pa12 + zQ11Q12) 
    13941397 
    13951398      ! Shear faulting 
     
    13971400         pmresult11 = 0.0_wp 
    13981401         pmresult12 = 0.0_wp 
    1399       ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 
    1400          pmresult11 = - kfrac * (pa11 - zQ12Q12) 
    1401          pmresult12 = - kfrac * (pa12 + zQ11Q12) 
     1402      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= ppthreshold)) THEN 
     1403         pmresult11 = - ppkfrac * (pa11 - zQ12Q12) 
     1404         pmresult12 = - ppkfrac * (pa12 + zQ11Q12) 
    14021405 
    14031406      ! Horizontal spalling 
    1404       ELSE  
     1407      ELSE 
    14051408         pmresult11 = 0.0_wp 
    14061409         pmresult12 = 0.0_wp 
     
    14131416      !!--------------------------------------------------------------------- 
    14141417      !!                   ***  ROUTINE rhg_eap_rst  *** 
    1415       !!                      
     1418      !! 
    14161419      !! ** Purpose :   Read or write RHG file in restart file 
    1417       !!       
     1420      !! 
    14181421      !! ** Method  :   use of IOM library 
    14191422      !!---------------------------------------------------------------------- 
     
    14241427      INTEGER  ::   id1, id2, id3, id4, id5   ! local integers 
    14251428      INTEGER  ::   ix, iy, ip, iz, n, ia     ! local integers 
    1426      
     1429 
    14271430      INTEGER, PARAMETER            ::    nz = 100 
    1428        
     1431 
    14291432      REAL(wp) ::   ainit, xinit, yinit, pinit, zinit 
    14301433      REAL(wp) ::   da, dx, dy, dp, dz, a1 
     
    14321435      !!clem 
    14331436      REAL(wp) ::   zw1, zw2, zfac, ztemp 
    1434       REAL(wp) ::   idx, idy, idz 
     1437      REAL(wp) ::   zidx, zidy, zidz 
     1438      REAL(wp) ::   zsaak(6)                  ! temporary array 
    14351439 
    14361440      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp 
     
    15081512!!$                        s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    15091513!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1510 !!$                           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)  
     1514!!$                           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    15111515!!$                        s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    15121516!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
     
    15431547!!$         ENDDO 
    15441548 
    1545          !! faster but still very slow => to be improved          
     1549         !! faster but still very slow => to be improved 
    15461550         zfac = dz/sin(2._wp*pphi) 
    15471551         DO ia = 1, na_yield-1 
     
    15491553            zw2 = w2(ainit+ia*da) 
    15501554            DO iz = 1, nz 
    1551                idz = zinit+iz*dz 
     1555               zidz = zinit+iz*dz 
    15521556               ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz)) 
    15531557               DO iy = 1, ny_yield 
    1554                   idy = yinit+iy*dy 
     1558                  zidy = yinit+iy*dy 
    15551559                  DO ix = 1, nx_yield 
    1556                      idx = xinit+ix*dx 
    1557                      s11r(ix,iy,ia) = s11r(ix,iy,ia) + ztemp * s11kr(idx,idy,idz)*zfac 
    1558                      s12r(ix,iy,ia) = s12r(ix,iy,ia) + ztemp * s12kr(idx,idy,idz)*zfac 
    1559                      s22r(ix,iy,ia) = s22r(ix,iy,ia) + ztemp * s22kr(idx,idy,idz)*zfac  
    1560                      s11s(ix,iy,ia) = s11s(ix,iy,ia) + ztemp * s11ks(idx,idy,idz)*zfac 
    1561                      s12s(ix,iy,ia) = s12s(ix,iy,ia) + ztemp * s12ks(idx,idy,idz)*zfac 
    1562                      s22s(ix,iy,ia) = s22s(ix,iy,ia) + ztemp * s22ks(idx,idy,idz)*zfac 
     1560                     zidx = xinit+ix*dx 
     1561                     call all_skr_sks(zidx,zidy,zidz,zsaak) 
     1562                     zsaak = ztemp*zsaak*zfac 
     1563                     s11r(ix,iy,ia) = s11r(ix,iy,ia) + zsaak(1) 
     1564                     s12r(ix,iy,ia) = s12r(ix,iy,ia) + zsaak(2) 
     1565                     s22r(ix,iy,ia) = s22r(ix,iy,ia) + zsaak(3) 
     1566                     s11s(ix,iy,ia) = s11s(ix,iy,ia) + zsaak(4) 
     1567                     s12s(ix,iy,ia) = s12s(ix,iy,ia) + zsaak(5) 
     1568                     s22s(ix,iy,ia) = s22s(ix,iy,ia) + zsaak(6) 
    15631569                  END DO 
    15641570               END DO 
    15651571            END DO 
    15661572         END DO 
    1567  
    15681573         zfac = 1._wp/sin(2._wp*pphi) 
    15691574         ia = na_yield 
    15701575         DO iy = 1, ny_yield 
    1571             idy = yinit+iy*dy 
     1576            zidy = yinit+iy*dy 
    15721577            DO ix = 1, nx_yield 
    1573                idx = xinit+ix*dx                   
    1574                s11r(ix,iy,ia) = 0.5_wp*s11kr(idx,idy,0._wp)*zfac 
    1575                s12r(ix,iy,ia) = 0.5_wp*s12kr(idx,idy,0._wp)*zfac 
    1576                s22r(ix,iy,ia) = 0.5_wp*s22kr(idx,idy,0._wp)*zfac 
    1577                s11s(ix,iy,ia) = 0.5_wp*s11ks(idx,idy,0._wp)*zfac 
    1578                s12s(ix,iy,ia) = 0.5_wp*s12ks(idx,idy,0._wp)*zfac 
    1579                s22s(ix,iy,ia) = 0.5_wp*s22ks(idx,idy,0._wp)*zfac 
     1578               zidx = xinit+ix*dx 
     1579               call all_skr_sks(zidx,zidy,0._wp,zsaak) 
     1580               zsaak = 0.5_wp*zsaak*zfac 
     1581               s11r(ix,iy,ia) = zsaak(1) 
     1582               s12r(ix,iy,ia) = zsaak(2) 
     1583               s22r(ix,iy,ia) = zsaak(3) 
     1584               s11s(ix,iy,ia) = zsaak(4) 
     1585               s12s(ix,iy,ia) = zsaak(5) 
     1586               s22s(ix,iy,ia) = zsaak(6) 
    15801587            ENDDO 
    15811588         ENDDO 
     
    16111618      REAL(wp) ::   w1 
    16121619      !!------------------------------------------------------------------- 
    1613     
     1620 
    16141621      w1 = -   223.87569446_wp & 
    16151622       &   +  2361.21986630_wp*pa & 
     
    16201627       &   - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 
    16211628       &   +  3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 
    1622     
     1629 
    16231630   END FUNCTION w1 
    16241631 
     
    16311638      REAL(wp) ::   w2 
    16321639      !!------------------------------------------------------------------- 
    1633     
     1640 
    16341641      w2 = -    6670.68911883_wp & 
    16351642       &   +   70222.33061536_wp*pa & 
     
    16401647       &   -  493379.44906738_wp*pa*pa*pa*pa*pa*pa & 
    16411648       &   +  102356.55151800_wp*pa*pa*pa*pa*pa*pa*pa 
    1642     
     1649 
    16431650   END FUNCTION w2 
    16441651 
    1645    FUNCTION s11kr(px,py,pz)  
    1646       !!------------------------------------------------------------------- 
    1647       !! Function : s11kr 
    1648       !!------------------------------------------------------------------- 
     1652   SUBROUTINE all_skr_sks( px, py, pz, allsk ) 
    16491653      REAL(wp), INTENT(in   ) ::   px,py,pz 
    1650     
    1651       REAL(wp) ::   s11kr, zpih 
    1652     
     1654      REAL(wp), INTENT(out  ) ::   allsk(6) 
     1655 
     1656      REAL(wp) ::   zs12r0, zs21r0 
     1657      REAL(wp) ::   zs12s0, zs21s0 
     1658 
     1659      REAL(wp) ::   zpih 
    16531660      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    16541661      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     
    16591666      REAL(wp) ::   zHen1t2, zHen2t1 
    16601667      !!------------------------------------------------------------------- 
    1661     
     1668 
    16621669      zpih = 0.5_wp*rpi 
    1663     
     1670 
    16641671      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    16651672      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     
    16871694      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    16881695      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1689     
     1696 
    16901697      IF (-zIIn1t2>=rsmall) THEN 
    16911698      zHen1t2 = 1._wp 
     
    16931700      zHen1t2 = 0._wp 
    16941701      ENDIF 
    1695     
     1702 
    16961703      IF (-zIIn2t1>=rsmall) THEN 
    16971704      zHen2t1 = 1._wp 
     
    16991706      zHen2t1 = 0._wp 
    17001707      ENDIF 
    1701     
    1702       s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 
    1703  
    1704    END FUNCTION s11kr 
    1705  
    1706    FUNCTION s12kr(px,py,pz) 
     1708 
     1709      !!------------------------------------------------------------------- 
     1710      !! Function : s11kr 
     1711      !!------------------------------------------------------------------- 
     1712      allsk(1) = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 
    17071713      !!------------------------------------------------------------------- 
    17081714      !! Function : s12kr 
    17091715      !!------------------------------------------------------------------- 
    1710       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1711  
    1712       REAL(wp) ::   s12kr, zs12r0, zs21r0, zpih 
    1713  
    1714       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1715       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1716       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1717       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1718       REAL(wp) ::   zd11, zd12, zd22 
    1719       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1720       REAL(wp) ::   zHen1t2, zHen2t1 
    1721       !!------------------------------------------------------------------- 
    1722       zpih = 0.5_wp*rpi 
    1723     
    1724       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1725       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1726       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1727       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1728       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1729       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1730       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1731       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1732       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1733       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1734       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1735       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1736       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1737       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1738       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1739       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1740       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1741       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1742       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1743       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1744       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1745       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1746     
    1747       IF (-zIIn1t2>=rsmall) THEN 
    1748       zHen1t2 = 1._wp 
    1749       ELSE 
    1750       zHen1t2 = 0._wp 
    1751       ENDIF 
    1752     
    1753       IF (-zIIn2t1>=rsmall) THEN 
    1754       zHen2t1 = 1._wp 
    1755       ELSE 
    1756       zHen2t1 = 0._wp 
    1757       ENDIF 
    1758     
    17591716      zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 
    17601717      zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21) 
    1761       s12kr=0.5_wp*(zs12r0+zs21r0) 
    1762     
    1763    END FUNCTION s12kr 
    1764  
    1765    FUNCTION s22kr(px,py,pz) 
     1718      allsk(2)=0.5_wp*(zs12r0+zs21r0) 
    17661719      !!------------------------------------------------------------------- 
    17671720      !! Function : s22kr 
    17681721      !!------------------------------------------------------------------- 
    1769       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1770  
    1771       REAL(wp) ::   s22kr, zpih 
    1772  
    1773       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1774       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1775       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1776       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1777       REAL(wp) ::   zd11, zd12, zd22 
    1778       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1779       REAL(wp) ::   zHen1t2, zHen2t1 
    1780       !!------------------------------------------------------------------- 
    1781       zpih = 0.5_wp*rpi 
    1782  
    1783       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1784       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1785       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1786       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1787       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1788       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1789       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1790       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1791       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1792       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1793       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1794       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1795       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1796       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1797       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1798       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1799       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1800       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1801       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1802       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1803       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1804       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1805  
    1806       IF (-zIIn1t2>=rsmall) THEN 
    1807       zHen1t2 = 1._wp 
    1808       ELSE 
    1809       zHen1t2 = 0._wp 
    1810       ENDIF 
    1811  
    1812       IF (-zIIn2t1>=rsmall) THEN 
    1813       zHen2t1 = 1._wp 
    1814       ELSE 
    1815       zHen2t1 = 0._wp 
    1816       ENDIF 
    1817  
    1818       s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 
    1819  
    1820    END FUNCTION s22kr 
    1821  
    1822    FUNCTION s11ks(px,py,pz) 
     1722      allsk(3) = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 
    18231723      !!------------------------------------------------------------------- 
    18241724      !! Function : s11ks 
    18251725      !!------------------------------------------------------------------- 
    1826       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1827  
    1828       REAL(wp) ::   s11ks, zpih 
    1829  
    1830       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1831       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1832       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1833       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1834       REAL(wp) ::   zd11, zd12, zd22 
    1835       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1836       REAL(wp) ::   zHen1t2, zHen2t1 
    1837       !!------------------------------------------------------------------- 
    1838       zpih = 0.5_wp*rpi 
    1839  
    1840       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1841       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1842       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1843       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1844       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1845       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1846       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1847       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1848       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1849       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1850       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1851       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1852       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1853       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1854       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1855       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1856       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1857       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1858       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1859       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1860       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1861       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1862  
    1863       IF (-zIIn1t2>=rsmall) THEN 
    1864       zHen1t2 = 1._wp 
    1865       ELSE 
    1866       zHen1t2 = 0._wp 
    1867       ENDIF 
    1868  
    1869       IF (-zIIn2t1>=rsmall) THEN 
    1870       zHen2t1 = 1._wp 
    1871       ELSE 
    1872       zHen2t1 = 0._wp 
    1873       ENDIF 
    1874  
    1875       s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 
    1876  
    1877    END FUNCTION s11ks 
    1878  
    1879    FUNCTION s12ks(px,py,pz) 
     1726 
     1727      allsk(4) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 
    18801728      !!------------------------------------------------------------------- 
    18811729      !! Function : s12ks 
    18821730      !!------------------------------------------------------------------- 
    1883       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1884  
    1885       REAL(wp) ::   s12ks, zs12s0, zs21s0, zpih 
    1886  
    1887       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1888       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1889       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1890       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1891       REAL(wp) ::   zd11, zd12, zd22 
    1892       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1893       REAL(wp) ::   zHen1t2, zHen2t1 
    1894       !!------------------------------------------------------------------- 
    1895       zpih = 0.5_wp*rpi 
    1896  
    1897       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1898       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1899       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1900       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1901       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1902       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1903       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1904       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1905       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1906       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1907       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1908       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1909       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1910       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1911       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1912       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1913       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1914       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1915       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1916       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1917       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1918       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1919  
    1920       IF (-zIIn1t2>=rsmall) THEN 
    1921       zHen1t2 = 1._wp 
    1922       ELSE 
    1923       zHen1t2 = 0._wp 
    1924       ENDIF 
    1925  
    1926       IF (-zIIn2t1>=rsmall) THEN 
    1927       zHen2t1 = 1._wp 
    1928       ELSE 
    1929       zHen2t1 = 0._wp 
    1930       ENDIF 
    1931  
    19321731      zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 
    19331732      zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21) 
    1934       s12ks=0.5_wp*(zs12s0+zs21s0) 
    1935  
    1936    END FUNCTION s12ks 
    1937  
    1938    FUNCTION s22ks(px,py,pz)  
     1733      allsk(5)=0.5_wp*(zs12s0+zs21s0) 
    19391734      !!------------------------------------------------------------------- 
    19401735      !! Function : s22ks 
    19411736      !!------------------------------------------------------------------- 
    1942       REAL(wp), INTENT(in   ) ::   px,py,pz 
    1943  
    1944       REAL(wp) ::   s22ks, zpih 
    1945  
    1946       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
    1947       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
    1948       REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
    1949       REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
    1950       REAL(wp) ::   zd11, zd12, zd22 
    1951       REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
    1952       REAL(wp) ::   zHen1t2, zHen2t1 
    1953       !!------------------------------------------------------------------- 
    1954       zpih = 0.5_wp*rpi 
    1955  
    1956       zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
    1957       zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
    1958       zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
    1959       zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
    1960       zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
    1961       zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
    1962       zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
    1963       zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
    1964       zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
    1965       zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
    1966       zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
    1967       zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
    1968       zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
    1969       zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
    1970       zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
    1971       zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    1972       zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    1973       zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    1974       zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    1975       zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
    1976       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
    1977       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    1978  
    1979       IF (-zIIn1t2>=rsmall) THEN 
    1980       zHen1t2 = 1._wp 
    1981       ELSE 
    1982       zHen1t2 = 0._wp 
    1983       ENDIF 
    1984  
    1985       IF (-zIIn2t1>=rsmall) THEN 
    1986       zHen2t1 = 1._wp 
    1987       ELSE 
    1988       zHen2t1 = 0._wp 
    1989       ENDIF 
    1990  
    1991       s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 
    1992  
    1993    END FUNCTION s22ks 
     1737      allsk(6) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 
     1738   END SUBROUTINE all_skr_sks 
    19941739 
    19951740#else 
     
    19971742   !!   Default option         Empty module           NO SI3 sea-ice model 
    19981743   !!---------------------------------------------------------------------- 
     1744   USE par_kind 
     1745   USE lib_mpp 
     1746   CONTAINS 
     1747   SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12, prdg_conv ) 
     1748      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step 
     1749      INTEGER                 , INTENT(in   ) ::   Kmm                                   ! ocean time level index 
     1750      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pstress1_i, pstress2_i, pstress12_i   ! 
     1751      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
     1752      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   paniso_11 , paniso_12                 ! structure tensor components 
     1753      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   prdg_conv                             ! for ridging 
     1754      CALL ctl_stop('EAP rheology is currently dsabled due to issues with AGRIF preprocessing') 
     1755   END SUBROUTINE ice_dyn_rhg_eap 
     1756   SUBROUTINE rhg_eap_rst( cdrw, kt ) 
     1757      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     1758      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
     1759   END SUBROUTINE rhg_eap_rst 
    19991760#endif 
    20001761 
Note: See TracChangeset for help on using the changeset viewer.