- Timestamp:
- 2020-09-25T11:01:23+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r13493 r13522 1 1 ! TODOLIST 2 3 !4 ! Code residual norm of the linear system5 ! Verify convergence check6 2 ! 7 3 ! Define all symbols 8 ! Check uses 4 ! - Do viscosity smoothing with sum (differentiable rheology) 5 ! - Re-calculate internal force diagnostic (end of the routine) 6 ! - Do proper masking of output fileds with ice mass (zmsk00 see EVP routine) 7 ! 8 ! quality control - compare code to mitGCM 9 ! quality control - comparing EVP and VP codes 10 ! 11 ! 12 ! Commit 9 13 ! Compile 10 14 ! 11 ! C heck12 ! --- lbclnks --> what is the policy?13 ! --- boundary conditions --> how to enforce them - is the fmask strategy sufficient?14 ! --- change of independent term in the U-V solvers?15 ! --- loop ji then jj in the V-solver, is it well-written15 ! Clarify 16 ! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ? 17 ! --- Swap of independent term in the U-V solvers ? 18 ! --- Is stress tensor for restart needed ? 19 ! --- Is stress tensor calculated at the end of the time step 16 20 ! 17 ! Subroutinize tridiagonal systems ? 21 ! Test 22 ! --- Can we add the 15% mask in the convergence criteria ? 23 ! --- Try ADI for u-v solver 24 ! 25 ! Add 26 ! - Charge ellipse as an output (good one from Lemieux and Dupont) 27 ! - Bottom drag 28 ! - Tensile strength 29 ! - agrif 30 ! - bdy 31 ! 18 32 ! Write documentation 19 33 ! 20 21 ! Convergence quality criteria 22 ! - average KE (for outer loop iterations) (see Lemieux et al, 2009) 23 ! - residual of the linear system (for inner loop iterations) 24 ! - velocity difference ? 25 26 27 ! add bottom drag 28 ! add tensile strength 29 30 ! quality control - compare to mitGCM 31 ! quality control - comparing EVP and VP 34 ! Optimize 35 ! - subroutinize common parts (diagnostics) 36 ! - namelist: rationalize common parameters EVP/VP 37 32 38 33 39 MODULE icedyn_rhg_vp … … 75 81 PUBLIC rhg_vp_rst ! called by icedyn_rhg.F90 76 82 77 !! * Substitutions78 # include "vectopt_loop_substitute.h90"79 80 83 !! for convergence tests 81 84 INTEGER :: ncvgid ! netcdf file id 82 INTEGER :: nvarid ! netcdf variable id 85 INTEGER :: nvarid_ucvg ! netcdf variable id 86 INTEGER :: nvarid_ures 87 INTEGER :: nvarid_vres 88 INTEGER :: nvarid_velres 89 INTEGER :: nvarid_udif 90 INTEGER :: nvarid_vdif 91 INTEGER :: nvarid_mke 83 92 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 84 93 !!---------------------------------------------------------------------- … … 87 96 !! Software governed by the CeCILL license (see ./LICENSE) 88 97 !!---------------------------------------------------------------------- 98 89 99 CONTAINS 90 100 … … 95 105 !! VP-LSR-C-grid 96 106 !! 97 !! ** purpose : determines sea ice drift from wind stress, ice-ocean107 !! ** Purpose : determines sea ice drift from wind stress, ice-ocean 98 108 !! stress and sea-surface slope. Internal forces assume viscous-plastic rheology (Hibler, 1979) 109 !! 110 !! ** Method 99 111 !! 100 112 !! The resolution algorithm follows from Zhang and Hibler (1997) and Losch (2010) 113 !! with elements from Lemieux and Tremblay (2008) and Lemieux and Tremblay (2009) 101 114 !! 102 !! It is based on a linear slope relaxation algorithm (LSR) 103 !! 104 !! * Time stepping 105 !! 106 !! The idea is to arrange the momentum equations as follows: 115 !! The components of the momentum equations are arranged following the ideas of Zhang and Hibler (1997) 116 !! 107 117 !! f1(u) = g1(v) 108 118 !! f2(v) = g2(v) 109 119 !! 110 !! and use two sub-time-step levels (outer and inner) to solve them. 111 !! 112 !! The right-hand side (RHS, g1 & g2) are expressed explicitly, and the left-hand side (LHS, f1, f2) implicitly 113 !! 114 !! The inner loop makes ~1500 steps max to solve a linear system problem, targetting to solve the implicit terms 115 !! The outer loop makes ~10 steps and updates the explicit terms 116 !! Outer loop uses a modified euler time step, with uC = 0.5 * ( uC + u ) 120 !! The right-hand side (RHS) is explicit 121 !! The left-hand side (LHS) is implicit 122 !! Coriolis is part of explicit terms, whereas ice-ocean drag is implicit 123 !! 124 !! Two iteration levels (outer and inner loops) are used to solve the equations 125 !! 126 !! The outer loop (OL, typically 10 iterations) is there to deal with the (strong) non-linearities in the equation 127 !! 128 !! The inner loop (IL, typically 1500 iterations) is there to solve the linear problem with a line-successive-relaxation algorithm 129 !! 130 !! The velocity used in the non-linear terms uses a "modified euler time step" (not sure its the correct term), 131 !!! with uk = ( uk-1 + uk-2 ) / 2. 117 132 !! 118 133 !! * Spatial discretization … … 143 158 !! ** Notes : 144 159 !! 145 !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010. 160 !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010., Lemieux et al., 2008, 2009, ... 146 161 !! 147 162 !! 148 163 !!------------------------------------------------------------------- 149 164 !! 150 !!! INTEGER , INTENT(in ) :: kt ! time step 151 !!! REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i ! 152 !!! REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 153 !!! !! 154 !!! INTEGER :: ji, jj ! dummy loop indices 155 !!! INTEGER :: jter ! local integers 156 !!! ! 157 !!! REAL(wp) :: zrhoco ! rau0 * rn_cio 158 !!! REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling 159 !!! REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 160 !!! REAL(wp) :: zalph1, z1_alph1, zalph2, z1_alph2 ! alpha coef from Bouillon 2009 or Kimmritz 2017 161 !!! REAl(wp) :: zbetau, zbetav 162 !!! REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 163 !!! REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 164 !!! REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 165 !!! REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 166 !!! REAL(wp) :: zvCr ! critical ice volume above which ice is landfast 167 !!! ! 168 !!! REAL(wp) :: zintb, zintn ! dummy argument 169 !!! REAL(wp) :: zfac_x, zfac_y 170 !!! REAL(wp) :: zshear, zdum1, zdum2 171 !!! ! 172 !!! REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 173 !!! REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 174 !!! ! 175 !!! REAL(wp), DIMENSION(jpi,jpj) :: zdt_m ! (dt / ice-snow_mass) on T points 176 !!! REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 177 !!! REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! (ice-snow_mass / dt) on U/V points 178 !!! REAL(wp), DIMENSION(jpi,jpj) :: zmf ! coriolis parameter at T points 179 !!! REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 180 !!! ! 181 !!! REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 182 !!! REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 183 !!! REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 184 !!! ! ! ocean surface (ssh_m) if ice is not embedded 185 !!! ! ! ice bottom surface if ice is embedded 186 !!! REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses 187 !!! REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points 188 !!! REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array 189 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points 190 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi, ztauy_oi ! ice-ocean stress at U-V points 165 INTEGER , INTENT(in ) :: kt ! time step 166 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i ! 167 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 168 !! 169 LOGICAL :: ll_u_iterate, ll_viterate ! continue iteration or not 170 ! 171 INTEGER :: ji, jj, jn ! dummy loop indices 172 INTEGER :: jter, i_out, i_inn ! 173 INTEGER :: ji_min, jj_min ! 174 INTEGER :: nn_zebra_vp ! number of zebra steps 175 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp) 176 ! 177 REAL(wp) :: zrhoco ! rau0 * rn_cio 178 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 179 REAL(wp) :: zglob_area ! global ice area for diagnostics 180 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 181 REAL(wp) :: zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 182 REAL(wp) :: zdeltat, zdeltat_star, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 183 REAL(wp) :: zp_deltastar_f 184 REAL(wp) :: zfac, zfac1, zfac2, zfac3 185 REAL(wp) :: zt12U, zt11U, zt22U, zt21U, zt122U, zt121U 186 REAL(wp) :: zt12V, zt11V, zt22V, zt21V, zt122V, zt121V 187 REAL(wp) :: zAA3, zw, zuerr_max, zverr_max 188 ! 189 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice 190 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 191 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! Acceleration term contribution to RHS 192 REAL(wp), DIMENSION(jpi,jpj) :: zmassU_t, zmassV_t ! Mass per unit area divided by time step 193 ! 194 REAL(wp), DIMENSION(jpi,jpj) :: zp_deltastar_t ! P/delta at T points 195 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zet ! Viscosity pre-factors at T points 196 REAL(wp), DIMENSION(jpi,jpj) :: zef ! Viscosity pre-factor at F point 197 ! 198 REAL(wp), DIMENSION(jpi,jpj) :: zmt ! Mass per unit area at t-point 199 REAL(wp), DIMENSION(jpi,jpj) :: zmf ! Coriolis factor (m*f) at t-point 200 REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 201 REAL(wp), DIMENSION(jpi,jpj) :: zu_c, zv_c ! "current" ice velocity (m/s), average of previous two OL iterates 202 REAL(wp), DIMENSION(jpi,jpj) :: zu_cV, zv_cU ! "current" u (v) ice velocity at V (U) point (m/s) 203 REAL(wp), DIMENSION(jpi,jpj) :: zu_b, zv_b ! velocity at previous sub-iterate 204 REAL(wp), DIMENSION(jpi,jpj) :: zuerr, zverr ! absolute U/Vvelocity difference between current and previous sub-iterates 205 ! 206 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 207 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components 208 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 209 ! ! ocean surface (ssh_m) if ice is not embedded 210 ! ! ice bottom surface if ice is embedded 211 REAL(wp), DIMENSION(jpi,jpj) :: zCwU, zCwV ! ice-ocean drag pre-factor (rho*c*module(u)) 212 REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points 213 REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array 214 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_ai, ztauy_ai ! ice-atm. stress at U-V points 215 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_oi_rhsu, ztauy_oi_rhsv ! ice-ocean stress RHS contribution at U-V points 216 REAL(wp), DIMENSION(jpi,jpj) :: zs1_rhsu, zs2_rhsu, zs12_rhsu ! internal stress contributions to RHSU 217 REAL(wp), DIMENSION(jpi,jpj) :: zs1_rhsv, zs2_rhsv, zs12_rhsv ! internal stress contributions to RHSV 218 REAL(wp), DIMENSION(jpi,jpj) :: zf_rhsu, zf_rhsv ! U- and V- components of internal force RHS contributions 219 REAL(wp), DIMENSION(jpi,jpj) :: zrhsu, zrhsv ! U and V RHS 220 REAL(wp), DIMENSION(jpi,jpj) :: zAU, zBU, zCU, zDU, zEU ! Linear system coefficients, U equation 221 REAL(wp), DIMENSION(jpi,jpj) :: zAV, zBV, zCV, zDV, zEV, zFV ! Linear system coefficients, V equation 222 REAL(wp), DIMENSION(jpi,jpj) :: zFU, zFU_prime, zBU_prime ! Rearranged linear system coefficients, U equation 223 REAL(wp), DIMENSION(jpi,jpj) :: zFV, zFV_prime, zBV_prime ! Rearranged linear system coefficients, V equation 191 224 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) 192 225 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 193 !!! ! 194 !!! REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 195 !!! REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 196 197 !!! 198 !!! REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 199 !!! REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 200 !!! REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 201 !!! 202 !!! REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 203 !!! !! --- check convergence 204 !!! REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 205 !!! !! --- diags 206 207 !! --- check convergence 208 209 210 !! --- LSR arrays 211 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice 212 213 REAL(wp), DIMENSION(jpi,jpj) :: zaU , zaV ! ice fraction on U/V points 214 REAL(wp), DIMENSION(jpi,jpj) :: zmassU , zmassV ! ice mass per unit area at U/V points (kg/m2) 215 216 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice_C , zv_ice_C ! ice velocities averaged between before and mid outer sub-iteration 217 REAL(wp), DIMENSION(jpi,jpj) :: uTmp , vTmp ! uTmp = velocity at previous *inner* interation 218 219 REAL(wp), DIMENSION(jpi,jpj) :: z_etat , z_zetat ! VP viscosities at T points 220 REAL(wp), DIMENSION(jpi,jpj) :: z_etaf ! VP eta viscosity at F point 221 222 REAL(wp), DIMENSION(jpi,jpj) :: zs1_rhsu, zs2_rhsu, zs12_rhsu ! Internal stress contribution to RHS of U-equation 223 REAL(wp), DIMENSION(jpi,jpj) :: zs1_rhsv, zs2_rhsv, zs12_rhsv ! Internal stress contribution to RHS of V-equation 224 REAL(wp), DIMENSION(jpi,jpj) :: zf_rhsu, zf_rhsv ! U- and V- components of internal forces 225 226 !!! !! --- SIMIP diags 226 ! 227 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 228 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 229 ! 230 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 231 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 232 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 233 !! --- diags 234 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 235 !! --- SIMIP diags 227 236 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 228 237 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) … … 231 240 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) 232 241 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) 233 234 !! --- Convergence diags 235 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zdiag_meanke(:) 236 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zdiag_resnor(:) 237 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zdiag_veldif(:) 238 242 239 243 !!---------------------------------------------------------------------------------------------------------------------- 240 244 241 245 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 242 243 ! Namelist parameters to outsource (once we all are super happy) 244 nn_nouter_vp = 3 ! Number of outer iterations (mitGCM 10) 245 nn_ninner_vp = 1500 ! Number of inner-loop iterations (mitGCM 1500) 246 rn_deltamin = 2.0e-9 ! minimum delta value following Lemieux and Dupont, GMD 2020 247 ln_visc_VP_smooth = .FALSE. ! zeta viscosity smoothing (if yes, tanh function of Lemieux and Tremblay JGR 2009) 248 249 ! Linear system parameters 250 ln_zebra = .FALSE. ! do two passages for the linear system problem, one every odd j-band, one every even one (recommended by Martin Losch) 251 IF ( ln_zebra ) THEN; nn_zebra_step = 2; ELSE; nn_zebra_step = 1; ENDIF 252 253 rn_relaxfacU = 0.95 ! U-relaxation factor for linear problem (Losch 0.95, Zhang 1.8) 254 rn_relaxfacV = 0.95 ! V-relaxation factor 255 256 rn_uerr_max = 0.8_wp ! maximum error on velocity, above which a forcing error is considered 257 rn_uerr_min = 1.0d-4 ! velocity difference beyond which the linear system procedure is stopped 258 nn_lsr_cvgcheck = 5 ! iterations every each convergence is checked 259 246 260 247 !------------------------------------------------------------------------------! 261 ! 1) Initialization 248 ! 249 ! --- Initialization 250 ! 262 251 !------------------------------------------------------------------------------! 252 253 ! for diagnostics and convergence tests 254 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 255 DO jj = 1, jpj 256 DO ji = 1, jpi 257 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 258 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 259 END DO 260 END DO 261 262 IF ( ln_zebra_vp ) THEN; nn_nzebra_vp = 2; ELSE; nn_nzebra_vp = 1; ENDIF 263 264 nn_nvp = nn_out_vp * nn_inn_vp ! maximum number of iterations 265 263 266 zrhoco = rau0 * rn_cio 264 267 … … 272 275 zs12(:,:) = pstress12_i(:,:) 273 276 274 ! Convergence checks 275 IF( ln_rhg_chkcvg ) THEN 276 277 ! Initialise convergence checks 278 IF( ln_rhg_chkcvg ) THEN 279 280 ! ice area for global mean kinetic energy 277 281 zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) ) ! global ice area (km2) 278 z1_32 = 1._wp / 32._wp279 280 ! --- mean kinetic energy281 ALLOCATE( zdiag_meanke(nn_outer_vp*nn_inner_vp) ) ! mean kinetic energy282 ALLOCATE( zdiag_resnor(nn_outer_vp*nn_inner_vp) ) ! residual norm of the linear system283 ALLOCATE( zdiag_u_idif(nn_outer_vp*nn_inner_vp) ) ! velocity difference between two iterations284 ALLOCATE( zdiag_v_idif(nn_outer_vp*nn_inner_vp) ) ! velocity difference between two iterations285 286 zdiag_resnor(:) = 0._wp287 zdiag_meanke(:) = 0._wp288 zdiag_u_idif(:) = 0._wp289 zdiag_v_idif(:) = 0._wp290 282 291 283 ENDIF 292 284 293 ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 285 ! Landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 286 ! MV: Not working yet... 294 287 IF( ln_landfast_L16 ) THEN ; zkt = rn_lf_tensile 295 288 ELSE ; zkt = 0._wp … … 297 290 298 291 !------------------------------------------------------------------------------! 299 ! 2) Time-independent quantities 292 ! 293 ! --- Time-independent quantities 294 ! 300 295 !------------------------------------------------------------------------------! 301 296 302 297 CALL ice_strength ! strength at T points 303 298 304 !-------- 305 ! F-mask306 !-------- 299 !------------------------------ 300 ! -- F-mask (code from EVP) 301 !------------------------------ 307 302 ! MartinV: 308 303 ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip) 309 ! I don't know exactly where it should be applied inhere310 311 ! ocean/land mask312 DO jj = 1, jpj m1313 DO ji = 1, jpi m1 ! NO vector opt.304 ! I am not sure the same recipe applies here 305 306 ! - ocean/land mask 307 DO jj = 1, jpj - 1 308 DO ji = 1, jpi - 1 314 309 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 315 310 END DO … … 318 313 319 314 ! Lateral boundary conditions on velocity (modify zfmask) 320 DO jj = 2, jpjm1 321 DO ji = fs_2, fs_jpim1 ! vector opt. 315 ! Can be computed once for all, at first time step, for all rheologies 316 DO jj = 2, jpj - 1 317 DO ji = 2, jpi - 1 ! vector opt. 322 318 IF( zfmask(ji,jj) == 0._wp ) THEN 323 319 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & … … 326 322 END DO 327 323 END DO 328 DO jj = 2, jpj m1324 DO jj = 2, jpj - 1 329 325 IF( zfmask(1,jj) == 0._wp ) THEN 330 326 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 331 327 ENDIF 332 328 IF( zfmask(jpi,jj) == 0._wp ) THEN 333 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpi m1,jj,1), umask(jpi,jj-1,1) ) )329 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpi - 1,jj,1), umask(jpi,jj-1,1) ) ) 334 330 ENDIF 335 331 END DO 336 DO ji = 2, jpi m1332 DO ji = 2, jpi - 1 337 333 IF( zfmask(ji,1) == 0._wp ) THEN 338 334 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 339 335 ENDIF 340 336 IF( zfmask(ji,jpj) == 0._wp ) THEN 341 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpj m1,1) ) )337 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpj - 1,1) ) ) 342 338 ENDIF 343 339 END DO … … 345 341 CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp ) 346 342 347 !--------------------------------------------------------------------------------------------------------- 348 ! 3) Time-independent contributions of acceleration, ocean drag, coriolis, atmospheric drag, surface tilt 349 !--------------------------------------------------------------------------------------------------------- 350 ! 351 ! --- Compute all terms & factors independent of velocities, or only depending on velocities at previous time step 352 353 DO jj = 2, jpjm1 354 DO ji = fs_2, fs_jpim1 343 !---------------------------------------------------------------------------------------------------------- 344 ! -- Time-independent pre-factors for acceleration, ocean drag, coriolis, atmospheric drag, surface tilt 345 !---------------------------------------------------------------------------------------------------------- 346 ! Compute all terms & factors independent of velocities, or only depending on velocities at previous time step 347 348 ! sea surface height 349 ! embedded sea ice: compute representative ice top surface 350 ! non-embedded sea ice: use ocean surface for slope calculation 351 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 352 353 DO jj = 2, jpj - 1 354 DO ji = 2, jpi - 1 355 355 356 356 ! Ice fraction at U-V points … … 359 359 360 360 ! Snow and ice mass at U-V points 361 zm 1= ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) )361 zmt(ji,jj) = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) ) 362 362 zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) ) 363 363 zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) ) 364 364 365 zmassU = 0.5_wp * ( zm 1* e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)366 zmassV = 0.5_wp * ( zm 1* e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)365 zmassU = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 366 zmassV = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 367 367 368 ! Mass per unit time368 ! Mass per unit area divided by time step 369 369 zmassU_t(ji,jj) = zmassU * r1_rdtice 370 370 zmassV_t(ji,jj) = zmassV * r1_rdtice … … 379 379 380 380 ! Coriolis factor at T points (m*f) 381 zmf(ji,jj) = zm 1* ff_t(ji,jj)381 zmf(ji,jj) = zmt(ji,jj) * ff_t(ji,jj) 382 382 383 383 ! Wind stress … … 393 393 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 394 394 395 ! switches 395 ! switches 396 396 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 397 397 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF … … 403 403 404 404 !------------------------------------------------------------------------------! 405 ! 4) Start outer loop, update velocities 405 ! 406 ! --- Start outer loop 407 ! 406 408 !------------------------------------------------------------------------------! 407 409 408 zu_ice_C(:,:) = u_ice(:,:) 409 zv_ice_C(:,:) = v_ice(:,:) 410 411 i_iter = 0 412 413 DO i_out = 1, nn_nouter_vp 414 415 i_iter = i_iter + 1 416 417 ! EVP CODE 418 ! l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 419 ! 420 ! convergence test 421 IF( ln_rhg_chkcvg ) THEN 422 DO jj = 1, jpj 423 DO ji = 1, jpi 424 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 425 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 426 END DO 427 END DO 428 ENDIF 429 430 ! Zhang and Hibler 97 time stepping 431 ! identified as the most rapid test with mitGCM and recommendation by Martin Losch 432 zu_ice_C(:,:) = 0.5_wp * ( u_ice(:,:) + zu_ice_C(:,:) ) 433 zv_ice_C(:,:) = 0.5_wp * ( v_ice(:,:) + zv_ice_C(:,:) ) 434 435 !------------------------------------------------------------------------------! 436 ! 5) Right-hand side of linear problem 437 !------------------------------------------------------------------------------! 438 439 ! 410 zu_c(:,:) = u_ice(:,:) 411 zv_c(:,:) = v_ice(:,:) 412 413 jter = 0 414 415 DO i_out = 1, nn_nout_vp 416 417 ! Velocities used in the non linear terms are the average of the past two iterates 418 ! u_it = 0.5 * ( u_{it-1} + u_{it-2}) 419 ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009) 420 zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) ) 421 zv_c(:,:) = 0.5_wp * ( v_ice(:,:) + zv_c(:,:) ) 422 423 !------------------------------------------------------------------------------! 424 ! 425 ! --- Right-hand side (RHS) of the linear problem 426 ! 427 !------------------------------------------------------------------------------! 440 428 ! In the outer loop, one needs to update all RHS terms 441 429 ! with explicit velocity dependencies (viscosities, coriolis, ocean stress) 442 ! 443 444 !-------------------------------------------------------- 445 ! 5.1 Strain rates, viscosities and replacement pressure 446 !-------------------------------------------------------- 430 ! as a function of uc 431 ! 432 433 !------------------------------------------ 434 ! -- Strain rates, viscosities and P/Delta 435 !------------------------------------------ 447 436 448 437 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 449 DO jj = 1, jpj m1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points450 DO ji = 1, jpi m1438 DO jj = 1, jpj - 1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 439 DO ji = 1, jpi - 1 451 440 452 441 ! shear at F points 453 zds(ji,jj) = ( ( zu_ ice_C(ji,jj+1) * r1_e1u(ji,jj+1) - zu_ice_C(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &454 & + ( zv_ ice_C(ji+1,jj) * r1_e2v(ji+1,jj) - zv_ice_C(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &442 zds(ji,jj) = ( ( zu_c(ji,jj+1) * r1_e1u(ji,jj+1) - zu_c(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 443 & + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 455 444 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 456 445 … … 458 447 END DO 459 448 460 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) 461 462 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12463 DO ji = 2, jpi ! no vector loop449 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 450 451 DO jj = 2, jpj - 1 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 452 DO ji = 2, jpi - 1 ! 464 453 465 454 ! shear**2 at T points (doc eq. A16) … … 469 458 470 459 ! divergence at T points 471 zdiv = ( e2u(ji,jj) * zu_ ice_C(ji,jj) - e2u(ji-1,jj) * zu_ice_C(ji-1,jj) &472 & + e1v(ji,jj) * zv_ ice_C(ji,jj) - e1v(ji,jj-1) * zv_ice_C(ji,jj-1) &460 zdiv = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) & 461 & + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) & 473 462 & ) * r1_e1e2t(ji,jj) 474 463 zdiv2 = zdiv * zdiv 475 464 476 465 ! tension at T points 477 zdt = ( ( zu_ ice_C(ji,jj) * r1_e2u(ji,jj) - zu_ice_C(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &478 & - ( zv_ ice_C(ji,jj) * r1_e1v(ji,jj) - zv_ice_C(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &466 zdt = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 467 & - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 479 468 & ) * r1_e1e2t(ji,jj) 480 469 zdt2 = zdt * zdt … … 486 475 zdeltat_star = MAX( zdeltat, rn_delta_min ) 487 476 488 ! IF ( ln_visc_VP_smooth ) THEN ! to code, other regularization from Lemieux and Tremblay (2009) 489 ! zdeltat_star = ... 490 ! ENDIF 477 IF ( ln_smooth_vp ) zdelta_star = zdeltat + rn_delta_min 491 478 492 479 ! P/delta at T-points … … 494 481 495 482 ! Temporary zzt and zet factors at T-points 496 zzt(ji,jj) = zp_deltastar_t * r1_e1e2t(ji,jj)483 zzt(ji,jj) = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj) 497 484 zet(ji,jj) = zzt(ji,jj) * z1_ecc2 498 485 … … 500 487 END DO 501 488 502 CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. ) 503 CALL lbc_lnk( 'icedyn_rhg_vp', zzt , 'T', 1. ) 504 CALL lbc_lnk( 'icedyn_rhg_vp', zet, 'T', 1. ) 505 506 DO jj = 1, jpjm1 507 DO ji = 1, jpim1 489 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 490 491 DO jj = 1, jpj - 1 492 DO ji = 1, jpi - 1 508 493 509 494 ! P/delta* at F points … … 519 504 520 505 !--------------------------------------------------- 521 ! 5.2Ocean-ice drag and Coriolis RHS contributions506 ! -- Ocean-ice drag and Coriolis RHS contributions 522 507 !--------------------------------------------------- 523 508 524 DO jj = 2, jpj m1525 DO ji = fs_2, fs_jpim1509 DO jj = 2, jpj - 1 510 DO ji = 2, jpi - 1 526 511 527 512 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 528 zu_ ice_CV(ji,jj) = 0.25_wp * ( zu_ice_C(ji,jj) + zu_ice_C(ji-1,jj) + zu_ice_C(ji,jj+1) + zu_ice_C(ji-1,jj+1) ) * vmask(ji,jj,1)529 zv_ ice_CU(ji,jj) = 0.25_wp * ( zv_ice_C(ji,jj) + zv_ice_C(ji,jj-1) + zv_ice_C(ji+1,jj) + zv_ice_C(ji+1,jj-1) ) * umask(ji,jj,1)530 513 zu_cV = 0.25_wp * ( zu_c(ji,jj) + zu_c(ji-1,jj) + zu_c(ji,jj+1) + zu_c(ji-1,jj+1) ) * vmask(ji,jj,1) 514 zv_cU = 0.25_wp * ( zv_c(ji,jj) + zv_c(ji,jj-1) + zv_c(ji+1,jj) + zv_c(ji+1,jj-1) ) * umask(ji,jj,1) 515 531 516 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 532 zCwU(ji,jj) = zaU(ji,jj) * zrhoco * SQRT( ( zu_ ice_C (ji,jj) - u_oce (ji,jj) ) * ( zu_ice_C(ji,jj) - u_oce (ji,jj) ) &533 & + ( zv_ ice_CU(ji,jj) - v_oceU(ji,jj) ) * ( zv_ice_CU(ji,jj)- v_oceU(ji,jj) ) )534 zCwV(ji,jj) = zaV(ji,jj) * zrhoco * SQRT( ( zv_ ice_C (ji,jj) - v_oce (ji,jj) ) * ( zv_ice_C(ji,jj) - v_oce (ji,jj) ) &535 & + ( zu_ ice_CV(ji,jj) - u_oceV(ji,jj) ) * ( zu_ice_CV(ji,jj)- u_oceV(ji,jj) ) )517 zCwU(ji,jj) = zaU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) ) & 518 & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 519 zCwV(ji,jj) = zaV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) ) & 520 & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 536 521 537 522 !--- Ocean-ice drag contributions to RHS … … 542 527 ! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done 543 528 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 544 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * zv_ ice_C(ji ,jj) + e1v(ji ,jj-1) * zv_ice_C(ji ,jj-1) ) &545 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_ ice_C(ji+1,jj) + e1v(ji+1,jj-1) * zv_ice_C(ji+1,jj-1) ) )529 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * zv_c(ji ,jj) + e1v(ji ,jj-1) * zv_c(ji ,jj-1) ) & 530 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) ) 546 531 547 532 ! --- V-component of Coriolis Force (energy conserving formulation) 548 533 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 549 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * zu_ ice_C(ji,jj ) + e2u(ji-1,jj ) * zu_ice_C(ji-1,jj ) ) &550 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_ ice_C(ji,jj+1) + e2u(ji-1,jj+1) * zu_ice_C(ji-1,jj+1) ) )534 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * zu_c(ji,jj ) + e2u(ji-1,jj ) * zu_c(ji-1,jj ) ) & 535 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 551 536 552 537 END DO 553 538 END DO 554 539 555 ! ??? lbclnk ??? 556 557 !-------------------------------------- 558 ! 5.3 Internal stress RHS contribution 559 !-------------------------------------- 560 561 ! --- Stress contributions at t-points 540 ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system, 541 ! so there is no need for lbclnk on drag and coriolis 542 543 !------------------------------------- 544 ! -- Internal stress RHS contribution 545 !------------------------------------- 546 547 ! --- Stress contributions at T-points 562 548 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 563 DO ji = 2, jpi ! no vector loop549 DO ji = 2, jpi ! 564 550 565 551 ! sig1 contribution to RHS of U-equation at T-points 566 zs1_rhsu(ji,jj) = zzt(ji,jj) * ( e1v(ji,jj) * zv_ ice_C(ji,jj) - e1v(ji,jj-1) * zv_ice_C(ji,jj-1) - 1.0_wp )552 zs1_rhsu(ji,jj) = zzt(ji,jj) * ( e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) - 1.0_wp ) 567 553 568 554 ! sig2 contribution to RHS of U-equation at T-points 569 zs2_rhsu(ji,jj) = - zet(ji,jj) * ( r1_e1v(ji,jj) * zv_ ice_C(ji,jj) - r1_e1v(ji,jj-1) * zv_ice_C(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)555 zs2_rhsu(ji,jj) = - zet(ji,jj) * ( r1_e1v(ji,jj) * zv_c(ji,jj) - r1_e1v(ji,jj-1) * zv_c(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) 570 556 571 557 ! sig1 contribution to RHS of V-equation at T-points 572 zs1_rhsv(ji,jj) = zzt(ji,jj) * ( e2u(ji,jj) * zu_ice_C(ji,jj) - e2u(ji-1,jj) * zu_ice_C(ji-1,jj) - 1.0_wp )558 zs1_rhsv(ji,jj) = zzt(ji,jj) * ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) - 1.0_wp ) 573 559 574 560 ! sig2 contribution to RHS of V-equation at T-points 575 zs2_rhsv(ji,jj) = zet(ji,jj) * ( r1_e2u(ji,jj) * zu_ ice_C(ji,jj) - r1_e2u(ji-1,jj) * zu_ice_C(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)561 zs2_rhsv(ji,jj) = zet(ji,jj) * ( r1_e2u(ji,jj) * zu_c(ji,jj) - r1_e2u(ji-1,jj) * zu_c(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) 576 562 577 563 END DO 578 564 END DO 579 565 580 ! ??? lbclnk ???566 ! a priori, no lbclnk, because rhsu is only used in the inner domain 581 567 582 568 ! --- Stress contributions at f-points 583 ! Note from MartinV: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing569 ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing 584 570 ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 585 571 586 DO jj = 1, jpj m1587 DO ji = 1, jpi m1572 DO jj = 1, jpj - 1 573 DO ji = 1, jpi - 1 588 574 589 575 ! sig12 contribution to RHS of U equation at F-points 590 zs12_rhsu(ji,jj) = - zef(ji,jj) * ( r1_e2v(ji+1,jj) * zv_ ice_C(ji+1,jj) - r1_e2v(ji,jj) * zv_ice_C(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj)576 zs12_rhsu(ji,jj) = - zef(ji,jj) * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) - r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj) 591 577 592 578 ! sig12 contribution to RHS of V equation at F-points 593 zs12_rhsv(ji,jj) = zef(ji,jj) * ( r1_e1u(ji,jj+1) * zu_ ice_C(ji,jj+1) - r1_e1u(ji,jj) * zu_ice_C(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj)579 zs12_rhsv(ji,jj) = zef(ji,jj) * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) - r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj) 594 580 595 581 END DO 596 582 END DO 597 583 598 ! ??? lbcblnk 599 600 ! --- Internal forces contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) 601 DO jj = 2, jpjm1 602 DO ji = fs_2, fs_jpim1 584 ! a priori, no lbclnk, because rhsu are only used in the inner domain 585 586 ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) 587 ! OPT: merge with next loop and use intermediate scalars for zf_rhsu 588 589 DO jj = 2, jpj - 1 590 DO ji = 2, jpi - 1 603 591 ! --- U component of internal force contribution to RHS at U points 604 592 zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & … … 606 594 & + 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) ) & 607 595 & + 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) ) 596 608 597 ! --- V component of internal force contribution to RHS at V points 609 598 zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * & … … 611 600 & + 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) ) & 612 601 & + 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) ) 602 613 603 END DO 614 604 END DO 615 605 616 ! ??? lbclnk ???617 618 606 !--------------------------- 619 ! 5.4Sum RHS contributions607 ! -- Sum RHS contributions 620 608 !--------------------------- 621 DO jj = 2, jpjm1 622 DO ji = fs_2, fs_jpim1 609 ! 610 ! OPT: could use intermediate scalars to reduce memory access 611 DO jj = 2, jpj - 1 612 DO ji = 2, jpi - 1 623 613 624 614 ! still miss ice ocean stress and acceleration contribution … … 629 619 END DO 630 620 631 ! ??? lbclnk ???621 ! inner domain calculations -> no lbclnk 632 622 633 !---------------------------------------------------------------------------------------! 634 ! 6) Linear system matrix 635 !---------------------------------------------------------------------------------------! 623 !---------------------------------------------------------------------------------------! 624 ! 625 ! --- Linear system matrix 626 ! 627 !---------------------------------------------------------------------------------------! 636 628 637 629 ! Linear system matrix contains all implicit contributions … … 641 633 ! AU * u_{i-1,j} + BU * u_{i,j} + CU * u_{i+1,j} 642 634 ! = DU * u_{i,j-1} + EU * u_{i,j+1} + RHS (! my convention, not the same as ZH97 ) 643 ! 644 !-------------------------------------- 645 ! 6.1 Internal forces LHS contribution 646 !-------------------------------------- 647 ! 648 ! Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ? 649 ! Note 2: "T" factor calculations can be optimized by putting things out of the loop 635 636 ! MV Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ? 637 ! MV Note 2: "T" factor calculations can be optimized by putting things out of the loop 650 638 ! only zzt and zet are iteration-dependent, other only depend on scale factors 651 ! 652 ! --- U-component 653 ! 654 ! "T" factors (intermediate results) 655 ! 656 zfac = 0.5_wp * r1_e1e2u(ji,jj) 657 zfac1 = zfac * e2u(ji,jj) 658 zfac2 = zfac * r1_e2u(ji,jj) 659 zfac3 = 2._wp * zfac * r1_e1u(ji,jj) 660 661 zt12U = - zfac1 * zzt(ji+1,jj) 662 zt11U = zfac1 * zzt(ji,jj) 663 664 zt22U = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 665 zt21U = zfac2 * zet(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) 666 667 zt122U = - zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) 668 zt121U = zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 669 670 ! 671 ! Linear system coefficients 672 ! 673 zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 674 zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj) 675 zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 676 677 zDU(ji,jj) = zt121U * r1_e1u(ji,jj-1) 678 zEU(ji,jj) = zt122U * r1_e1u(ji,jj+1) 679 680 ! 681 ! --- V-component 682 ! 683 ! "T" factors (intermediate results) 684 ! 685 zfac = 0.5_wp * r1_e1e2v(ji,jj) 686 zfac1 = zfac * e2v(ji,jj) 687 zfac2 = zfac * r1_e1v(ji,jj) 688 zfac3 = 2._wp * zfac * r1_e2v(ji,jj) 689 690 zt12V = - zfac1 * zzt(ji,jj+1) 691 zt11V = zfac1 * zzt(ji,jj) 692 693 zt22V = zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 694 zt21V = - zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) 695 696 zt122V = zfac3 * zef(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) 697 zt121V = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 698 699 ! 700 ! Linear system coefficients 701 ! 702 zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1) 703 zBV(ji,jj) = ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj) 704 zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1) 705 706 zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here 707 zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj) 708 709 !--------------------------------- 710 ! Ocean-ice drag LHS contribution 711 !--------------------------------- 712 zBU(i,j) = zBU(i,j) + zCwU(ji,jj) 713 714 !------------------------------- 715 ! Acceleration LHS contribution 716 !------------------------------- 717 zBU(i,j) = zBU(i,j) + zmassU_t(ji,jj) 639 640 DO ji = 2, jpi - 1 ! internal domain do loop 641 DO jj = 2, jpj - 1 642 643 !------------------------------------- 644 ! -- Internal forces LHS contribution 645 !------------------------------------- 646 ! 647 ! --- U-component 648 ! 649 ! "T" factors (intermediate results) 650 ! 651 zfac = 0.5_wp * r1_e1e2u(ji,jj) 652 zfac1 = zfac * e2u(ji,jj) 653 zfac2 = zfac * r1_e2u(ji,jj) 654 zfac3 = 2._wp * zfac * r1_e1u(ji,jj) 655 656 zt12U = - zfac1 * zzt(ji+1,jj) 657 zt11U = zfac1 * zzt(ji,jj) 658 659 zt22U = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 660 zt21U = zfac2 * zet(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) 661 662 zt122U = - zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) 663 zt121U = zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 664 665 ! 666 ! Linear system coefficients 667 ! 668 zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 669 zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj) 670 zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 671 672 zDU(ji,jj) = zt121U * r1_e1u(ji,jj-1) 673 zEU(ji,jj) = zt122U * r1_e1u(ji,jj+1) 674 675 ! 676 ! --- V-component 677 ! 678 ! "T" factors (intermediate results) 679 ! 680 zfac = 0.5_wp * r1_e1e2v(ji,jj) 681 zfac1 = zfac * e2v(ji,jj) 682 zfac2 = zfac * r1_e1v(ji,jj) 683 zfac3 = 2._wp * zfac * r1_e2v(ji,jj) 684 685 zt12V = - zfac1 * zzt(ji,jj+1) 686 zt11V = zfac1 * zzt(ji,jj) 687 688 zt22V = zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 689 zt21V = - zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) 690 691 zt122V = zfac3 * zef(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) 692 zt121V = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 693 694 ! 695 ! Linear system coefficients 696 ! 697 zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1) 698 zBV(ji,jj) = ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj) 699 zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1) 700 701 zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here 702 zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj) 703 704 !----------------------------------------------------- 705 ! -- Ocean-ice drag and acceleration LHS contribution 706 !----------------------------------------------------- 707 zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 708 zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 709 710 END DO 711 END DO 718 712 719 713 !------------------------------------------------------------------------------! 720 ! 6) inner loop: solve linear system, check convergence 714 ! 715 ! --- Inner loop: solve linear system, check convergence 716 ! 721 717 !------------------------------------------------------------------------------! 722 718 723 719 ! Inner loop solves the linear problem .. requires 1500 iterations 724 doIterate4u = .TRUE. 725 doIterate4v = .TRUE. 726 727 DO i_inner = 1, nn_ninner_vp ! inner loop iterations 728 729 ! initialize residual ... 730 731 ... 732 733 IF ( doIterate4u .OR. doIterate4v ) THEN 734 735 zu_p(:,:) = u_ice(:,:) ! previous sub-iteration value 736 zv_p(:,:) = v_ice(:,:) 737 738 ! ---------------------------- ! 739 IF ( doIterate4u ) THEN ! --- Solve for u-velocity --- ! 740 ! ---------------------------- ! 720 ll_u_iterate = .TRUE. 721 ll_v_iterate = .TRUE. 722 723 DO i_inn = 1, nn_ninn_vp ! inner loop iterations 724 725 ! mitgcm computes initial value of residual 726 jter = jter + 1 727 l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 728 729 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN 730 731 732 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 733 zv_b(:,:) = v_ice(:,:) 734 735 ! ---------------------------- ! 736 IF ( ll_u_iterate ) THEN ! --- Solve for u-velocity --- ! 737 ! ---------------------------- ! 741 738 742 739 ! what follows could be subroutinized... 743 740 744 DO k = 1, nn_zebra_step ! "zebra" loop 745 746 IF ( k == 1 ); THEN; j_min = 2; ELSE; j_min = 3 747 748 zFU(:,:) = 0._wp 749 750 DO jj = j_min, jpjm1, nn_zebra_step 741 DO jn = 1, nn_nzebra_vp ! "zebra" loop (! red-black-sor!!! ) 742 743 ! OPT: could be even better optimized with a true red-black SOR 744 745 IF ( jn == 1 ) THEN ; jj_min = 2 746 ELSE ; jj_min = 3 747 ENDIF 748 749 zFU(:,:) = 0._wp 750 zFU_prime(:,:) = 0._wp 751 zBU_prime(:,:) = 0._wp 752 753 DO jj = jj_min, jpj - 1, nn_nzebra_vp 751 754 752 755 !------------------------------------------- 753 ! 1)Tridiagonal system solver for each row756 ! -- Tridiagonal system solver for each row 754 757 !------------------------------------------- 758 ! 759 ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan 760 ! 755 761 ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 756 762 763 ! - Right-hand side of tridiagonal system (zFU) 764 DO ji = 2, jpi - 1 765 766 ! boundary condition substitution 767 ! see Zhang and Hibler, 1997, Appendix B 768 ! MV NOTE possibly not fully appropriate 769 zAA3 = 0._wp 770 IF ( ji == 2 ) zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 771 IF ( ji == jpi - 1 ) zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj) 772 773 ! right hand side 774 zFU(ji,jj) = ( zrhsu(ji,jj) & ! right-hand side terms 775 & + zAA3 ! boundary condition translation 776 & + zDU(ji,jj) * u_ice(ji,jj-1) ! internal force, j-1 777 & + zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1 778 779 END DO 780 781 ! - Thomas Algorithm 782 ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 783 ! Forward sweep 784 DO ji = 3, jpi - 1 785 zw = zAU(ji,jj) / zBU(ji-1,jj) 786 zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 787 zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 788 END DO 789 790 ! Backward sweep 791 ! MV I don't see how this will be reproducible 792 u_ice(jpi-1,jj) = zFU_prime(jpi-1,jj) / zBU_prime(jpi-1,jj) * umask(jpi-1,jj,1) ! do last row first 793 DO ji = jpi-2 , 2, -1 ! all other rows ! MV OPT: could be turned into forward loop (by substituting ji) 794 u_ice(ji,jj) = zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji,jj+1) * umask(ji,jj,1) / zBU_prime(ji,jj) ! 795 END DO 796 797 !--------------- 798 ! -- Relaxation 799 !--------------- 800 DO ji = 2, jpi - 1 801 u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) 802 END DO 803 804 END DO ! jj 805 806 END DO ! zebra loop 807 808 ENDIF ! ll_u_iterate 809 810 ! ! ---------------------------- ! 811 IF ( ll_v_iterate ) THEN ! --- Solve for V-velocity --- ! 812 ! ! ---------------------------- ! 813 814 ! MV OPT: what follows could be subroutinized... 815 816 DO jn = 1, nn_nzebra_vp ! "zebra" loop 817 818 IF ( jn == 1 ) THEN ; ji_min = 2 819 ELSE ; ji_min = 3 820 ENDIF 821 822 zFV(:,:) = 0._wp 823 zFV_prime(:,:) = 0._wp 824 zBV_prime(:,:) = 0._wp 825 826 DO ji = ji_min, jpi - 1, nn_nzebra_vp 827 828 !!! It is intentional to have a ji then jj loop for V-velocity 829 !!! ZH97 explain it is critical for convergence speed 830 831 !------------------------------------------- 832 ! -- Tridiagonal system solver for each row 833 !------------------------------------------- 834 ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F 835 757 836 ! --- Right-hand side of tridiagonal system (zFU) 758 DO j i = fs_2, fs_jpim1759 760 ! boundary condition substitution (check whetheris correctly applied !!!)837 DO jj = 2, jpj - 1 838 839 ! boundary condition substitution (check it is correctly applied !!!) 761 840 ! see Zhang and Hibler, 1997, Appendix B 762 841 zAA3 = 0._wp 763 IF ( j i .EQ. fs_2 ) zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj)764 IF ( j i .EQ. fs_jpim1 ) zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj)842 IF ( jj .EQ. 2 ) zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1) 843 IF ( jj .EQ. jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1) 765 844 766 845 ! right hand side 767 zF U(ji,jj) = zrhsu(ji,jj) &! right-hand side terms768 & + zAA3 ! boundary condition translation769 & + zDU(ji,jj) * u_ice(ji,jj-1) ! internal force, j-1770 & + zEU(ji,jj) * u_ice(ji,jj+1) ! internal force, j+1846 zFV(ji,jj) = ( zrhsv(ji,jj) & ! right-hand side terms 847 & + zAA3 ! boundary condition translation 848 & + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1 849 & + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1 771 850 772 zFU(ji,jj) = zFU(ji,jj) * umask(ji,jj) ! mask - useful ?773 774 851 END DO 775 852 … … 777 854 ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 778 855 ! Forward sweep 779 DO j i = fs_2 + 1, fs_jpim1780 zw = zA U(ji,jj) / zBU(ji-1,jj)781 zB U(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj)782 zF U(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj)856 DO jj = 3, jpj - 1 857 zw = zAV(ji,jj) / zBV(ji,jj-1) 858 zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 859 zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 783 860 END DO 784 861 785 862 ! Backward sweep 786 u_ice(fs_jpim1,jj) = zFU(ji) / zBU(ji)! last row787 DO j i = fs_jpim - 1, fs_2, -1 ! can be turned into forward row by substituting jiif needed788 u_ice(ji,jj) = zFU(ji) - zCU(ji,jj) * u_ice(ji,jj+1) / zBU(ji,jj)863 v_ice(ji,jpj-1) = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1) ! last row 864 DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed 865 v_ice(ji,jj) = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj) / zBV_prime(ji,jj) 789 866 END DO 790 867 791 868 !--------------- 792 ! 2)Relaxation869 ! -- Relaxation 793 870 !--------------- 794 DO j i = fs_2, fs_jpim1795 u_ice(ji,jj) = zu_p(ji,jj) + rn_relaxfacU * ( u_ice(ji,jj) - zu_p(ji,jj) )871 DO jj = 2, jpj - 1 872 v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) ) 796 873 END DO 797 874 798 END DO ! j j875 END DO ! ji 799 876 800 877 END DO ! zebra loop 801 802 CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1. ) 803 804 ENDIF ! doIterate4u 805 806 ! ---------------------------- ! 807 IF ( doIterate4v ) THEN ! --- Solve for V-velocity --- ! 808 ! ---------------------------- ! 809 810 ! what follows could be subroutinized... 811 812 DO k = 1, nn_zebra_step ! "zebra" loop 813 814 IF ( k == 1 ); THEN; i_min = 2; ELSE; i_min = 3 815 816 zFV(:,:) = 0._wp 817 818 DO ji = i_min, fs_jpim1, nn_zebra_step 819 820 !!! HERE WE HAVE ji then jj loops - I don't know what this implies 821 !!! But ZH explain it is critical for convergence speed 822 823 !------------------------------------------- 824 ! 1) Tridiagonal system solver for each row 825 !------------------------------------------- 826 ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F 827 828 ! --- Right-hand side of tridiagonal system (zFU) 829 DO jj = 2, jpjm1 830 831 ! boundary condition substitution (check it is correctly applied !!!) 832 ! see Zhang and Hibler, 1997, Appendix B 833 zAA3 = 0._wp 834 IF ( jj .EQ. 2 ) zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1) 835 IF ( jj .EQ. jpjm1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1) 836 837 ! right hand side 838 zFV(ji,jj) = zrhsv(ji,jj) & ! right-hand side terms 839 & + zAA3 ! boundary condition translation 840 & + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1 841 & + zEV(ji,jj) * v_ice(ji+1,jj) ! internal force, j+1 842 843 zFV(ji,jj) = zFV(ji,jj) * vmask(ji,jj) ! mask - useful ? 844 845 END DO 846 847 ! --- Thomas Algorithm 848 ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 849 ! Forward sweep 850 DO jj = 3, jpjm1 851 zw = zAV(ji,jj) / zBV(ji,jj-1) 852 zBV(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 853 zFV(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 854 END DO 855 856 ! Backward sweep 857 v_ice(ji,jpjm1) = zFV(ji) / zBV(ji) ! last row 858 DO jj = jpjm1 - 1, 2, -1 ! can be turned into forward row by substituting jj if needed 859 v_ice(ji,jj) = zFV(ji) - zCV(ji,jj) * v_ice(ji,jj+1) / zBV(ji,jj) 860 END DO 861 862 !--------------- 863 ! 2) Relaxation 864 !--------------- 865 DO jj = 2, jpjm1 866 v_ice(ji,jj) = zv_p(ji,jj) + rn_relaxfacV * ( v_ice(ji,jj) - zv_p(ji,jj) ) 867 END DO 868 869 END DO ! ji 870 871 END DO ! zebra loop 872 873 CALL lbc_lnk( 'icedyn_rhg_vp', v_ice, 'V', -1. ) 874 875 ENDIF ! doIterate4v 878 879 ENDIF ! ll_v_iterate 880 881 IF ( ( ll_u_iterate .OR. ll_v_iterate ) .OR. jter == nn_nvp ) CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 876 882 877 !-------------------------------------------------- 878 ! Continue the loop - based on velocity difference879 !-------------------------------------------------- 883 !-------------------------------------------------------------------------------------- 884 ! -- Check convergence based on maximum velocity difference, continue or stop the loop 885 !-------------------------------------------------------------------------------------- 880 886 881 887 !------ 882 888 ! on U 883 889 !------ 884 IF ( doIterate4u .AND. MOD ( i_inner, nn_lsr_cvgcheck ) == 0 ) THEN 885 886 zuerr(:,:) = 0._wp 887 DO jj = 2, jpjm1 888 DO ji = 2, jpim1 889 zuerr(ji,jj) = ( u_ice(ji,jj) - zu_p(ji,jj) ) * umask(ji,jj) 890 ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 891 ! then we must optimize the use of the mpp_max, which is prohibitive 892 893 IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN 894 895 ! - Maximum U-velocity difference 896 zuerr(:,:) = 0._wp 897 DO jj = 2, jpj - 1 898 DO ji = 2, jpi - 1 899 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) ! * zmask15 ---> MV TEST mask at 15% concentration 890 900 END DO 891 901 END DO 892 893 902 zuerr_max = MAXVAL( zuerr ) 894 CALL mpp_max( 'icedyn_rhg_evp', zuerr_max ) ! max over the global domain 895 896 897 IF ( ln_rhg_chkcvg ) zdiag_u_idif(i_iter) = zuerr_max 898 899 ! "safeguard against bad forcing" -> stop relaxation 900 IF ( ( i_inner > 1 ) && zuerr_max > rn_uerr_max ) THEN 901 !---> multiply relaxation factor above to stop relaxation 902 zrelaxation_switchU = 0._wp 903 CALL mpp_max( 'icedyn_rhg_evp', zuerr_max ) ! max over the global domain - damned! 904 905 ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 906 IF ( i_inn > 1 && zuerr_max > rn_uerr_max_vp ) THEN 907 IF ( lwp ) " VP rheology error was too large : ", zu_err_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 908 ll_u_iterate = .FALSE. 903 909 ENDIF 904 910 905 IF ( zuerr_max .LT. rn_uerr_min ) THEN 906 doIterate4u = .FALSE. 907 IF ( lwp ) " VP rheology finished in outer U-iteration ", i_outer, " after ", i_inner, " iterations " 911 ! - Stop if error small enough 912 IF ( zuerr_max < rn_uerr_min_vp ) THEN 913 IF ( lwp ) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! " 914 ll_u_iterate = .FALSE. 908 915 ENDIF 909 910 ENDIF 916 917 ENDIF ! ll_u_iterate 911 918 912 919 !------ … … 914 921 !------ 915 922 916 IF ( doIterate4v .AND. MOD ( i_inner, nn_lsr_cvgcheck ) == 0 ) THEN 917 923 IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN 924 925 ! - Maximum V-velocity difference 918 926 zverr(:,:) = 0._wp 919 DO jj = 2, jpj m1920 DO ji = 2, jpi m1921 zverr(ji,jj) = ( v_ice(ji,jj) - zv_p(ji,jj) ) * vmask(ji,jj)927 DO jj = 2, jpj - 1 928 DO ji = 2, jpi - 1 929 zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1) 922 930 END DO 923 931 END DO 924 932 925 933 zverr_max = MAXVAL( zverr ) 926 CALL mpp_max( 'icedyn_rhg_evp', zverr_max ) ! max over the global domain 927 928 IF ( ln_rhg_chkcvg ) zdiag_v_idif(i_iter) = zverr_max 929 930 ! "safeguard against bad forcing" -> stop relaxation 931 IF ( ( i_inner > 1 ) && zverr_max > rn_uerr_max ) THEN 932 !---> multiply relaxation factor above to stop relaxation 933 zrelaxation_switchV = 0._wp 934 CALL mpp_max( 'icedyn_rhg_evp', zverr_max ) ! max over the global domain - damned! 935 936 ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 937 IF ( i_inn > 1 && zverr_max > rn_uerr_max_vp ) THEN 938 IF ( lwp ) " VP rheology error was too large : ", zv_err_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 939 ll_v_iterate = .FALSE. 934 940 ENDIF 935 941 936 IF ( zverr_max .LT. rn_verr_min ) THEN 937 doIterate4v = .FALSE. 938 IF ( lwp ) " VP rheology finished in outer V-iteration ", i_outer, " after ", i_inner, " iterations " 942 ! - Stop if error small enough 943 IF ( zverr_max < rn_verr_min ) THEN 944 IF ( lwp ) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! " 945 ll_v_iterate = .FALSE. 939 946 ENDIF 940 947 941 ENDIF 942 943 !----------------------------------------------------------------------------- 944 ! Convergence diagnostics: residual norm and mean kinetic energy at inner iteration 945 !----------------------------------------------------------------------------- 946 947 IF( ln_rhg_chkcvg ) THEN 948 949 ! residual norm, velocity difference 950 ! ... zdiag_resnor 951 952 ! mean kinetic energy 953 ztotke = 0._wp 954 DO jj = 2, jpjm1 955 DO ji = 2, jpim1 956 ztotke = ztotke + ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * at_i(ji,jj) & 957 & ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * z1_32 958 END DO 959 END DO 960 961 zdiag_meanke(i_iter) = ztotke / zglob_area 962 963 ENDIF 964 965 ENDIF ! --- end doIterate4u or doIterate4v 948 ENDIF ! ll_v_iterate 949 950 !--------------------------------------------------------------------------------------- 951 ! 952 ! --- Calculate extra convergence diagnostics and save them 953 ! 954 !--------------------------------------------------------------------------------------- 955 IF( ln_rhg_chkcvg ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 956 & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 957 958 959 ENDIF ! --- end ll_u_iterate or ll_v_iterate 966 960 967 961 END DO ! i_inn, end of inner loop 968 962 969 963 !------------------------------------------------------------------------------! 964 ! 970 965 ! 6) Mask final velocities 966 ! 971 967 !------------------------------------------------------------------------------! 972 968 973 969 END ! End of outer loop (i_out) ============================================================================================= 974 970 975 !976 971 !------------------------------------------------------------------------------! 977 ! 7) Recompute delta, shear and div (inputs for mechanical redistribution) 972 ! 973 ! --- Recompute delta, shear and div (inputs for mechanical redistribution) 974 ! 978 975 !------------------------------------------------------------------------------! 979 DO jj = 1, jpjm1 980 DO ji = 1, jpim1 976 ! 977 ! OPT: subroutinize ? 978 979 DO jj = 1, jpj - 1 980 DO ji = 1, jpi - 1 981 981 982 982 ! shear at F points … … 988 988 END DO 989 989 990 DO jj = 2, jpj m1991 DO ji = 2, jpi m1 ! no vector loop990 DO jj = 2, jpj - 1 991 DO ji = 2, jpi - 1 ! 992 992 993 993 ! tension**2 at T points … … 1017 1017 END DO 1018 1018 END DO 1019 1019 1020 CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 1021 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 1020 1022 1021 1023 ! --- Store the stress tensor for the next time step --- ! 1022 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 1024 ! MV OPT: Are they computed at the end of the tucking fime step ??? 1025 ! MV Is stress at previous time step needed for VP, normally no, because they equation is not tucking fime dependent!!! 1026 ! 1023 1027 pstress1_i (:,:) = zs1 (:,:) 1024 1028 pstress2_i (:,:) = zs2 (:,:) … … 1027 1031 1028 1032 !------------------------------------------------------------------------------! 1029 ! 8) diagnostics 1033 ! 1034 ! --- Diagnostics 1035 ! 1030 1036 !------------------------------------------------------------------------------! 1031 ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 1037 ! MV OPT: subroutinize 1038 ! 1039 ! --- Ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 1032 1040 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 1033 1041 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN … … 1046 1054 ENDIF 1047 1055 1048 ! --- divergence, shear and strength --- !1056 ! --- Divergence, shear and strength --- ! 1049 1057 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 1050 1058 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 1051 1059 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 1052 1060 1053 ! --- stress tensor --- !1061 ! --- Stress tensor --- ! 1054 1062 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 1055 1063 ! 1056 1064 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 1057 1065 ! 1058 DO jj = 2, jpjm1 1059 DO ji = 2, jpim1 1066 DO jj = 2, jpj - 1 1067 DO ji = 2, jpi - 1 1068 1060 1069 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 1061 1070 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & … … 1066 1075 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 1067 1076 1068 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)1069 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)1070 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress1071 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))1072 1077 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 1073 1078 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 1074 1079 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 1080 1075 1081 END DO 1076 1082 END DO 1083 1077 1084 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 1078 1085 ! … … 1086 1093 1087 1094 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 1095 1088 1096 ENDIF 1089 1097 1090 ! --- SIMIP --- !1098 ! --- SIMIP, terms of tendency for momentum equation --- ! 1091 1099 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 1092 1100 & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1093 1101 ! 1102 !!!!!!!!! ATTENTION LA FORCE INTERNE DOIT ETTRE RECALCIULEEE ICCI !!!!!!!!!!!!!!!! 1094 1103 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & 1095 1104 & zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) … … 1109 1118 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 1110 1119 ! 1111 DO jj = 2, jpj m11112 DO ji = 2, jpi m11120 DO jj = 2, jpj - 1 1121 DO ji = 2, jpi - 1 1113 1122 ! 2D ice mass, snow mass, area transport arrays (X, Y) 1114 1123 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) … … 1143 1152 ENDIF 1144 1153 1145 ! 1146 ! --- convergence diagnostics --- ! 1154 ! --- Convergence diagnostics --- ! 1147 1155 IF( ln_rhg_chkcvg ) THEN 1148 1156 1149 ! still miss residual norm 1150 1151 IF ( iom_use('mean_tke_ice') ) THEN 1152 CALL iom_put( 'mean_tke_ice', zdiag_meanke ) 1153 DEALLOCATE( zdiag_meanke ) 1154 ENDIF 1155 1156 IF ( iom_use('u_veldif') ) THEN 1157 CALL iom_put ( zdiag_u_idif(i_iter) ) 1158 DEALLOCATE( zdiag_u_idif ) 1159 ENDIF 1160 1161 IF ( iom_use('v_veldif') ) THEN 1162 CALL iom_put ( zdiag_v_idif(i_iter) ) 1163 DEALLOCATE( zdiag_v_idif ) 1164 ENDIF 1165 1157 IF( iom_use('uice_cvg') ) THEN 1158 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & 1159 & ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 1160 ENDIF 1161 1166 1162 ENDIF ! ln_rhg_chkcvg 1167 1163 1168 ! EVP convergence diag1169 ! IF( iom_use('uice_cvg') ) THEN1170 !1171 !1172 ! ! EVP-diag1173 ! IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )1174 ! CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &1175 ! & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )1176 ! ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )1177 ! CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &1178 ! & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )1179 ! ENDIF1180 1181 ENDIF1182 ENDIF1183 1164 ! 1184 1165 DEALLOCATE( zmsk00, zmsk15 ) 1185 1186 1166 1187 1167 END SUBROUTINE ice_dyn_rhg_vp 1188 1168 1169 1170 1171 SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, & 1172 & prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 1173 1174 ! CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 1175 ! & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 1176 !!---------------------------------------------------------------------- 1177 !! *** ROUTINE rhg_cvg_vp *** 1178 !! 1179 !! ** Purpose : check convergence of VP ice rheology 1180 !! 1181 !! ** Method : create a file ice_cvg.nc containing a few convergence diagnostics 1182 !! This routine is called every sub-iteration, so it is cpu expensive 1183 !! 1184 !! Calculates / stores 1185 !! - maximum absolute U-V difference (uice_cvg, u_dif, v_dif, m/s) 1186 !! - residuals in U, V and UV-mean taken as square-root of area-weighted mean square residual (u_res, v_res, vel_res, N/m2) 1187 !! - mean kinetic energy (mke_ice, J/m2) 1188 !! 1189 !! 1190 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 1191 !!---------------------------------------------------------------------- 1192 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1193 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pmt ! now velocity and mass per unit area 1194 REAL(wp), INTENT(in) :: puerr_max, pverr_max ! absolute mean velocity difference 1195 REAL(wp), INTENT(in) :: pglob_area ! global ice area 1196 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients 1197 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV 1198 !! 1199 INTEGER :: it, idtime, istatus 1200 INTEGER :: ji, jj ! dummy loop indices 1201 REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, zmke, zu, zv ! local scalars 1202 REAL(wp), DIMENSION(jpi,jpj) :: zu_res(:,:), zv_res(:,:), zvel2(:,:) ! local arrays 1203 1204 CHARACTER(len=20) :: clname 1205 :: zres ! check convergence 1206 !!---------------------------------------------------------------------- 1207 1208 ! create file 1209 IF( kt == nit000 .AND. kiter == 1 ) THEN 1210 ! 1211 IF( lwp ) THEN 1212 WRITE(numout,*) 1213 WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 1214 WRITE(numout,*) '~~~~~~~~~~~' 1215 ENDIF 1216 ! 1217 IF( lwm ) THEN 1218 clname = 'ice_cvg.nc' 1219 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 1220 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 1221 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 1222 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid_ucvg ) 1223 istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1224 istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1225 istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1226 istatus = NF90_DEF_VAR( ncvgid, 'u_dif', NF90_DOUBLE , (/ idtime /), nvarid_udif ) 1227 istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE , (/ idtime /), nvarid_vdif ) 1228 istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke ) 1229 istatus = NF90_ENDDEF(ncvgid) 1230 ENDIF 1231 ! 1232 ENDIF 1233 1234 ! time 1235 it = ( kt - 1 ) * kitermax + kiter 1236 1237 ! --- Max absolute velocity difference with previous iterate (zveldif) 1238 ! EVP code 1239 ! IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 1240 ! zveldif = 0._wp 1241 ! ELSE 1242 ! DO jj = 1, jpj 1243 ! DO ji = 1, jpi 1244 ! zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1245 ! & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 1246 ! END DO 1247 ! END DO 1248 ! zveldif = MAXVAL( zres ) 1249 ! CALL mpp_max( 'icedyn_rhg_vp', zveldif ) ! max over the global domain 1250 ! ENDIF 1251 ! VP code 1252 zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code 1253 ! if puerrmask and pverrmax are masked at 15% (TEST) 1254 1255 ! -- Mean residual (N/m^2), zu_res_mean 1256 ! Here we take the residual of the linear system (N/m^2), 1257 ! We define it as in mitgcm: square-root of area-weighted mean square residual 1258 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1259 ! i.e., how close we are to a solution 1260 DO jj = 2, jpj - 1 1261 DO ji = 2, jpi - 1 1262 zu_res(ji,jj) = zrhsu(ji,jj) + zDU(ji,jj) * pu(ji,jj-1) + zEU(ji,jj) * pu(ji,jj+1) & 1263 & - zAU(ji,jj) * pu(ji-1,jj) - zBU(ji,jj) * pu(ji,jj) - zCU(ji,jj) * pu(ji+1,jj) 1264 1265 zv_res(ji,jj) = zrhsv(ji,jj) + zDV(ji,jj) * pv(ji-1,jj) + zEV(ji,jj) * pv(ji+1,jj) & 1266 & - zAV(ji,jj) * pv(ji,jj-1) - zBV(ji,jj) * pv(ji,jj) - zCV(ji,jj) * pv(ji,jj+1) 1267 END DO 1268 END DO 1269 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:) ) 1270 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:) ) 1271 zu_res_mean = SQRT( zu_resmean / pglob_area ) 1272 zv_res_mean = SQRT( zv_resmean / pglob_area ) 1273 zvelres = MEAN( zu_res_mean, zv_res_mean ) 1274 1275 ! -- Global mean kinetic energy per unit area (J/m2) 1276 DO jj = 2, jpj - 1 1277 DO ji = 2, jpi - 1 1278 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1279 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1280 zvel2(:,:) = zu*zu + zv*zv ! square of ice velocity at T-point 1281 END DO 1282 END DO 1283 1284 zmke = 0.5_wp * globsum( 'ice_rhg_vp', zmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) 1285 1286 ! ! ==================== ! 1287 1288 IF( lwm ) THEN 1289 ! write variables 1290 istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) ) 1291 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) 1292 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) 1293 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) 1294 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) 1295 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) 1296 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) 1297 ! close file 1298 IF( kt == nitend ) istatus = NF90_CLOSE( ncvgid ) 1299 ENDIF 1300 1301 END SUBROUTINE rhg_cvg_vp 1302 1303 1189 1304 1190 1305 SUBROUTINE rhg_vp_rst( cdrw, kt )
Note: See TracChangeset
for help on using the changeset viewer.
