Changeset 13544
- Timestamp:
- 2020-09-30T09:03:49+02:00 (3 years ago)
- Location:
- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE
- Files:
-
- 3 edited
-
icedyn_rhg.F90 (modified) (3 diffs)
-
icedyn_rhg_evp.F90 (modified) (1 diff)
-
icedyn_rhg_vp.F90 (modified) (33 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90
r13536 r13544 17 17 USE ice ! sea-ice: variables 18 18 USE icedyn_rhg_evp ! sea-ice: EVP rheology 19 USE icedyn_rhg_vp ! sea-ice: VP rheology 19 20 USE icectl ! sea-ice: control prints 20 21 ! … … 82 83 ! !---------------------------------! 83 84 CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 84 !85 END SELECT86 87 85 ! !---------------------------------------------! 88 86 CASE( np_rhgVP ) ! Viscous-Plastic rheology ! … … 166 164 ! 167 165 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'READ' ) !* read or initialize all required files 168 IF( ln_rhg_VP ) CALL rhg_lsr_rst( 'READ' ) !* read or initialize all required files169 166 ! 170 167 END SUBROUTINE ice_dyn_rhg_init -
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90
r13536 r13544 169 169 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 170 170 !! --- diags 171 REAL(wp) :: zsig1, zsig2, zsig12 171 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 172 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i 173 173 !! --- SIMIP diags -
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r13536 r13544 1 ! TODOLOIST2 !3 ! Commit4 ! Compile5 !6 ! - Re-calculate internal force diagnostic (currently incorrect, see end of the routine)7 ! - QC:compare code to mitGCM8 ! - QC: comparing EVP and VP codes9 !10 ! Clarify11 ! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ?12 ! --- Swap of independent term in the U-V solvers ? -> JF Lemieux says maybe13 ! --- Is stress tensor for restart needed ? No!!!14 !15 ! Test16 ! --- Try ADI for u-v solver17 !18 ! Add19 ! - Charge ellipse as an output (good one from Lemieux and Dupont)20 ! - Bottom drag (landfast ice)21 ! - Tensile strength22 ! - agrif23 ! - bdy24 !25 ! Write documentation26 !27 ! Optimize28 ! - subroutinize common parts (diagnostics)29 ! - namelist: rationalize common parameters EVP/VP30 31 32 1 MODULE icedyn_rhg_vp 33 2 !!====================================================================== … … 73 42 74 43 !! for convergence tests 75 INTEGER :: ncvgid ! netcdf file id44 INTEGER :: ncvgid ! netcdf file id 76 45 INTEGER :: nvarid_ucvg ! netcdf variable id 77 46 INTEGER :: nvarid_ures … … 82 51 INTEGER :: nvarid_mke 83 52 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 53 84 54 !!---------------------------------------------------------------------- 85 55 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 157 127 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 158 128 !! 159 LOGICAL :: ll_u_iterate, ll_v iterate ! continue iteration or not129 LOGICAL :: ll_u_iterate, ll_v_iterate ! continue iteration or not 160 130 ! 161 131 INTEGER :: ji, jj, jn ! dummy loop indices … … 171 141 REAL(wp) :: zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 172 142 REAL(wp) :: zdeltat, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 173 REAL(wp) :: zp_deltastar_f 143 REAL(wp) :: zp_deltastar_f ! 144 REAL(wp) :: zu_cV, zv_cU ! 174 145 REAL(wp) :: zfac, zfac1, zfac2, zfac3 175 146 REAL(wp) :: zt12U, zt11U, zt22U, zt21U, zt122U, zt121U … … 178 149 ! 179 150 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice 180 REAL(wp), DIMENSION(jpi,jpj) :: za U , zaV ! ice fraction on U/V points151 REAL(wp), DIMENSION(jpi,jpj) :: za_iU , za_iV ! ice fraction on U/V points 181 152 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! Acceleration term contribution to RHS 182 153 REAL(wp), DIMENSION(jpi,jpj) :: zmassU_t, zmassV_t ! Mass per unit area divided by time step … … 191 162 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 192 163 REAL(wp), DIMENSION(jpi,jpj) :: zu_c, zv_c ! "current" ice velocity (m/s), average of previous two OL iterates 193 REAL(wp), DIMENSION(jpi,jpj) :: zu_cV, zv_cU ! "current" u (v) ice velocity at V (U) point (m/s)194 164 REAL(wp), DIMENSION(jpi,jpj) :: zu_b, zv_b ! velocity at previous sub-iterate 195 165 REAL(wp), DIMENSION(jpi,jpj) :: zuerr, zverr ! absolute U/Vvelocity difference between current and previous sub-iterates … … 209 179 REAL(wp), DIMENSION(jpi,jpj) :: zrhsu, zrhsv ! U and V RHS 210 180 REAL(wp), DIMENSION(jpi,jpj) :: zAU, zBU, zCU, zDU, zEU ! Linear system coefficients, U equation 211 REAL(wp), DIMENSION(jpi,jpj) :: zAV, zBV, zCV, zDV, zEV , zFV! Linear system coefficients, V equation181 REAL(wp), DIMENSION(jpi,jpj) :: zAV, zBV, zCV, zDV, zEV ! Linear system coefficients, V equation 212 182 REAL(wp), DIMENSION(jpi,jpj) :: zFU, zFU_prime, zBU_prime ! Rearranged linear system coefficients, U equation 213 183 REAL(wp), DIMENSION(jpi,jpj) :: zFV, zFV_prime, zBV_prime ! Rearranged linear system coefficients, V equation … … 222 192 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 223 193 !! --- diags 224 REAL(wp) :: zsig1, zsig2, zsig12 194 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 225 195 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12, zs12f, zten_i ! stress tensor components 226 196 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 227 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s, SIMIP) 228 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s, SIMIP) 229 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s, SIMIP) 230 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s, SIMIP) 231 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s, SIMIP) 232 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s, SIMIP) 197 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztaux_oi, ztauy_oi 198 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice, zdiag_ymtrp_ice ! X/Y-component of ice mass transport (kg/s, SIMIP) 199 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw, zdiag_ymtrp_snw ! X/Y-component of snow mass transport (kg/s, SIMIP) 200 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp, zdiag_yatrp ! X/Y-component of area transport (m2/s, SIMIP) 233 201 234 202 !!---------------------------------------------------------------------------------------------------------------------- … … 251 219 END DO 252 220 253 IF ( ln_zebra_vp ) THEN; nn_ nzebra_vp = 2; ELSE; nn_nzebra_vp = 1; ENDIF254 255 nn_nvp = nn_ out_vp * nn_inn_vp ! maximum number of iterations221 IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2; ELSE; nn_zebra_vp = 1; ENDIF 222 223 nn_nvp = nn_nout_vp * nn_ninn_vp ! maximum number of iterations 256 224 257 225 zrhoco = rau0 * rn_cio … … 341 309 342 310 ! Ice fraction at U-V points 343 za U(ji,jj)= 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)344 za V(ji,jj)= 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)311 za_iU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 312 za_iV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 345 313 346 314 ! Snow and ice mass at U-V points … … 368 336 369 337 ! Wind stress 370 ztaux_ai(ji,jj) = za U(ji,jj) * utau_ice(ji,jj)371 ztauy_ai(ji,jj) = za V(ji,jj) * vtau_ice(ji,jj)338 ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj) 339 ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj) 372 340 373 341 ! Force due to sea surface tilt(- m*g*GRAD(ssh)) … … 380 348 381 349 ! Mask for lots of ice (1) or little ice (0) 382 IF( zmassU <= zmmin .AND. za U(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp350 IF( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 383 351 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 384 IF( zmassV <= zmmin .AND. za V(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp352 IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 385 353 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 386 354 … … 499 467 500 468 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 501 zCwU(ji,jj) = za U(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) ) &502 & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )503 zCwV(ji,jj) = za V(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) ) &504 & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )469 zCwU(ji,jj) = za_iU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) ) & 470 & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 471 zCwV(ji,jj) = za_iV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) ) & 472 & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 505 473 506 474 !--- Ocean-ice drag contributions to RHS … … 577 545 ( e2u(ji,jj) * ( zs1_rhsu(ji+1,jj) - zs1_rhsu(ji,jj) ) & 578 546 & + r1_e2u(ji,jj) * ( e2t(ji+1,jj) * e2t(ji+1,jj) * zs2_rhsu(ji+1,jj) - e2t(ji,jj) * e2t(ji,jj) * zs2_rhsu(ji,jj) ) & 579 & + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj) * e1f(ji,jj) * zs12_rhsu(ji,jj) - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) ) 547 & + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj) * e1f(ji,jj) * zs12_rhsu(ji,jj) - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) ) ) 580 548 581 549 ! --- V component of internal force contribution to RHS at V points … … 583 551 & ( e1v(ji,jj) * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) ) & 584 552 & + r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj) * zs2_rhsv(ji,jj) ) & 585 & + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj) * e2f(ji,jj) * zs12_rhsv(ji,jj) - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) 553 & + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj) * e2f(ji,jj) * zs12_rhsv(ji,jj) - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) ) 586 554 587 555 END DO … … 723 691 ! what follows could be subroutinized... 724 692 725 DO jn = 1, nn_ nzebra_vp ! "zebra" loop (! red-black-sor!!! )693 DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) 726 694 727 695 ! OPT: could be even better optimized with a true red-black SOR … … 735 703 zBU_prime(:,:) = 0._wp 736 704 737 DO jj = jj_min, jpj - 1, nn_ nzebra_vp705 DO jj = jj_min, jpj - 1, nn_zebra_vp 738 706 739 707 !------------------------------------------- … … 756 724 757 725 ! right hand side 758 zFU(ji,jj) = ( zrhsu(ji,jj) & ! right-hand side terms759 & + zAA3 ! boundary condition translation760 & + zDU(ji,jj) * u_ice(ji,jj-1) ! internal force, j-1761 & + zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1726 zFU(ji,jj) = ( zrhsu(ji,jj) & ! right-hand side terms 727 & + zAA3 & ! boundary condition translation 728 & + zDU(ji,jj) * u_ice(ji,jj-1) & ! internal force, j-1 729 & + zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1 762 730 763 731 END DO … … 806 774 ! MV OPT: what follows could be subroutinized... 807 775 808 DO jn = 1, nn_ nzebra_vp ! "zebra" loop776 DO jn = 1, nn_zebra_vp ! "zebra" loop 809 777 810 778 IF ( jn == 1 ) THEN ; ji_min = 2 … … 816 784 zBV_prime(:,:) = 0._wp 817 785 818 DO ji = ji_min, jpi - 1, nn_ nzebra_vp786 DO ji = ji_min, jpi - 1, nn_zebra_vp 819 787 820 788 !!! It is intentional to have a ji then jj loop for V-velocity … … 836 804 837 805 ! right hand side 838 zFV(ji,jj) = ( zrhsv(ji,jj) & ! right-hand side terms839 & + zAA3 ! boundary condition translation840 & + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1841 & + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1806 zFV(ji,jj) = ( zrhsv(ji,jj) & ! right-hand side terms 807 & + zAA3 & ! boundary condition translation 808 & + zDV(ji,jj) * v_ice(ji-1,jj) & ! internal force, j-1 809 & + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1 842 810 843 811 END DO … … 855 823 v_ice(ji,jpj-1) = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1) ! last row 856 824 DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed 857 v_ice(ji,jj) = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj ) / zBV_prime(ji,jj)825 v_ice(ji,jj) = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) / zBV_prime(ji,jj) 858 826 END DO 859 827 … … 902 870 903 871 ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 904 IF ( i_inn > 1 &&zuerr_max > rn_uerr_max_vp ) THEN905 IF ( lwp ) " VP rheology error was too large : ", zu_err_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped "872 IF ( i_inn > 1 .AND. zuerr_max > rn_uerr_max_vp ) THEN 873 IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 906 874 ll_u_iterate = .FALSE. 907 875 ENDIF … … 909 877 ! - Stop if error small enough 910 878 IF ( zuerr_max < rn_uerr_min_vp ) THEN 911 IF ( lwp ) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! "879 IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! " 912 880 ll_u_iterate = .FALSE. 913 881 ENDIF … … 933 901 934 902 ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 935 IF ( i_inn > 1 &&zverr_max > rn_uerr_max_vp ) THEN936 IF ( lwp ) " VP rheology error was too large : ", zv_err_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped "903 IF ( i_inn > 1 .AND. zverr_max > rn_uerr_max_vp ) THEN 904 IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zverr_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 937 905 ll_v_iterate = .FALSE. 938 906 ENDIF 939 907 940 908 ! - Stop if error small enough 941 IF ( zverr_max < rn_ verr_min) THEN942 IF ( lwp ) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! "909 IF ( zverr_max < rn_uerr_min_vp ) THEN 910 IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! " 943 911 ll_v_iterate = .FALSE. 944 912 ENDIF … … 965 933 !------------------------------------------------------------------------------! 966 934 967 END ! End of outer loop (i_out) =============================================================================================935 END DO ! End of outer loop (i_out) ============================================================================================= 968 936 969 937 !------------------------------------------------------------------------------! … … 1077 1045 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 1078 1046 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 1079 ! 1080 CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 1081 & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 1047 1048 ALLOCATE( ztaux_oi(jpi,jpj) , ztauy_oi(jpi,jpj) ) 1049 1050 !--- Recalculate oceanic stress at last inner iteration 1051 DO jj = 2, jpj - 1 1052 DO ji = 2, jpi - 1 1053 1054 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 1055 zu_cV = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 1056 zv_cU = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 1057 1058 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 1059 zCwU(ji,jj) = za_iU(ji,jj) * zrhoco * SQRT( ( u_ice(ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) ) & 1060 & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 1061 zCwV(ji,jj) = za_iV(ji,jj) * zrhoco * SQRT( ( v_ice(ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) ) & 1062 & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 1063 1064 !--- Ocean-ice stress 1065 ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) ) 1066 ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) ) 1067 1068 END DO 1069 END DO 1070 1071 ! 1072 CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, & 1073 ! & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 1082 1074 ! 1083 1075 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) … … 1085 1077 CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 1086 1078 CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 1087 CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 1088 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 1079 ! CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 1080 ! CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 1081 1082 DEALLOCATE( ztaux_oi , ztauy_oi ) 1089 1083 1090 1084 ENDIF … … 1167 1161 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 1168 1162 & iom_use('corstrx') .OR. iom_use('corstry') ) THEN 1163 1164 ! --- Recalculate Coriolis stress at last inner iteration 1165 DO jj = 2, jpj - 1 1166 DO ji = 2, jpi - 1 1167 ! --- U-component 1168 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 1169 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 1170 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 1171 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 1172 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 1173 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 1174 END DO 1175 END DO 1169 1176 ! 1170 1177 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & 1171 1178 & zCorU, 'U', -1., zCorV, 'V', -1. ) 1172 1179 ! 1173 1180 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) 1174 1181 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) 1175 1182 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) 1176 1183 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) 1184 1177 1185 ENDIF 1178 1186 1179 1187 IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1180 1188 1181 1182 ! Recalculate internal forces (divergence of stress tensor) 1189 ! Recalculate internal forces (divergence of stress tensor) at last inner iteration 1183 1190 DO jj = 2, jpj - 1 1184 1191 DO ji = 2, jpi - 1 … … 1291 1298 INTEGER :: it, idtime, istatus 1292 1299 INTEGER :: ji, jj ! dummy loop indices 1293 REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, z mke, zu, zv ! local scalars1294 REAL(wp), DIMENSION(jpi,jpj) :: zu_res(:,:), zv_res(:,:), zvel2(:,:)! local arrays1300 REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 1301 REAL(wp), DIMENSION(jpi,jpj) :: zu_res, zv_res, zvel2 ! local arrays 1295 1302 1296 1303 CHARACTER(len=20) :: clname 1297 :: zres ! check convergence1298 1304 !!---------------------------------------------------------------------- 1299 1305 … … 1352 1358 DO jj = 2, jpj - 1 1353 1359 DO ji = 2, jpi - 1 1354 zu_res(ji,jj) = zrhsu(ji,jj) + zDU(ji,jj) * pu(ji,jj-1) + zEU(ji,jj) * pu(ji,jj+1) &1355 & - zAU(ji,jj) * pu(ji-1,jj) - zBU(ji,jj) * pu(ji,jj) - zCU(ji,jj) * pu(ji+1,jj)1360 zu_res(ji,jj) = prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1361 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) 1356 1362 1357 zv_res(ji,jj) = zrhsv(ji,jj) + zDV(ji,jj) * pv(ji-1,jj) + zEV(ji,jj) * pv(ji+1,jj) &1358 & - zAV(ji,jj) * pv(ji,jj-1) - zBV(ji,jj) * pv(ji,jj) - zCV(ji,jj) * pv(ji,jj+1)1363 zv_res(ji,jj) = prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1364 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) 1359 1365 END DO 1360 1366 END DO 1361 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,: ) )1362 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,: ) )1363 zu_res_mean = SQRT( zu_res mean / pglob_area )1364 zv_res_mean = SQRT( zv_res mean / pglob_area )1365 zvelres = MEAN( zu_res_mean,zv_res_mean )1367 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:,1) ) 1368 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:,1) ) 1369 zu_res_mean = SQRT( zu_res_mean / pglob_area ) 1370 zv_res_mean = SQRT( zv_res_mean / pglob_area ) 1371 zvelres = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1366 1372 1367 1373 ! -- Global mean kinetic energy per unit area (J/m2) … … 1374 1380 END DO 1375 1381 1376 zmke = 0.5_wp * glob sum( 'ice_rhg_vp', zmt(:,:) * e1e2t(:,:) * zvel2(:,:) )1382 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) 1377 1383 1378 1384 ! ! ==================== !
Note: See TracChangeset
for help on using the changeset viewer.
