Changeset 14117 for NEMO/trunk/tests
- Timestamp:
- 2020-12-07T11:21:42+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
r14021 r14117 6 6 !! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code 7 7 !! 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 9 9 !! 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 11 11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 12 12 !! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + mEVP (Bouillon 2013) … … 14 14 !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] 15 15 !! ! 2019 (S. Rynders, Y. Aksenov, C. Rousset) change into eap rheology from 16 !! CICE code (Tsamados, Heorton) 16 !! CICE code (Tsamados, Heorton) 17 17 !!---------------------------------------------------------------------- 18 18 #if defined key_si3 … … 30 30 USE icevar ! ice_var_sshdyn 31 31 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 34 34 #if defined key_agrif 35 35 USE agrif_ice_interp … … 66 66 INTEGER :: ncvgid ! netcdf file id 67 67 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 69 70 !!---------------------------------------------------------------------- 70 71 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 80 81 !! 81 82 !! ** 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 . 85 86 !! 86 87 !! The points in the C-grid look like this, dear reader … … 90 91 !! | 91 92 !! (ji-1,jj) | (ji,jj) 92 !! --------- 93 !! --------- 93 94 !! | | 94 95 !! | (ji,jj) |------(ji,jj) 95 96 !! | | 96 !! --------- 97 !! --------- 97 98 !! (ji-1,jj-1) (ji,jj-1) 98 99 !! … … 101 102 !! snow total volume (vt_s) per unit area 102 103 !! 103 !! ** Action : - compute u_ice, v_ice : the components of the 104 !! ** Action : - compute u_ice, v_ice : the components of the 104 105 !! sea-ice velocity vector 105 106 !! - compute delta_i, shear_i, divu_i, which are inputs … … 107 108 !! 108 109 !! ** Steps : 0) compute mask at F point 109 !! 1) Compute ice snow mass, ice strength 110 !! 1) Compute ice snow mass, ice strength 110 111 !! 2) Compute wind, oceanic stresses, mass terms and 111 112 !! coriolis terms of the momentum equation … … 166 167 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 167 168 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 169 170 ! 170 171 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear … … 172 173 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 173 174 ! ! 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 175 176 REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses 176 177 REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points … … 192 193 !! --- diags 193 194 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 195 196 !! --- SIMIP diags 196 197 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 199 200 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 200 201 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) 202 203 !!------------------------------------------------------------------- 203 204 204 205 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 205 206 ! 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 ! 208 215 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 211 217 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 212 223 ! 213 224 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... … … 249 260 ! 1) define some variables and initialize arrays 250 261 !------------------------------------------------------------------------------! 251 zrhoco = rho0 * rn_cio 262 zrhoco = rho0 * rn_cio 252 263 !extra code for test case boundary conditions 253 264 zinvw=1._wp/(zrhoco*0.5_wp) … … 270 281 ENDIF 271 282 z1_dtevp = 1._wp / zdtevp 272 273 ! Initialise stress tensor 274 zs1 (:,:) = pstress1_i (:,:) 283 284 ! Initialise stress tensor 285 zs1 (:,:) = pstress1_i (:,:) 275 286 zs2 (:,:) = pstress2_i (:,:) 276 287 zs12(:,:) = pstress12_i(:,:) … … 319 330 ! dt/m at T points (for alpha and beta coefficients) 320 331 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) 321 332 322 333 ! m/dt 323 334 zmU_t(ji,jj) = zmassU * z1_dtevp 324 335 zmV_t(ji,jj) = zmassV * z1_dtevp 325 336 326 337 ! Drag ice-atm. 327 338 ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) … … 353 364 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) 354 365 ! 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 356 367 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 357 368 ! 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 359 370 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 360 371 ! 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 362 373 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 363 374 END_2D … … 377 388 ! ! ==================== ! 378 389 DO jter = 1 , nn_nevp ! loop over jter ! 379 ! ! ==================== ! 390 ! ! ==================== ! 380 391 l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 381 392 ! … … 404 415 & + 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) & 405 416 & ) * 0.25_wp * r1_e1e2t(ji,jj) 406 417 407 418 ! divergence at T points 408 419 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & … … 410 421 & ) * r1_e1e2t(ji,jj) 411 422 zdiv2 = zdiv * zdiv 412 423 413 424 ! tension at T points 414 425 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) & … … 418 429 419 430 ! delta at T points 420 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 431 zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 421 432 422 433 END_2D 423 434 CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp ) 424 435 425 436 ! P/delta at T points 426 437 DO_2D( 1, 1, 1, 1 ) … … 430 441 DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 431 442 432 ! shear at T points 443 ! shear at T points 433 444 zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 434 445 & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 435 446 & ) * 0.25_wp * r1_e1e2t(ji,jj) 436 447 437 448 ! divergence at T points (duplication to avoid communications) 438 449 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 439 450 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 440 451 & ) * r1_e1e2t(ji,jj) 441 452 442 453 ! tension at T points (duplication to avoid communications) 443 454 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) & … … 459 470 zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp) 460 471 zyield12(ji,jj) = zstress12tmp(ji,jj) 461 prdg_conv(ji,jj) = -min( zalphar, 0._wp) 472 prdg_conv(ji,jj) = -min( zalphar, 0._wp) 462 473 ENDIF 463 474 … … 491 502 492 503 DO_2D( 1, 0, 1, 0 ) 493 ! stress12tmp at F points 504 ! stress12tmp at F points 494 505 zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) & 495 506 & + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) & … … 504 515 ! zalph2 = zalph2 - 1._wp 505 516 ENDIF 506 517 507 518 ! stress at F points (zkt/=0 if landfast) 508 519 zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1 … … 570 581 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 571 582 & ) / 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) & 573 584 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 574 585 & ) / ( zbetav + 1._wp ) & … … 630 641 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 631 642 & ) / ( 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 633 644 & ) * zmsk00x(ji,jj) 634 645 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) … … 637 648 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 638 649 & + ( 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 640 651 & ) * zmsk00x(ji,jj) 641 652 ENDIF … … 689 700 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 690 701 & ) / ( 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 692 703 & ) * zmsk00x(ji,jj) 693 704 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) … … 744 755 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 745 756 & + 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 ) & 747 758 & ) * 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 748 759 & ) * zmsk00y(ji,jj) … … 781 792 ! 782 793 !------------------------------------------------------------------------------! 783 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 794 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 784 795 !------------------------------------------------------------------------------! 785 796 DO_2D( 1, 0, 1, 0 ) … … 791 802 792 803 END_2D 793 804 794 805 DO_2D( 0, 0, 0, 0 ) 795 806 796 807 ! tension**2 at T points 797 808 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) & … … 799 810 & ) * r1_e1e2t(ji,jj) 800 811 zdt2 = zdt * zdt 801 812 802 813 zten_i(ji,jj) = zdt 803 814 … … 806 817 & + 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) & 807 818 & ) * 0.25_wp * r1_e1e2t(ji,jj) 808 819 809 820 ! shear at T points 810 821 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) … … 814 825 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 815 826 & ) * r1_e1e2t(ji,jj) 816 827 817 828 ! 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 819 830 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 820 831 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl … … 824 835 & zten_i, 'T', 1.0_wp, zs1 , 'T', 1.0_wp, zs2 , 'T', 1.0_wp, & 825 836 & zs12, 'F', 1.0_wp ) 826 837 827 838 ! --- Store the stress tensor for the next time step --- ! 828 839 pstress1_i (:,:) = zs1 (:,:) … … 841 852 & ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 842 853 ! 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 851 862 ! --- divergence, shear and strength --- ! 852 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence853 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear854 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta855 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength863 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 856 867 857 868 ! --- Stress tensor invariants (SIMIP diags) --- ! … … 859 870 ! 860 871 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 861 ! 872 ! 862 873 DO_2D( 1, 1, 1, 1 ) 863 874 864 875 ! Ice stresses 865 876 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 866 877 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 867 878 ! 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 ) 869 880 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 870 881 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 871 882 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 872 883 873 884 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 874 885 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 875 886 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 876 887 877 888 END_2D 878 889 ! 879 890 ! 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 stress881 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress882 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 883 894 DEALLOCATE ( zsig_I, zsig_II ) 884 895 885 896 ENDIF 886 897 … … 891 902 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 892 903 ! 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 ! 895 906 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 898 909 ! and **deformations** at current iterates 899 910 ! following Lemieux & Dupont (2020) … … 902 913 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 903 914 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 904 915 905 916 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 906 917 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 907 918 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 908 919 909 920 ! Normalized principal stresses (used to display the ellipse) 910 921 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) … … 913 924 END_2D 914 925 ! 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 918 929 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 919 930 920 931 ENDIF 921 932 … … 925 936 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp ) 926 937 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 ) 930 941 ENDIF 931 942 932 943 ! --- 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 938 949 ! --- SIMIP --- ! 939 950 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 944 955 & zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp ) 945 956 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) 952 963 ENDIF 953 964 … … 960 971 DO_2D( 0, 0, 0, 0 ) 961 972 ! 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) 964 975 965 976 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component … … 979 990 980 991 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 982 993 CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s) 983 994 CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport … … 995 1006 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 996 1007 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(:,:) ) 998 1009 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 999 1010 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(:,:) ) 1001 1012 ENDIF 1002 1013 ENDIF 1003 ENDIF 1004 ! 1005 DEALLOCATE( zmsk00, zmsk15 ) 1014 ENDIF 1006 1015 ! 1007 1016 END SUBROUTINE ice_dyn_rhg_eap 1008 1017 1009 1018 1010 1019 SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 1011 1020 !!---------------------------------------------------------------------- 1012 1021 !! *** ROUTINE rhg_cvg *** 1013 !! 1022 !! 1014 1023 !! ** Purpose : check convergence of oce rheology 1015 1024 !! … … 1019 1028 !! This routine is called every sub-iteration, so it is cpu expensive 1020 1029 !! 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) 1022 1031 !!---------------------------------------------------------------------- 1023 1032 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index … … 1026 1035 INTEGER :: it, idtime, istatus 1027 1036 INTEGER :: ji, jj ! dummy loop indices 1028 REAL(wp) :: zresm ! local real 1037 REAL(wp) :: zresm ! local real 1029 1038 CHARACTER(len=20) :: clname 1030 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence1031 1039 !!---------------------------------------------------------------------- 1032 1040 … … 1053 1061 ! time 1054 1062 it = ( kt - 1 ) * kitermax + kiter 1055 1063 1056 1064 ! convergence 1057 1065 IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) … … 1059 1067 ELSE 1060 1068 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 1063 1075 END_2D 1064 1076 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 ) 1070 1078 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1071 1079 ENDIF … … 1077 1085 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) 1078 1086 ENDIF 1079 1087 1080 1088 END SUBROUTINE rhg_cvg 1081 1089 … … 1085 1093 !!--------------------------------------------------------------------- 1086 1094 !! *** SUBROUTINE update_stress_rdg *** 1087 !! 1095 !! 1088 1096 !! ** Purpose : Updates the stress depending on values of strain rate and structure 1089 1097 !! tensor and for the last subcycle step it computes closing and sliding rate … … 1098 1106 INTEGER :: kx ,ky, ka 1099 1107 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 1106 1114 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 1118 1124 ! Factor to maintain the same stress as in EVP (see Section 3) 1119 1125 ! Can be set to 1 otherwise 1120 ! zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction)1121 zinvstressconviso = 1._wp1122 1123 zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso1124 !now uses phi as set in higher code1125 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 1126 1132 ! compute eigenvalues, eigenvectors and angles for structure tensor, strain 1127 1133 ! rates … … 1132 1138 zQ12Q12 = rsmall 1133 1139 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 1138 1144 zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + & 1139 1145 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 1142 1148 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 1146 1152 ! 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 1149 1155 ! make first principal value the largest 1150 1156 zatempprime = max(zatempprime, 1.0_wp - zatempprime) 1151 1157 1152 1158 ! 2) strain rate 1153 1159 zdtemp11 = 0.5_wp*(pdivu + ptension) … … 1155 1161 zdtemp22 = 0.5_wp*(pdivu - ptension) 1156 1162 1157 ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 1163 ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 1158 1164 1159 1165 zQd11Qd11 = 1.0_wp … … 1166 1172 ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12) 1167 1173 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 1169 1175 zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2 1170 1176 zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin … … 1177 1183 IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN 1178 1184 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 1182 1188 IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp 1183 1189 1184 1190 ! y: angle between major principal axis of strain rate and structure 1185 1191 ! tensor 1186 ! y = gamma - alpha 1192 ! y = gamma - alpha 1187 1193 ! Expressesed componently with 1188 1194 ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 1189 1195 1190 1196 zTany_1 = zQ11Q12 - zQd11Qd12 1191 1197 zTany_2 = zQ11Q11 - zQd12Qd12 1192 1198 1193 1199 zy = 0._wp 1194 1200 1195 1201 IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN 1196 1202 zy = atan2(zTany_1,zTany_2) 1197 1203 ENDIF 1198 1204 1199 1205 ! to make sure y is between 0 and pi 1200 1206 IF (zy > rpi) zy = zy - rpi 1201 1207 IF (zy < 0) zy = zy + rpi 1202 1208 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) 1210 1213 1211 1214 ! % need 8 coords and 8 weights … … 1213 1216 kx = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 1214 1217 !!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 1217 1220 ky = int(zy*zinvdy) + 1 1218 1221 !!clem ky = MAX( 1, MIN( ny_yield-1, INT(zy*zinvdy) + 1 ) ) 1219 kyw = ky - zy*zinvdy 1220 1222 kyw = ky - zy*zinvdy 1223 1221 1224 ka = int((zatempprime-0.5_wp)*zinvda) + 1 1222 1225 !!clem ka = MAX( 1, MIN( na_yield-1, INT((zatempprime-0.5_wp)*zinvda) + 1 ) ) 1223 1226 kaw = ka - (zatempprime-0.5_wp)*zinvda 1224 1227 1225 1228 ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 1226 1229 !!$ zstemp11r = zkxw * kyw * kaw * s11r(kx ,ky ,ka ) & … … 1248 1251 !!$ & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) & 1249 1252 !!$ & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 1250 !!$ 1253 !!$ 1251 1254 !!$ zstemp11s = zkxw * kyw * kaw * s11s(kx ,ky ,ka ) & 1252 1255 !!$ & + (1._wp-zkxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & … … 1280 1283 zstemp12s = s12s(kx,ky,ka) 1281 1284 zstemp22s = s22s(kx,ky,ka) 1282 1283 1285 1286 1284 1287 ! Calculate mean ice stress over a collection of floes (Equation 3 in 1285 1288 ! Tsamados 2013) 1286 1289 1287 zsig11 = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin1288 zsig12 = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin1289 zsig22 = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin1290 zsig11 = pstrength*(zstemp11r + ppkfriction*zstemp11s) * ppinvsin 1291 zsig12 = pstrength*(zstemp12r + ppkfriction*zstemp12s) * ppinvsin 1292 zsig22 = pstrength*(zstemp22r + ppkfriction*zstemp22s) * ppinvsin 1290 1293 1291 1294 ! Back - rotation of the stress from principal axes into general coordinates … … 1300 1303 pstressm = zsgprm11 - zsgprm22 1301 1304 1302 ! Compute ridging and sliding functions in general coordinates 1305 ! Compute ridging and sliding functions in general coordinates 1303 1306 ! (Equation 11 in Tsamados 2013) 1304 1307 IF (ksub == kndte) THEN … … 1307 1310 zrotstemp12r = zQ11Q11*zstemp12r + zQ11Q12*(zstemp11r-zstemp22r) & 1308 1311 - zQ12Q12*zstemp12r 1309 zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 1312 zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 1310 1313 + zQ11Q11*zstemp22r 1311 1314 … … 1314 1317 zrotstemp12s = zQ11Q11*zstemp12s + zQ11Q12*(zstemp11s-zstemp22s) & 1315 1318 - zQ12Q12*zstemp12s 1316 zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 1319 zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 1317 1320 + zQ11Q11*zstemp22s 1318 1321 … … 1322 1325 + zrotstemp22s*zdtemp22 1323 1326 ENDIF 1324 END SUBROUTINE update_stress_rdg 1327 END SUBROUTINE update_stress_rdg 1325 1328 1326 1329 !======================================================================= … … 1331 1334 !!--------------------------------------------------------------------- 1332 1335 !! *** ROUTINE calc_ffrac *** 1333 !! 1336 !! 1334 1337 !! ** Purpose : Computes term in evolution equation for structure tensor 1335 1338 !! which determines the ice floe re-orientation due to fracture … … 1346 1349 REAL (wp) :: zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 1347 1350 1348 !!$ REAL (wp), PARAMETER :: kfrac = 0.0001_wp ! rate of fracture formation1349 REAL (wp), PARAMETER :: kfrac = 1.e-3_wp ! rate of fracture formation1350 REAL (wp), PARAMETER :: threshold = 0.3_wp ! critical confinement ratio1351 !!$ 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 1351 1354 !!--------------------------------------------------------------- 1352 1355 ! 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) 1356 1359 1357 1360 ! Here's the change - no longer calculate gamma, 1358 1361 ! use trig formulation, angle signs are kept correct, don't worry 1359 1362 1360 1363 ! 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 1364 1367 ! error to the calculated angles < rsmall 1365 1368 … … 1373 1376 zsigma22 ) + 4.0_wp*zsigma12*zsigma12) 1374 1377 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 1376 1379 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 1378 1381 ENDIF 1379 1382 … … 1390 1393 ! which leads to the loss of their shape, so we again model it through diffusion 1391 1394 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) 1394 1397 1395 1398 ! Shear faulting … … 1397 1400 pmresult11 = 0.0_wp 1398 1401 pmresult12 = 0.0_wp 1399 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN1400 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) 1402 1405 1403 1406 ! Horizontal spalling 1404 ELSE 1407 ELSE 1405 1408 pmresult11 = 0.0_wp 1406 1409 pmresult12 = 0.0_wp … … 1413 1416 !!--------------------------------------------------------------------- 1414 1417 !! *** ROUTINE rhg_eap_rst *** 1415 !! 1418 !! 1416 1419 !! ** Purpose : Read or write RHG file in restart file 1417 !! 1420 !! 1418 1421 !! ** Method : use of IOM library 1419 1422 !!---------------------------------------------------------------------- … … 1424 1427 INTEGER :: id1, id2, id3, id4, id5 ! local integers 1425 1428 INTEGER :: ix, iy, ip, iz, n, ia ! local integers 1426 1429 1427 1430 INTEGER, PARAMETER :: nz = 100 1428 1431 1429 1432 REAL(wp) :: ainit, xinit, yinit, pinit, zinit 1430 1433 REAL(wp) :: da, dx, dy, dp, dz, a1 … … 1432 1435 !!clem 1433 1436 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 1435 1439 1436 1440 REAL(wp), PARAMETER :: eps6 = 1.0e-6_wp … … 1508 1512 !!$ s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1509 1513 !!$ 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) 1511 1515 !!$ s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1512 1516 !!$ exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & … … 1543 1547 !!$ ENDDO 1544 1548 1545 !! faster but still very slow => to be improved 1549 !! faster but still very slow => to be improved 1546 1550 zfac = dz/sin(2._wp*pphi) 1547 1551 DO ia = 1, na_yield-1 … … 1549 1553 zw2 = w2(ainit+ia*da) 1550 1554 DO iz = 1, nz 1551 idz = zinit+iz*dz1555 zidz = zinit+iz*dz 1552 1556 ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz)) 1553 1557 DO iy = 1, ny_yield 1554 idy = yinit+iy*dy1558 zidy = yinit+iy*dy 1555 1559 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) 1563 1569 END DO 1564 1570 END DO 1565 1571 END DO 1566 1572 END DO 1567 1568 1573 zfac = 1._wp/sin(2._wp*pphi) 1569 1574 ia = na_yield 1570 1575 DO iy = 1, ny_yield 1571 idy = yinit+iy*dy1576 zidy = yinit+iy*dy 1572 1577 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) 1580 1587 ENDDO 1581 1588 ENDDO … … 1611 1618 REAL(wp) :: w1 1612 1619 !!------------------------------------------------------------------- 1613 1620 1614 1621 w1 = - 223.87569446_wp & 1615 1622 & + 2361.21986630_wp*pa & … … 1620 1627 & - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 1621 1628 & + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 1622 1629 1623 1630 END FUNCTION w1 1624 1631 … … 1631 1638 REAL(wp) :: w2 1632 1639 !!------------------------------------------------------------------- 1633 1640 1634 1641 w2 = - 6670.68911883_wp & 1635 1642 & + 70222.33061536_wp*pa & … … 1640 1647 & - 493379.44906738_wp*pa*pa*pa*pa*pa*pa & 1641 1648 & + 102356.55151800_wp*pa*pa*pa*pa*pa*pa*pa 1642 1649 1643 1650 END FUNCTION w2 1644 1651 1645 FUNCTION s11kr(px,py,pz) 1646 !!------------------------------------------------------------------- 1647 !! Function : s11kr 1648 !!------------------------------------------------------------------- 1652 SUBROUTINE all_skr_sks( px, py, pz, allsk ) 1649 1653 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 1653 1660 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1654 1661 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 … … 1659 1666 REAL(wp) :: zHen1t2, zHen2t1 1660 1667 !!------------------------------------------------------------------- 1661 1668 1662 1669 zpih = 0.5_wp*rpi 1663 1670 1664 1671 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1665 1672 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) … … 1687 1694 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1688 1695 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1689 1696 1690 1697 IF (-zIIn1t2>=rsmall) THEN 1691 1698 zHen1t2 = 1._wp … … 1693 1700 zHen1t2 = 0._wp 1694 1701 ENDIF 1695 1702 1696 1703 IF (-zIIn2t1>=rsmall) THEN 1697 1704 zHen2t1 = 1._wp … … 1699 1706 zHen2t1 = 0._wp 1700 1707 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) 1707 1713 !!------------------------------------------------------------------- 1708 1714 !! Function : s12kr 1709 1715 !!------------------------------------------------------------------- 1710 REAL(wp), INTENT(in ) :: px,py,pz1711 1712 REAL(wp) :: s12kr, zs12r0, zs21r0, zpih1713 1714 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i221715 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i221716 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i221717 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i221718 REAL(wp) :: zd11, zd12, zd221719 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t21720 REAL(wp) :: zHen1t2, zHen2t11721 !!-------------------------------------------------------------------1722 zpih = 0.5_wp*rpi1723 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 * zd221744 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd221745 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd221746 1747 IF (-zIIn1t2>=rsmall) THEN1748 zHen1t2 = 1._wp1749 ELSE1750 zHen1t2 = 0._wp1751 ENDIF1752 1753 IF (-zIIn2t1>=rsmall) THEN1754 zHen2t1 = 1._wp1755 ELSE1756 zHen2t1 = 0._wp1757 ENDIF1758 1759 1716 zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 1760 1717 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) 1766 1719 !!------------------------------------------------------------------- 1767 1720 !! Function : s22kr 1768 1721 !!------------------------------------------------------------------- 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) 1823 1723 !!------------------------------------------------------------------- 1824 1724 !! Function : s11ks 1825 1725 !!------------------------------------------------------------------- 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) 1880 1728 !!------------------------------------------------------------------- 1881 1729 !! Function : s12ks 1882 1730 !!------------------------------------------------------------------- 1883 REAL(wp), INTENT(in ) :: px,py,pz1884 1885 REAL(wp) :: s12ks, zs12s0, zs21s0, zpih1886 1887 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i221888 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i221889 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i221890 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i221891 REAL(wp) :: zd11, zd12, zd221892 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t21893 REAL(wp) :: zHen1t2, zHen2t11894 !!-------------------------------------------------------------------1895 zpih = 0.5_wp*rpi1896 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 * zd221917 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd221918 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd221919 1920 IF (-zIIn1t2>=rsmall) THEN1921 zHen1t2 = 1._wp1922 ELSE1923 zHen1t2 = 0._wp1924 ENDIF1925 1926 IF (-zIIn2t1>=rsmall) THEN1927 zHen2t1 = 1._wp1928 ELSE1929 zHen2t1 = 0._wp1930 ENDIF1931 1932 1731 zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 1933 1732 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) 1939 1734 !!------------------------------------------------------------------- 1940 1735 !! Function : s22ks 1941 1736 !!------------------------------------------------------------------- 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 1994 1739 1995 1740 #else … … 1997 1742 !! Default option Empty module NO SI3 sea-ice model 1998 1743 !!---------------------------------------------------------------------- 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 1999 1760 #endif 2000 1761
Note: See TracChangeset
for help on using the changeset viewer.
