Changeset 15531
- Timestamp:
- 2021-11-23T19:09:14+01:00 (2 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_rhg_vp.F90 (modified) (70 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_rhg_vp.F90
r15519 r15531 41 41 PUBLIC ice_dyn_rhg_vp ! called by icedyn_rhg.F90 42 42 43 44 LOGICAL :: lp_zebra_vp =.TRUE. 43 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp) 44 LOGICAL :: lp_zebra_vp =.TRUE. ! activate zebra (solve the linear system problem every odd j-band, then one every even one) 45 45 REAL(wp) :: zrelaxu_vp=0.95 ! U-relaxation factor (MV: can probably be merged with V-factor once ok) 46 46 REAL(wp) :: zrelaxv_vp=0.95 ! V-relaxation factor 47 47 REAL(wp) :: zuerr_max_vp=0.80 ! maximum velocity error, above which a forcing error is considered and solver is stopped 48 REAL(wp) :: zuerr_min_vp=1.e-04 48 REAL(wp) :: zuerr_min_vp=1.e-04 ! minimum velocity error, beyond which convergence is assumed 49 49 50 50 !! for convergence tests 51 51 INTEGER :: ncvgid ! netcdf file id 52 INTEGER :: nvarid_ures 53 INTEGER :: nvarid_vres 54 INTEGER :: nvarid_velres 55 INTEGER :: nvarid_udif 56 INTEGER :: nvarid_vdif 57 INTEGER :: nvarid_veldif 52 INTEGER :: nvarid_ures, nvarid_vres, nvarid_velres 53 INTEGER :: nvarid_uerr_max, nvarid_verr_max, nvarid_velerr_max 54 INTEGER :: nvarid_umad, nvarid_vmad, nvarid_velmad 55 INTEGER :: nvarid_umad_outer, nvarid_vmad_outer, nvarid_velmad_outer 58 56 INTEGER :: nvarid_mke 59 INTEGER :: nvarid_ures_xy, nvarid_vres_xy60 57 61 58 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice … … 88 85 !! 89 86 !! f1(u) = g1(v) 90 !! f2(v) = g2( v)87 !! f2(v) = g2(u) 91 88 !! 92 89 !! The right-hand side (RHS) is explicit … … 141 138 ! 142 139 INTEGER :: ji, ji2, jj, jj2, jn ! dummy loop indices 143 INTEGER :: jter, i_out, i_inn!140 INTEGER :: i_out, i_inn, i_inn_tot ! 144 141 INTEGER :: ji_min, jj_min ! 145 142 INTEGER :: nn_zebra_vp ! number of zebra steps 146 143 147 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp)148 144 ! 149 145 REAL(wp) :: zrhoco ! rho0 * rn_cio … … 152 148 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 153 149 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 154 REAL(wp) :: zd eltat, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars155 REAL(wp) :: zp_del tastar_f!150 REAL(wp) :: zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 151 REAL(wp) :: zp_delstar_f ! 156 152 REAL(wp) :: zu_cV, zv_cU ! 157 153 REAL(wp) :: zfac, zfac1, zfac2, zfac3 … … 164 160 REAL(wp), DIMENSION(jpi,jpj) :: zmassU_t, zmassV_t ! Mass per unit area divided by time step 165 161 ! 166 REAL(wp), DIMENSION(jpi,jpj) :: zdelta star_t !Delta* at T-points167 REAL(wp), DIMENSION(jpi,jpj) :: zten _i ! Tension168 REAL(wp), DIMENSION(jpi,jpj) :: zp_del tastar_t! P/delta* at T points162 REAL(wp), DIMENSION(jpi,jpj) :: zdeltat, zdelstar_t ! Delta & Delta* at T-points 163 REAL(wp), DIMENSION(jpi,jpj) :: ztens, zshear ! Tension, shear 164 REAL(wp), DIMENSION(jpi,jpj) :: zp_delstar_t ! P/delta* at T points 169 165 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zet ! Viscosity pre-factors at T points 170 166 REAL(wp), DIMENSION(jpi,jpj) :: zef ! Viscosity pre-factor at F point … … 194 190 REAL(wp), DIMENSION(jpi,jpj) :: zFU, zFU_prime, zBU_prime ! Rearranged linear system coefficients, U equation 195 191 REAL(wp), DIMENSION(jpi,jpj) :: zFV, zFV_prime, zBV_prime ! Rearranged linear system coefficients, V equation 196 REAL(wp), DIMENSION(jpi,jpj) :: zCU_prime, zCV_prime ! Rearranged linear system coefficients, V equation197 192 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) 198 193 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 199 194 ! 200 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 , zmsk15195 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 201 196 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! mask for lots of ice (1), little ice (0) 202 197 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence (1), no ice (0) … … 206 201 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 207 202 !! --- diags 208 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y203 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 209 204 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12, zs12f ! stress tensor components 210 205 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p … … 214 209 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp, zdiag_yatrp ! X/Y-component of area transport (m2/s, SIMIP) 215 210 216 217 CALL ctl_stop( 'STOP', 'icedyn_rhg_vp: stop because vp rheology is an ongoing work and should not be used' ) 211 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zvel_res ! Residual of the linear system at last iteration 212 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zvel_diff ! Absolute velocity difference @last outer iteration 213 218 214 219 215 !!---------------------------------------------------------------------------------------------------------------------- 220 !!$ ! DEBUG put all forcing terms to zero221 !!$ ! air-ice drag222 !!$ utau_ice(:,:) = 0._wp223 !!$ vtau_ice(:,:) = 0._wp224 !!$ ! coriolis225 !!$ ff_t(:,:) = 0._wp226 !!$ ! ice-ocean drag227 !!$ rn_cio = 0._wp228 !!$ ! ssh229 !!$ ! done line 330 !!! dont forget to act there230 !!$ ! END DEBUG231 216 232 217 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' … … 243 228 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 244 229 END_2D 245 IF( nn_rhg_chkcvg > 0 ) THEN246 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )247 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less248 END_2D249 ENDIF250 230 251 231 IF ( lp_zebra_vp ) THEN; nn_zebra_vp = 2 … … 267 247 IF( nn_rhg_chkcvg /= 0 ) THEN 268 248 269 ! ice area for global mean kinetic energy 270 zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) ) ! global ice area (km2)249 ! ice area for global mean kinetic energy (m2) 250 zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 271 251 272 252 ENDIF … … 279 259 280 260 zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 281 zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp; 282 zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp; 283 zrhsu(:,:) = 0._wp; zrhsv(:,:) = 0._wp 284 zf_rhsu(:,:) = 0._wp; zf_rhsv(:,:) = 0._wp 261 zrhsu (:,:) = 0._wp; zrhsv (:,:) = 0._wp; zf_rhsu(:,:) = 0._wp; zf_rhsv(:,:) = 0._wp 262 zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 263 zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp 285 264 286 265 !------------------------------------------------------------------------------! … … 292 271 CALL ice_strength ! strength at T points 293 272 294 !--------------------------- ---295 ! -- F-mask 296 !--------------------------- ---273 !--------------------------- 274 ! -- F-mask (code from EVP) 275 !--------------------------- 297 276 IF( kt == nit000 ) THEN 298 277 ! MartinV: … … 328 307 ! embedded sea ice: compute representative ice top surface 329 308 ! non-embedded sea ice: use ocean surface for slope calculation 330 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 331 !!$ zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 309 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 332 310 333 311 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 334 zm 1 = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) ) ! Ice/snow mass at U-V points335 zmf (ji,jj) = zm1 * ff_t(ji,jj) ! Coriolisat T points (m*f)312 zmt(ji,jj) = rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) ! Snow and ice mass at T-point 313 zmf(ji,jj) = zmt(ji,jj) * ff_t(ji,jj) ! Coriolis factor at T points (m*f) 336 314 END_2D 337 315 338 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 316 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 339 317 340 318 ! Ice fraction at U-V points 341 319 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) 342 320 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) 343 344 ! Ice/snowmass at U-V points345 zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ))346 zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ))347 zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) )321 322 ! Snow and ice mass at U-V points 323 zm1 = zmt(ji,jj) 324 zm2 = zmt(ji+1,jj) 325 zm3 = zmt(ji,jj+1) 348 326 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 349 327 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) … … 379 357 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 380 358 381 ! MV TEST DEBUG 382 !!$ IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji+1,jj) <= zmmin ) .AND. & 383 !!$ & ( at_i(ji,jj) <= zamin .OR. at_i(ji+1,jj) <= zamin ) ) THEN ; zmsk01x(ji,jj) = 0._wp 384 !!$ ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 385 !!$ 386 !!$ IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji,jj+1) <= zmmin ) .AND. & 387 !!$ & ( at_i(ji,jj) <= zamin .OR. at_i(ji,jj+1) <= zamin ) ) THEN ; zmsk01y(ji,jj) = 0._wp 388 !!$ ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 389 ! END MV TEST DEBUG 390 391 END_2D 392 393 !!$ CALL iom_put( 'zmsk00x' , zmsk00x ) ! MV DEBUG 394 !!$ CALL iom_put( 'zmsk00y' , zmsk00y ) ! MV DEBUG 395 !!$ CALL iom_put( 'zmsk01x' , zmsk01x ) ! MV DEBUG 396 !!$ CALL iom_put( 'zmsk01y' , zmsk01y ) ! MV DEBUG 397 !!$ CALL iom_put( 'ztaux_ai' , ztaux_ai ) ! MV DEBUG 398 !!$ CALL iom_put( 'ztauy_ai' , ztauy_ai ) ! MV DEBUG 399 !!$ CALL iom_put( 'zspgU' , zspgU ) ! MV DEBUG 400 !!$ CALL iom_put( 'zspgV' , zspgV ) ! MV DEBUG 401 !!$ 359 END_2D 360 402 361 !------------------------------------------------------------------------------! 403 362 ! … … 409 368 zv_c(:,:) = v_ice(:,:) 410 369 411 jter= 0370 i_inn_tot = 0 412 371 413 372 DO i_out = 1, nn_vp_nout 414 373 415 IF( lwp ) WRITE(numout,*) ' outer loop i_out : ', i_out416 417 374 ! 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} )375 ! u_it = 0.5 * ( u_{it-1} + u_{it-2} ) 419 376 ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009) 420 377 zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) ) … … 428 385 ! In the outer loop, one needs to update all RHS terms 429 386 ! with explicit velocity dependencies (viscosities, coriolis, ocean stress) 430 ! as a function of uc 431 ! 387 ! as a function of "current" velocities (uc, vc) 432 388 433 389 !------------------------------------------ … … 436 392 437 393 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 438 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 439 394 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! 1->jpi-1 395 396 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 440 397 ! shear at F points 441 398 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) & … … 445 402 END_2D 446 403 447 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 448 !!$ CALL iom_put( 'zds' , zds ) ! MV DEBUG 449 450 IF( lwp ) WRITE(numout,*) ' outer loop 1a i_out : ', i_out 451 452 !DO jj = 2, jpj - 1 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 453 ! DO ji = 2, jpi - 1 ! 454 455 ! MV DEBUG 456 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 457 DO ji = 2, jpi ! 458 ! END MV DEBUG 459 460 ! shear**2 at T points (doc eq. A16) 461 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 462 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 463 & ) * 0.25_wp * r1_e1e2t(ji,jj) 404 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! necessary, zds2 uses jpi/jpj values for zds 405 406 DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! 2 -> jpj; 2,jpi !!! CHECK !!! 407 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 408 409 ! shear**2 at T points (doc eq. A16) 410 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 411 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 412 & ) * 0.25_wp * r1_e1e2t(ji,jj) 464 413 465 466 467 468 469 414 ! divergence at T points 415 zdiv = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) & 416 & + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) & 417 & ) * r1_e1e2t(ji,jj) 418 zdiv2 = zdiv * zdiv 470 419 471 472 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) &473 &- ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &474 &) * r1_e1e2t(ji,jj)475 zdt2= zdt * zdt420 ! tension at T points 421 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) & 422 & - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 423 & ) * r1_e1e2t(ji,jj) 424 zdt2 = zdt * zdt 476 425 477 478 zdeltat= SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )426 ! delta at T points 427 zdeltat(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 479 428 480 481 zdeltastar_t(ji,jj) = zdeltat + rn_creepl429 ! delta* at T points (following Lemieux and Dupont, GMD 2020) 430 zdelstar_t(ji,jj) = zdeltat(ji,jj) + rn_creepl ! OPT zdelstar_t can be totally removed and put into next line directly. Could change results 482 431 483 ! P/deltaat T-points484 zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj)432 ! P/delta* at T-points 433 zp_delstar_t(ji,jj) = strength(ji,jj) / zdelstar_t(ji,jj) 485 434 486 487 zzt(ji,jj) = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)488 zet(ji,jj)= zzt(ji,jj) * z1_ecc2435 ! Temporary zzt and zet factors at T-points 436 zzt(ji,jj) = zp_delstar_t(ji,jj) * r1_e1e2t(ji,jj) 437 zet(ji,jj) = zzt(ji,jj) * z1_ecc2 489 438 490 END DO 491 END DO 492 493 CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 494 495 !!$ CALL iom_put( 'zzt' , zzt ) ! MV DEBUG 496 !!$ CALL iom_put( 'zet' , zet ) ! MV DEBUG 497 !!$ CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 498 499 IF( lwp ) WRITE(numout,*) ' outer loop 1b i_out : ', i_out 500 501 DO jj = 1, jpj - 1 502 DO ji = 1, jpi - 1 503 439 END_2D 440 441 CALL lbc_lnk( 'icedyn_rhg_vp', zp_delstar_t , 'T', 1. ) ! necessary, used for ji = 1 and jj = 1 442 443 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )! 1-> jpj-1; 1->jpi-1 444 504 445 ! P/delta* at F points 505 zp_del tastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) )446 zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) ) 506 447 507 448 ! Temporary zef factor at F-point 508 zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * fimask(ji,jj) 509 510 END DO 511 END DO 512 513 CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 514 !!$ CALL iom_put( 'zef' , zef ) ! MV DEBUG 515 IF( lwp ) WRITE(numout,*) ' outer loop 1c i_out : ', i_out 516 449 zef(ji,jj) = zp_delstar_f * r1_e1e2f(ji,jj) * z1_ecc2 * fimask(ji,jj) * 0.5_wp 450 451 END_2D 452 517 453 !--------------------------------------------------- 518 454 ! -- Ocean-ice drag and Coriolis RHS contributions 519 455 !--------------------------------------------------- 520 456 521 DO jj = 2, jpj - 1 522 DO ji = 2, jpi - 1 457 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 458 459 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 460 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) 461 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) 523 462 524 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 525 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) 526 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) 463 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 464 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) ) & 465 & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 466 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) ) & 467 & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 468 469 !--- Ocean-ice drag contributions to RHS 470 ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj) 471 ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj) 527 472 528 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 529 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) ) & 530 & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 531 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) ) & 532 & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 533 534 !--- Ocean-ice drag contributions to RHS 535 ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj) 536 ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj) 537 538 ! --- U-component of Coriolis Force (energy conserving formulation) 539 ! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done 540 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 541 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * zv_c(ji ,jj) + e1v(ji ,jj-1) * zv_c(ji ,jj-1) ) & 542 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) ) 473 !--- U-component of Coriolis Force (energy conserving formulation) 474 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 475 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * zv_c(ji ,jj) + e1v(ji ,jj-1) * zv_c(ji ,jj-1) ) & 476 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) ) 543 477 544 ! --- V-component of Coriolis Force (energy conserving formulation) 545 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 546 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * zu_c(ji,jj ) + e2u(ji-1,jj ) * zu_c(ji-1,jj ) ) & 547 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 548 549 END DO 550 END DO 551 552 IF( lwp ) WRITE(numout,*) ' outer loop 1d i_out : ', i_out 553 554 CALL lbc_lnk( 'icedyn_rhg_vp', zCwU , 'U', -1., zCwV, 'V', -1. ) 555 CALL lbc_lnk( 'icedyn_rhg_vp', zCorU, 'U', -1., zCorV, 'V', -1. ) 556 557 !!$ CALL iom_put( 'zCwU' , zCwU ) ! MV DEBUG 558 !!$ CALL iom_put( 'zCwV' , zCwV ) ! MV DEBUG 559 !!$ CALL iom_put( 'zCorU' , zCorU ) ! MV DEBUG 560 !!$ CALL iom_put( 'zCorV' , zCorV ) ! MV DEBUG 561 562 IF( lwp ) WRITE(numout,*) ' outer loop 1f i_out : ', i_out 563 564 ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system, 565 ! so there is no need for lbclnk on drag and coriolis 566 478 !--- V-component of Coriolis Force (energy conserving formulation) 479 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 480 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * zu_c(ji,jj ) + e2u(ji-1,jj ) * zu_c(ji-1,jj ) ) & 481 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 482 483 END_2D 484 567 485 !------------------------------------- 568 486 ! -- Internal stress RHS contribution 569 487 !------------------------------------- 570 488 571 ! --- Stress contributions at T-points 572 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 573 DO ji = 2, jpi ! 489 ! --- Stress contributions at T-points 490 DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! 2 -> jpj; 2,jpi !!! CHECK !!! 491 492 ! loop to jpi,jpj to avoid making a communication for zs1 & zs2 574 493 575 ! sig1 contribution to RHS of U-equation at T-points 576 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 ) 494 ! sig1 contribution to RHS of U-equation at T-points 495 zs1_rhsu(ji,jj) = zzt(ji,jj) * ( e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) ) & 496 & - zp_delstar_t(ji,jj) * zdeltat(ji,jj) 577 497 578 ! sig2 contribution to RHS of U-equation at T-points 579 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) 580 581 ! sig1 contribution to RHS of V-equation at T-points 582 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 ) 583 584 ! sig2 contribution to RHS of V-equation at T-points 585 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) 586 587 END DO 588 END DO 589 590 !!$ CALL iom_put( 'zs1_rhsu' , zs1_rhsu ) ! MV DEBUG 591 !!$ CALL iom_put( 'zs2_rhsu' , zs2_rhsu ) ! MV DEBUG 592 !!$ CALL iom_put( 'zs1_rhsv' , zs1_rhsv ) ! MV DEBUG 593 !!$ CALL iom_put( 'zs2_rhsv' , zs2_rhsv ) ! MV DEBUG 594 595 ! a priori, no lbclnk, because rhsu is only used in the inner domain 596 597 ! --- Stress contributions at f-points 498 ! sig2 contribution to RHS of U-equation at T-points 499 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) 500 501 ! sig1 contribution to RHS of V-equation at T-points 502 zs1_rhsv(ji,jj) = zzt(ji,jj) * ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) ) & 503 & - zp_delstar_t(ji,jj) * zdeltat(ji,jj) 504 505 ! sig2 contribution to RHS of V-equation at T-points 506 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) 507 508 END_2D 509 510 ! --- Stress contributions at F-points 598 511 ! MV NOTE: I applied fimask on zds, by mimetism on EVP, but without deep understanding of what I was doing 599 512 ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 600 513 601 IF( lwp ) WRITE(numout,*) ' outer loop 2 i_out : ', i_out 602 603 DO jj = 1, jpj - 1 604 DO ji = 1, jpi - 1 605 606 ! sig12 contribution to RHS of U equation at F-points 607 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) * fimask(ji,jj) 608 609 ! sig12 contribution to RHS of V equation at F-points 610 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) * fimask(ji,jj) 611 612 END DO 613 END DO 614 615 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 616 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) 617 618 !!$ CALL iom_put( 'zs12_rhsu' , zs12_rhsu ) ! MV DEBUG 619 !!$ CALL iom_put( 'zs12_rhsv' , zs12_rhsv ) ! MV DEBUG 620 621 ! a priori, no lbclnk, because rhsu are only used in the inner domain 622 514 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! 1->jpi-1 515 516 ! sig12 contribution to RHS of U equation at F-points 517 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) * fimask(ji,jj) 518 519 ! sig12 contribution to RHS of V equation at F-points 520 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) * fimask(ji,jj) 521 522 END_2D 523 623 524 ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) 624 525 ! OPT: merge with next loop and use intermediate scalars for zf_rhsu 625 626 DO jj = 2, jpj - 1 627 DO ji = 2, jpi - 1 526 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 527 628 528 ! --- U component of internal force contribution to RHS at U points 629 529 zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & … … 635 535 zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * & 636 536 & ( e1v(ji,jj) * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) ) & 637 & + 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) ) &537 & - 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) ) & 638 538 & + 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) ) ) 639 640 END DO 641 END DO 642 643 !!$ CALL iom_put( 'zf_rhsu' , zf_rhsu ) ! MV DEBUG 644 !!$ CALL iom_put( 'zf_rhsv' , zf_rhsv ) ! MV DEBUG 539 540 END_2D 645 541 646 542 !--------------------------- … … 649 545 ! 650 546 ! OPT: could use intermediate scalars to reduce memory access 651 DO jj = 2, jpj - 1 652 DO ji = 2, jpi - 1 653 654 ! still miss ice ocean stress and acceleration contribution 655 zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj) 656 zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsu(ji,jj) 657 658 END DO 659 END DO 660 661 CALL lbc_lnk( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V', -1.) 662 CALL lbc_lnk( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V', -1.) 663 CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V', -1.) 664 665 !!$ CALL iom_put( 'zmU_t' , zmU_t ) ! MV DEBUG 666 !!$ CALL iom_put( 'zmV_t' , zmV_t ) ! MV DEBUG 667 !!$ CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu ) ! MV DEBUG 668 !!$ CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv ) ! MV DEBUG 669 !!$ CALL iom_put( 'zrhsu' , zrhsu ) ! MV DEBUG 670 !!$ CALL iom_put( 'zrhsv' , zrhsv ) ! MV DEBUG 671 !!$ 672 ! inner domain calculations -> no lbclnk 673 674 IF( lwp ) WRITE(numout,*) ' outer loop 4 i_out : ', i_out 675 547 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 548 549 zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj) 550 zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsv(ji,jj) 551 552 END_2D 553 676 554 !---------------------------------------------------------------------------------------! 677 555 ! … … 691 569 ! only zzt and zet are iteration-dependent, other only depend on scale factors 692 570 693 DO ji = 2, jpi - 1 ! internal domain do loop 694 DO jj = 2, jpj - 1 695 696 !------------------------------------- 697 ! -- Internal forces LHS contribution 698 !------------------------------------- 699 ! 700 ! --- U-component 701 ! 702 ! "T" factors (intermediate results) 703 ! 704 zfac = 0.5_wp * r1_e1e2u(ji,jj) 705 zfac1 = zfac * e2u(ji,jj) 706 zfac2 = zfac * r1_e2u(ji,jj) 707 zfac3 = 2._wp * zfac * r1_e1u(ji,jj) 708 709 zt12U = - zfac1 * zzt(ji+1,jj) 710 zt11U = zfac1 * zzt(ji,jj) 711 712 zt22U = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 713 zt21U = zfac2 * zet(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) 714 715 zt122U = - zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) 716 zt121U = zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 571 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 572 573 !------------------------------------- 574 ! -- Internal forces LHS contribution 575 !------------------------------------- 576 ! 577 ! --- U-component 578 ! 579 ! "T" factors (intermediate results) 580 ! 581 zfac = 0.5_wp * r1_e1e2u(ji,jj) 582 zfac1 = zfac * e2u(ji,jj) 583 zfac2 = zfac * r1_e2u(ji,jj) 584 zfac3 = 2._wp * zfac * r1_e1u(ji,jj) 585 586 zt11U = zfac1 * zzt(ji,jj) 587 zt12U = zfac1 * zzt(ji+1,jj) 588 589 zt21U = zfac2 * zet(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) 590 zt22U = zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 591 592 zt121U = zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 593 zt122U = zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) 717 594 718 719 720 721 zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U* r1_e2u(ji-1,jj)722 zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj)723 zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U* r1_e2u(ji+1,jj)724 725 zDU(ji,jj) =zt121U * r1_e1u(ji,jj-1)726 zEU(ji,jj) =zt122U * r1_e1u(ji,jj+1)595 ! 596 ! Linear system coefficients 597 ! 598 zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 599 zBU(ji,jj) = ( zt11U + zt12U ) * e2u(ji,jj) + ( zt21U + zt22U ) * r1_e2u(ji,jj) + ( zt121U + zt122U ) * r1_e1u(ji,jj) 600 zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 601 602 zDU(ji,jj) = zt121U * r1_e1u(ji,jj-1) 603 zEU(ji,jj) = zt122U * r1_e1u(ji,jj+1) 727 604 728 729 730 731 732 733 734 zfac1 = zfac * e2v(ji,jj)735 736 737 738 zt12V = - zfac1 * zzt(ji,jj+1)739 zt11V = zfac1 * zzt(ji,jj)740 741 zt22V = zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1)742 zt21V = - zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj)743 744 zt122V = zfac3 * zef(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj)745 zt121V = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)746 747 748 749 750 zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V* r1_e1v(ji,jj-1)751 zBV(ji,jj) = ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) -( zt122V + zt121V ) * r1_e2v(ji,jj)752 zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V* r1_e1v(ji,jj+1)753 754 zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here755 zEV(ji,jj) = -zt122V * r1_e2v(ji+1,jj)605 ! 606 ! --- V-component 607 ! 608 ! "T" factors (intermediate results) 609 ! 610 zfac = 0.5_wp * r1_e1e2v(ji,jj) 611 zfac1 = zfac * e1v(ji,jj) 612 zfac2 = zfac * r1_e1v(ji,jj) 613 zfac3 = 2._wp * zfac * r1_e2v(ji,jj) 614 615 zt11V = zfac1 * zzt(ji,jj) 616 zt12V = zfac1 * zzt(ji,jj+1) 617 618 zt21V = zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) 619 zt22V = zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 620 621 zt121V = zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 622 zt122V = zfac3 * zef(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) 623 624 ! 625 ! Linear system coefficients 626 ! 627 zAV(ji,jj) = - zt11V * e1v(ji,jj-1) - zt21V * r1_e1v(ji,jj-1) 628 zBV(ji,jj) = ( zt11V + zt12V ) * e1v(ji,jj) + ( zt21V + zt22V ) * r1_e1v(ji,jj) + ( zt122V + zt121V ) * r1_e2v(ji,jj) 629 zCV(ji,jj) = - zt12V * e1v(ji,jj+1) - zt22V * r1_e1v(ji,jj+1) 630 631 zDV(ji,jj) = zt121V * r1_e2v(ji-1,jj) 632 zEV(ji,jj) = zt122V * r1_e2v(ji+1,jj) 756 633 757 !----------------------------------------------------- 758 ! -- Ocean-ice drag and acceleration LHS contribution 759 !----------------------------------------------------- 760 zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 761 zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 762 763 END DO 764 END DO 765 766 CALL lbc_lnk( 'icedyn_rhg_vp', zAU , 'U', 1., zAV , 'V', 1. ) 767 CALL lbc_lnk( 'icedyn_rhg_vp', zBU , 'U', 1., zBV , 'V', 1. ) 768 CALL lbc_lnk( 'icedyn_rhg_vp', zCU , 'U', 1., zCV , 'V', 1. ) 769 CALL lbc_lnk( 'icedyn_rhg_vp', zDU , 'U', 1., zDV , 'V', 1. ) 770 CALL lbc_lnk( 'icedyn_rhg_vp', zEU , 'U', 1., zEV , 'V', 1. ) 771 772 !!$ CALL iom_put( 'zAU' , zAU ) ! MV DEBUG 773 !!$ CALL iom_put( 'zBU' , zBU ) ! MV DEBUG 774 !!$ CALL iom_put( 'zCU' , zCU ) ! MV DEBUG 775 !!$ CALL iom_put( 'zDU' , zDU ) ! MV DEBUG 776 !!$ CALL iom_put( 'zEU' , zEU ) ! MV DEBUG 777 !!$ CALL iom_put( 'zAV' , zAV ) ! MV DEBUG 778 !!$ CALL iom_put( 'zBV' , zBV ) ! MV DEBUG 779 !!$ CALL iom_put( 'zCV' , zCV ) ! MV DEBUG 780 !!$ CALL iom_put( 'zDV' , zDV ) ! MV DEBUG 781 !!$ CALL iom_put( 'zEV' , zEV ) ! MV DEBUG 782 634 !----------------------------------------------------- 635 ! -- Ocean-ice drag and acceleration LHS contribution 636 !----------------------------------------------------- 637 zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 638 zBV(ji,jj) = zBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 639 640 END_2D 641 783 642 !------------------------------------------------------------------------------! 784 643 ! … … 793 652 DO i_inn = 1, nn_vp_ninn ! inner loop iterations 794 653 795 IF( lwp ) WRITE(numout,*) ' inner loop 1 i_inn : ', i_inn796 797 654 !--- mitgcm computes initial value of residual here... 798 655 799 jter = jter + 1 800 ! l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 801 802 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 803 zv_b(:,:) = v_ice(:,:) 804 805 ! zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 806 ! zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp 656 i_inn_tot = i_inn_tot + 1 657 ! l_full_nf_update = i_inn_tot == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 658 659 zu_b(:,:) = u_ice(:,:) ! velocity at previous inner-iterate 660 zv_b(:,:) = v_ice(:,:) 807 661 808 662 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN … … 817 671 ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 818 672 819 zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; zCU_prime(:,:) = 0._wp673 zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; 820 674 821 675 DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) … … 826 680 ELSE ; jj_min = 3 827 681 ENDIF 828 829 IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min830 682 831 683 DO jj = jj_min, jpj - 1, nn_zebra_vp … … 835 687 !------------------------ 836 688 DO ji = 2, jpi - 1 837 838 ! boundary condition substitution 689 ! note: these are key lines linking information between processors 690 ! u_ice/v_ice need to be lbc_linked 691 692 ! sub-domain boundary condition substitution 839 693 ! see Zhang and Hibler, 1997, Appendix B 840 694 zAA3 = 0._wp … … 852 706 END DO 853 707 854 CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. )855 856 708 !--------------- 857 709 ! Forward sweep … … 859 711 DO jj = jj_min, jpj - 1, nn_zebra_vp 860 712 713 zBU_prime(2,jj) = zBU(2,jj) 714 zFU_prime(2,jj) = zFU(2,jj) 715 861 716 DO ji = 3, jpi - 1 862 717 … … 869 724 870 725 END DO 871 872 CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime, 'U', 1., zBU_prime, 'U', 1. ) 873 726 874 727 !----------------------------- 875 728 ! Backward sweep & relaxation … … 879 732 880 733 ! --- Backward sweep 734 881 735 ! last row 882 736 zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 883 737 u_ice(jpi-1,jj) = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 884 738 & * umask(jpi-1,jj,1) 885 DO ji = jpi-2 , 2, -1 ! all other rows ! ---> original backward loop 739 740 DO ji = jpi - 2 , 2, -1 ! all other rows ! ---> original backward loop 886 741 zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 887 742 u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) & … … 889 744 END DO 890 745 891 !--- Relaxation 892 ! and velocity masking for little-ice and no-ice cases 746 !--- Relaxation and masking (for low-ice/no-ice cases) 893 747 DO ji = 2, jpi - 1 894 748 … … 902 756 903 757 END DO ! jj 758 759 CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1. ) 904 760 905 761 END DO ! zebra loop … … 917 773 !!! ZH97 explain it is critical for convergence speed 918 774 919 zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp775 zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; 920 776 921 777 DO jn = 1, nn_zebra_vp ! "zebra" loop … … 925 781 ENDIF 926 782 927 IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min928 929 783 DO ji = ji_min, jpi - 1, nn_zebra_vp 930 784 … … 934 788 DO jj = 2, jpj - 1 935 789 936 ! boundary condition substitution (check it is correctly applied !!!)790 ! subdomain boundary condition substitution (check it is correctly applied !!!) 937 791 ! see Zhang and Hibler, 1997, Appendix B 938 792 zAA3 = 0._wp … … 950 804 END DO 951 805 952 CALL lbc_lnk( 'icedyn_rhg_vp', zFV, 'V', 1.)953 954 806 !--------------- 955 807 ! Forward sweep … … 957 809 DO ji = ji_min, jpi - 1, nn_zebra_vp 958 810 959 DO jj = 3, jpj - 1 811 zBV_prime(ji,2) = zBV(ji,2) 812 zFV_prime(ji,2) = zFV(ji,2) 813 814 DO jj = 3, jpj - 1 960 815 961 816 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) … … 968 823 END DO 969 824 970 CALL lbc_lnk( 'icedyn_rhg_vp', zFV_prime, 'V', 1., zBV_prime, 'V', 1. )971 972 825 !----------------------------- 973 826 ! Backward sweep & relaxation … … 988 841 END DO 989 842 990 ! --- Relaxation & masking (should it be now or later)843 ! --- Relaxation & masking 991 844 DO jj = 2, jpj - 1 992 845 … … 1000 853 1001 854 END DO ! ji 855 856 CALL lbc_lnk( 'icedyn_rhg_vp', v_ice, 'V', -1. ) 1002 857 1003 858 END DO ! zebra loop … … 1005 860 ENDIF ! ll_v_iterate 1006 861 1007 CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )862 ! I suspect the communication should go into the zebra loop if we want reproducibility 1008 863 1009 864 !-------------------------------------------------------------------------------------- … … 1016 871 ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 1017 872 ! then we must optimize the use of the mpp_max, which is prohibitive 1018 zuerr_max = 0._wp873 zuerr_max = 0._wp 1019 874 1020 875 IF ( ll_u_iterate .AND. MOD ( i_inn, nn_vp_chkcvg ) == 0 ) THEN … … 1022 877 ! - Maximum U-velocity difference 1023 878 zuerr(:,:) = 0._wp 1024 DO jj = 2, jpj - 1 1025 DO ji = 2, jpi - 1 1026 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 1027 END DO 1028 END DO 879 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 880 881 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 882 883 END_2D 884 1029 885 zuerr_max = MAXVAL( zuerr ) 1030 886 CALL mpp_max( 'icedyn_rhg_evp', zuerr_max ) ! max over the global domain - damned! 1031 1032 ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)887 888 ! - Stop if max error is too large ("safeguard against bad forcing" of original Zhang routine) 1033 889 IF ( i_inn > 1 .AND. zuerr_max > zuerr_max_vp ) THEN 1034 890 IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " … … 1053 909 ! - Maximum V-velocity difference 1054 910 zverr(:,:) = 0._wp 1055 DO jj = 2, jpj -11056 DO ji = 2, jpi - 1911 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 912 1057 913 zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1) 1058 END DO1059 END DO1060 914 915 END_2D 916 1061 917 zverr_max = MAXVAL( zverr ) 1062 918 CALL mpp_max( 'icedyn_rhg_evp', zverr_max ) ! max over the global domain - damned! … … 1083 939 ! 1084 940 !--------------------------------------------------------------------------------------- 1085 1086 IF( nn_rhg_chkcvg/=0 .AND. MOD ( i_inn - 1, nn_vp_chkcvg ) == 0 ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 1087 & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 1088 1089 IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 941 IF( nn_rhg_chkcvg/=0 .AND. MOD ( i_inn - 1, nn_vp_chkcvg ) == 0 ) THEN 942 943 CALL rhg_cvg_vp( kt, i_out, i_inn, i_inn_tot, nn_vp_nout, nn_vp_ninn, nn_nvp, & 944 & u_ice, v_ice, zu_b, zv_b, zu_c, zv_c, & 945 & zmt, za_iU, za_iV, zuerr_max, zverr_max, zglob_area, & 946 & zrhsu, zAU, zBU, zCU, zDU, zEU, zFU, & 947 & zrhsv, zAV, zBV, zCV, zDV, zEV, zFV, & 948 zvel_res, zvel_diff ) 949 950 ENDIF 1090 951 1091 952 END DO ! i_inn, end of inner loop … … 1093 954 END DO ! End of outer loop (i_out) ============================================================================================= 1094 955 1095 IF ( lwp ) WRITE(numout,*) ' We are out of outer loop ' 1096 1097 CALL lbc_lnk( 'icedyn_rhg_vp', zFU , 'U', 1., zFV , 'V', 1. ) 1098 CALL lbc_lnk( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zBV_prime , 'V', 1. ) 1099 CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime , 'U', 1., zFV_prime , 'V', 1. ) 1100 CALL lbc_lnk( 'icedyn_rhg_vp', zCU_prime , 'U', 1., zCV_prime , 'V', 1. ) 1101 1102 !!$ CALL iom_put( 'zFU' , zFU ) ! MV DEBUG 1103 !!$ CALL iom_put( 'zBU_prime' , zBU_prime ) ! MV DEBUG 1104 !!$ CALL iom_put( 'zCU_prime' , zCU_prime ) ! MV DEBUG 1105 !!$ CALL iom_put( 'zFU_prime' , zFU_prime ) ! MV DEBUG 1106 !!$ 1107 !!$ CALL iom_put( 'zFV' , zFV ) ! MV DEBUG 1108 !!$ CALL iom_put( 'zBV_prime' , zBV_prime ) ! MV DEBUG 1109 !!$ CALL iom_put( 'zCV_prime' , zCV_prime ) ! MV DEBUG 1110 !!$ CALL iom_put( 'zFV_prime' , zFV_prime ) ! MV DEBUG 1111 1112 CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1113 1114 !!$ IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg ' 1115 !!$ IF( iom_use('uice_dbg' ) ) CALL iom_put( 'uice_dbg' , u_ice ) ! ice velocity u after solver 1116 !!$ IF( iom_use('vice_dbg' ) ) CALL iom_put( 'vice_dbg' , v_ice ) ! ice velocity v after solver 1117 956 IF( nn_rhg_chkcvg/=0 ) THEN 957 958 IF( iom_use('velo_res') ) CALL iom_put( 'velo_res', zvel_res ) ! linear system residual @last inner&outer iteration 959 IF( iom_use('velo_ero') ) CALL iom_put( 'velo_ero', zvel_diff ) ! abs velocity difference @last outer iteration 960 IF( iom_use('uice_eri') ) CALL iom_put( 'uice_eri', zuerr ) ! abs velocity difference @last inner iteration 961 IF( iom_use('vice_eri') ) CALL iom_put( 'vice_eri', zverr ) ! abs velocity difference @last inner iteration 962 963 DEALLOCATE( zvel_res , zvel_diff ) 964 965 ENDIF ! nn_rhg_chkcvg 966 1118 967 !------------------------------------------------------------------------------! 1119 968 ! 1120 ! --- Convergence diagnostics969 ! --- Recompute delta, shear and div (inputs for mechanical redistribution) 1121 970 ! 1122 971 !------------------------------------------------------------------------------! 1123 1124 IF( nn_rhg_chkcvg /= 0 ) THEN1125 1126 IF( iom_use('uice_cvg') ) THEN1127 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & ! ice velocity difference at last iteration1128 & ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )1129 ENDIF1130 1131 ENDIF1132 1133 !!$ ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead1134 !!$ DO jj = 2, jpj - 11135 !!$ DO ji = 2, jpi - 11136 !!$1137 !!$ u_ice(ji,jj) = zmsk00x(ji,jj) &1138 !!$ & * ( zmsk01x(ji,jj) * u_oce(ji,jj) * 0.01_wp &1139 !!$ + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp )1140 !!$1141 !!$ v_ice(ji,jj) = zmsk00y(ji,jj) &1142 !!$ & * ( zmsk01y(ji,jj) * v_oce(ji,jj) * 0.01_wp &1143 !!$ + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp )1144 !!$1145 !!$ END DO1146 !!$ END DO1147 !!$1148 !!$ CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )1149 !!$1150 !!$ IF ( lwp ) WRITE(numout,*) ' Velocity replaced '1151 1152 ! END DEBUG1153 1154 !------------------------------------------------------------------------------!1155 !1156 ! --- Recompute delta, shear and div (inputs for mechanical redistribution)1157 !1158 !------------------------------------------------------------------------------!1159 972 ! 1160 973 ! MV OPT: subroutinize ? 1161 1162 DO jj = 1, jpj - 1 1163 DO ji = 1, jpi - 1 974 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) ! 1->jpj-1; 1->jpi-1 1164 975 1165 976 ! shear at F points … … 1168 979 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 1169 980 1170 END DO 1171 END DO 1172 1173 DO jj = 2, jpj - 1 1174 DO ji = 2, jpi - 1 ! 981 END_2D 982 983 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1175 984 1176 985 ! tension**2 at T points … … 1180 989 zdt2 = zdt * zdt 1181 990 1182 zten _i(ji,jj)= zdt991 ztens(ji,jj) = zdt 1183 992 1184 993 ! shear**2 at T points (doc eq. A16) … … 1187 996 & ) * 0.25_wp * r1_e1e2t(ji,jj) 1188 997 1189 ! shear at T points 1190 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 998 ! maximum shear rate at T points (includees tension, output only) 999 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) ! i think this is maximum shear rate and not actual shear --- i'm not totally sure here 1000 1001 ! shear at T-points 1002 zshear(ji,jj) = SQRT( zds2 ) 1191 1003 1192 1004 ! divergence at T points … … 1194 1006 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 1195 1007 & ) * r1_e1e2t(ji,jj) 1008 1009 zdiv2 = pdivu_i(ji,jj) * pdivu_i(ji,jj) 1196 1010 1197 1011 ! delta at T points 1198 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )1012 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 1199 1013 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 1200 1201 !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 1202 pdelta_i(ji,jj) = zdelta + rn_creepl 1203 1204 END DO 1205 END DO 1206 1207 IF ( lwp ) WRITE(numout,*) ' Deformation recalculated ' 1014 1015 pdelta_i(ji,jj) = zdelta + rn_creepl ! * rswitch 1016 1017 END_2D 1208 1018 1209 1019 CALL lbc_lnk( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) … … 1222 1032 ! 1223 1033 ! ---- Sea ice stresses at T-points 1224 IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1225 1226 DO jj = 2, jpj - 1 1227 DO ji = 2, jpi - 1 1228 zp_deltastar_t(ji,jj) = strength(ji,jj) / pdelta_i(ji,jj) 1229 zfac = zp_deltastar_t(ji,jj) 1034 IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. & 1035 & iom_use('intstrx') .OR. iom_use('intstry') .OR. & 1036 & iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 1037 1038 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1039 1040 zp_delstar_t(ji,jj) = strength(ji,jj) / pdelta_i(ji,jj) 1041 zfac = zp_delstar_t(ji,jj) 1230 1042 zs1(ji,jj) = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 1231 zs2(ji,jj) = zfac * z1_ecc2 * zten _i(ji,jj)1232 zs12(ji,jj) = zfac * z1_ecc2 * pshear_i(ji,jj)1233 END DO 1234 END DO1043 zs2(ji,jj) = zfac * z1_ecc2 * ztens(ji,jj) 1044 zs12(ji,jj) = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bug 12 nov 1045 1046 END_2D 1235 1047 1236 1048 CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. ) … … 1241 1053 IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1242 1054 1243 DO jj = 1, jpj - 1 1244 DO ji = 1, jpi - 1 1055 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) ! 1->jpj-1; 1->jpi-1 1245 1056 1246 1057 ! P/delta* at F points 1247 zp_del tastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) )1058 zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) ) 1248 1059 1249 1060 ! s12 at F-points 1250 zs12f(ji,jj) = zp_del tastar_f * z1_ecc2 * zds(ji,jj)1061 zs12f(ji,jj) = zp_delstar_f * z1_ecc2 * zds(ji,jj) 1251 1062 1252 END DO 1253 END DO 1063 END_2D 1254 1064 1255 1065 CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. ) 1256 1066 1257 1067 ENDIF 1258 1259 IF ( lwp ) WRITE(numout,*) ' zs12f recalculated '1260 1068 1261 1069 ! … … 1271 1079 1272 1080 !--- Recalculate oceanic stress at last inner iteration 1273 DO jj = 2, jpj - 1 1274 DO ji = 2, jpi - 1 1081 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1275 1082 1276 1083 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) … … 1288 1095 ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) ) 1289 1096 1290 END DO 1291 END DO 1097 END_2D 1292 1098 1293 1099 ! 1294 1100 CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, & 1295 ! & 1101 ! & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 1296 1102 ! 1297 1103 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) … … 1308 1114 ! --- Divergence, shear and strength --- ! 1309 1115 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 1310 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear1116 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! maximum shear rate 1311 1117 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 1312 1118 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 1313 1119 1314 IF ( lwp ) WRITE(numout,*) 'Some terms recalculated '1315 1316 1120 ! --- Stress tensor invariants (SIMIP diags) --- ! 1317 1121 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN … … 1325 1129 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 1326 1130 ! 1327 DO jj = 2, jpj - 1 1328 DO ji = 2, jpi - 1 1131 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1329 1132 ! Stress invariants 1330 zsig_I(ji,jj) = zs1(ji,jj) * 0.5_wp ! 1st invariant, aka average normal stress aka negative pressure 1331 zsig_II(ji,jj) = SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) ) ! 2nd invariant, aka maximum shear stress 1332 END DO 1333 END DO 1133 zsig_I(ji,jj) = zs1(ji,jj) * 0.5_wp ! 1st invariant, aka average normal stress aka negative pressure 1134 zsig_II(ji,jj) = 0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) ) ! 2nd invariant, aka maximum shear stress 1135 END_2D 1334 1136 1335 1137 CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.) … … 1341 1143 1342 1144 ENDIF 1343 1344 IF ( lwp ) WRITE(numout,*) 'SIMIP work done'1345 1145 1346 1146 ! --- Normalized stress tensor principal components --- ! … … 1355 1155 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 1356 1156 ! 1357 DO jj = 2, jpj -11358 DO ji = 2, jpi - 11157 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1158 1359 1159 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 1360 1160 ! and **deformations** at current iterates 1361 1161 ! following Lemieux & Dupont (2020) 1362 zfac = zp_deltastar_t(ji,jj) 1363 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 1364 !!$ zsig1 = 0._wp !!! FUCKING DEBUG TEST !!! 1365 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 1366 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 1162 zfac = zp_delstar_t(ji,jj) 1163 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) ) 1164 zsig2 = zfac * z1_ecc2 * ztens(ji,jj) 1165 zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bugfix 12 Nov 1367 1166 1368 1167 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 1369 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st invariant1370 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp +zsig12 ) ! 2nd invariant1168 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st invariant 1169 zsig_II(ji,jj) = 0.5_wp * SQRT ( zsig2 * zsig2 + 4. *zsig12 * zsig12 ) ! 2nd invariant 1371 1170 1372 1171 ! Normalized principal stresses (used to display the ellipse) … … 1374 1173 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 1375 1174 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 1376 END DO 1377 END DO 1378 IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done' 1175 1176 END_2D 1379 1177 ! 1380 1178 CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 1381 !1382 IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll '1383 1179 ! 1384 1180 CALL iom_put( 'sig1_pnorm' , zsig1_p ) … … 1386 1182 1387 1183 DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 1388 1389 IF ( lwp ) WRITE(numout,*) ' So what ??? '1390 1184 1391 1185 ENDIF … … 1396 1190 1397 1191 ! --- Recalculate Coriolis stress at last inner iteration 1398 DO jj = 2, jpj - 1 1399 DO ji = 2, jpi - 1 1192 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1400 1193 ! --- U-component 1401 1194 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & … … 1405 1198 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 1406 1199 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 1407 END DO 1408 END DO 1200 END_2D 1409 1201 ! 1410 1202 CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & … … 1421 1213 1422 1214 ! Recalculate internal forces (divergence of stress tensor) at last inner iteration 1423 DO jj = 2, jpj -11424 DO ji = 2, jpi - 1 1215 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1216 1425 1217 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 1426 1218 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & … … 1429 1221 & ) * 2._wp * r1_e1u(ji,jj) & 1430 1222 & ) * r1_e1e2u(ji,jj) 1223 1431 1224 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 1432 1225 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & … … 1435 1228 & ) * 2._wp * r1_e2v(ji,jj) & 1436 1229 & ) * r1_e1e2v(ji,jj) 1437 END DO 1438 END DO1230 1231 END_2D 1439 1232 1440 1233 CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. ) … … 1452 1245 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 1453 1246 ! 1454 DO jj = 2, jpj - 1 1455 DO ji = 2, jpi - 1 1456 ! 2D ice mass, snow mass, area transport arrays (X, Y) 1247 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 ! 2D ice mass, snow mass, area transport arrays (X, Y) 1248 1457 1249 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 1458 1250 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) … … 1467 1259 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 1468 1260 1469 END DO 1470 END DO 1471 1261 END_2D 1262 1472 1263 CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 1473 1264 & zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & … … 1489 1280 1490 1281 1491 1492 SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, & 1493 & prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 1494 1282 SUBROUTINE rhg_cvg_vp( kt, kitout, kitinn, kitinntot, kitoutmax, kitinnmax, kitinntotmax , & 1283 & pu, pv, pub, pvb, pub_outer, pvb_outer , & 1284 & pmt, pat_iu, pat_iv, puerr_max, pverr_max, pglob_area , & 1285 & prhsu, pAU, pBU, pCU, pDU, pEU, pFU , & 1286 & prhsv, pAV, pBV, pCV, pDV, pEV, pFV , & 1287 & pvel_res, pvel_diff ) 1288 !! 1495 1289 !!---------------------------------------------------------------------- 1496 1290 !! *** ROUTINE rhg_cvg_vp *** … … 1507 1301 !! 1508 1302 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 1303 !! 1509 1304 !!---------------------------------------------------------------------- 1510 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1511 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pmt ! now velocity and mass per unit area 1512 REAL(wp), INTENT(in) :: puerr_max, pverr_max ! absolute mean velocity difference 1513 REAL(wp), INTENT(in) :: pglob_area ! global ice area 1514 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients 1515 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV 1516 !! 1517 INTEGER :: it, idtime, istatus, ix_dim, iy_dim 1305 !! 1306 INTEGER , INTENT(in) :: kt, kitout, kitinn, kitinntot ! ocean model iterate, outer, inner and total n-iterations 1307 INTEGER , INTENT(in) :: kitoutmax, kitinnmax ! max number of outer & inner iterations 1308 INTEGER , INTENT(in) :: kitinntotmax ! max number of total sub-iterations 1309 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now & sub-iter-before velocities 1310 REAL(wp), DIMENSION(:,:), INTENT(in) :: pub_outer, pvb_outer ! velocities @before outer iterations 1311 REAL(wp), DIMENSION(:,:), INTENT(in) :: pmt, pat_iu, pat_iv ! mass at T-point, ice concentration at U&V 1312 REAL(wp), INTENT(in) :: puerr_max, pverr_max ! absolute mean velocity difference 1313 REAL(wp), INTENT(in) :: pglob_area ! global ice area 1314 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsu, pAU, pBU, pCU, pDU, pEU, pFU ! linear system coefficients 1315 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV, pFV 1316 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pvel_res ! velocity residual @last inner iteration 1317 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pvel_diff ! velocity difference @last outer iteration 1318 !! 1319 1320 INTEGER :: idtime, istatus, ix_dim, iy_dim 1518 1321 INTEGER :: ji, jj ! dummy loop indices 1519 REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 1520 REAL(wp) :: z1_pglob_area 1322 INTEGER :: it_inn_file, it_out_file 1323 REAL(wp) :: zu_res_mean, zv_res_mean, zvel_res_mean ! mean residuals of the linear system 1324 REAL(wp) :: zu_mad, zv_mad, zvel_mad ! mean absolute deviation, sub-iterates 1325 REAL(wp) :: zu_mad_outer, zv_mad_outer, zvel_mad_outer ! mean absolute deviation, outer-iterates 1326 REAL(wp) :: zvel_err_max, zmke, zu, zv ! local scalars 1327 REAL(wp) :: z1_pglob_area ! inverse global ice area 1328 1521 1329 REAL(wp), DIMENSION(jpi,jpj) :: zu_res, zv_res, zvel2 ! local arrays 1330 REAL(wp), DIMENSION(jpi,jpj) :: zu_diff, zv_diff ! local arrays 1522 1331 1523 1332 CHARACTER(len=20) :: clname 1524 1333 !!---------------------------------------------------------------------- 1525 1334 1335 1526 1336 IF( lwp ) THEN 1337 1527 1338 WRITE(numout,*) 1528 1339 WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 1529 1340 WRITE(numout,*) '~~~~~~~~~~~' 1530 WRITE(numout,*) ' kiter = : ', kiter 1531 WRITE(numout,*) ' kitermax = : ', kitermax 1341 WRITE(numout,*) ' kt = : ', kt 1342 WRITE(numout,*) ' kitout = : ', kitout 1343 WRITE(numout,*) ' kitinn = : ', kitinn 1344 WRITE(numout,*) ' kitinntot = : ', kitinntot 1345 WRITE(numout,*) ' kitoutmax (nn_vp_nout) = ', kitoutmax 1346 WRITE(numout,*) ' kitinnmax (nn_vp_ninn) = ', kitinnmax 1347 WRITE(numout,*) ' kitinntotmax (nn_nvp) = ', kitinntotmax 1348 WRITE(numout,*) 1349 1532 1350 ENDIF 1533 1351 1352 z1_pglob_area = 1._wp / pglob_area ! inverse global ice area 1353 1534 1354 ! create file 1535 IF( kt == nit000 .AND. kit er== 1 ) THEN1355 IF( kt == nit000 .AND. kitinntot == 1 ) THEN 1536 1356 ! 1537 1357 IF( lwm ) THEN … … 1545 1365 istatus = NF90_DEF_DIM( ncvgid, 'y' , jpj, iy_dim ) 1546 1366 1547 ! i suggest vel_dif instead 1548 istatus = NF90_DEF_VAR( ncvgid, 'u_res' , NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1549 istatus = NF90_DEF_VAR( ncvgid, 'v_res' , NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1550 istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1551 istatus = NF90_DEF_VAR( ncvgid, 'u_dif' , NF90_DOUBLE , (/ idtime /), nvarid_udif ) 1552 istatus = NF90_DEF_VAR( ncvgid, 'v_dif' , NF90_DOUBLE , (/ idtime /), nvarid_vdif ) 1553 istatus = NF90_DEF_VAR( ncvgid, 'vel_dif', NF90_DOUBLE , (/ idtime /), nvarid_veldif ) 1367 istatus = NF90_DEF_VAR( ncvgid, 'u_res' , NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1368 istatus = NF90_DEF_VAR( ncvgid, 'v_res' , NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1369 istatus = NF90_DEF_VAR( ncvgid, 'vel_res' , NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1370 1371 istatus = NF90_DEF_VAR( ncvgid, 'uerr_max_sub' , NF90_DOUBLE , (/ idtime /), nvarid_uerr_max ) 1372 istatus = NF90_DEF_VAR( ncvgid, 'verr_max_sub' , NF90_DOUBLE , (/ idtime /), nvarid_verr_max ) 1373 istatus = NF90_DEF_VAR( ncvgid, 'velerr_max_sub', NF90_DOUBLE , (/ idtime /), nvarid_velerr_max ) 1374 1375 istatus = NF90_DEF_VAR( ncvgid, 'umad_sub' , NF90_DOUBLE , (/ idtime /), nvarid_umad ) 1376 istatus = NF90_DEF_VAR( ncvgid, 'vmad_sub' , NF90_DOUBLE , (/ idtime /), nvarid_vmad ) 1377 istatus = NF90_DEF_VAR( ncvgid, 'velmad_sub' , NF90_DOUBLE , (/ idtime /), nvarid_velmad ) 1378 1379 istatus = NF90_DEF_VAR( ncvgid, 'umad_outer' , NF90_DOUBLE , (/ idtime /), nvarid_umad_outer ) 1380 istatus = NF90_DEF_VAR( ncvgid, 'vmad_outer' , NF90_DOUBLE , (/ idtime /), nvarid_vmad_outer ) 1381 istatus = NF90_DEF_VAR( ncvgid, 'velmad_outer' , NF90_DOUBLE , (/ idtime /), nvarid_velmad_outer ) 1382 1554 1383 istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke ) 1555 1384 1556 istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy)1557 istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy)1558 1559 1385 istatus = NF90_ENDDEF(ncvgid) 1560 1386 … … 1563 1389 ENDIF 1564 1390 1565 IF ( lwp ) WRITE(numout,*) ' File created ' 1566 1567 ! --- Max absolute velocity difference with previous iterate (zveldif) 1568 zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code 1569 ! if puerrmask and pverrmax are masked at 15% (TEST) 1570 1571 ! --- Mean residual and kinetic energy 1572 IF ( kiter == 1 ) THEN 1573 1574 zu_res_mean = 0._wp 1575 zv_res_mean = 0._wp 1576 zvelres = 0._wp 1577 zmke = 0._wp 1391 !------------------------------------------------------------ 1392 ! 1393 ! Max absolute velocity difference with previous sub-iterate 1394 ! ( zvel_err_max ) 1395 ! 1396 !------------------------------------------------------------ 1397 ! 1398 ! This comes from the criterion used to stop the iterative procedure 1399 zvel_err_max = 0.5_wp * ( puerr_max + pverr_max ) ! average of U- and V- maximum error over the whole domain 1400 1401 !---------------------------------------------- 1402 ! 1403 ! Mean-absolute-deviation (sub-iterates) 1404 ! ( zu_mad, zv_mad, zvel_mad) 1405 ! 1406 !---------------------------------------------- 1407 ! 1408 ! U 1409 zu_diff(:,:) = 0._wp 1410 1411 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1412 1413 zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * z1_pglob_area 1414 1415 END_2D 1416 1417 ! V 1418 zv_diff(:,:) = 0._wp 1419 1420 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1421 1422 zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * z1_pglob_area 1423 1424 END_2D 1425 1426 ! global sum & U-V average 1427 CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_diff, 'U', 1., zv_diff , 'V', 1. ) 1428 zu_mad = glob_sum( 'icedyn_rhg_vp : ', zu_diff ) 1429 zv_mad = glob_sum( 'icedyn_rhg_vp : ', zv_diff ) 1430 1431 zvel_mad = 0.5_wp * ( zu_mad + zv_mad ) 1432 1433 !----------------------------------------------- 1434 ! 1435 ! Mean-absolute-deviation (outer-iterates) 1436 ! ( zu_mad_outer, zv_mad_outer, zvel_mad_outer) 1437 ! 1438 !----------------------------------------------- 1439 ! 1440 IF ( kitinn == kitinnmax ) THEN ! only work at the end of outer iterates 1441 1442 ! * U 1443 zu_diff(:,:) = 0._wp 1444 1445 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1446 1447 zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub_outer(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * & 1448 & z1_pglob_area 1449 1450 END_2D 1451 1452 ! * V 1453 zv_diff(:,:) = 0._wp 1454 1455 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1456 1457 zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb_outer(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * & 1458 & z1_pglob_area 1459 1460 END_2D 1461 1462 ! Global ice-concentration, grid-cell-area weighted mean 1463 CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_diff, 'U', 1., zv_diff , 'V', 1. ) ! abs behaves as a scalar no ? 1464 1465 zu_mad_outer = glob_sum( 'icedyn_rhg_vp : ', zu_diff ) 1466 zv_mad_outer = glob_sum( 'icedyn_rhg_vp : ', zv_diff ) 1467 1468 ! Average of both U & V 1469 zvel_mad_outer = 0.5_wp * ( zu_mad_outer + zv_mad_outer ) 1470 1471 ENDIF 1472 1473 ! --- Spatially-resolved absolute difference to send back to main routine 1474 ! (last iteration only, T-point) 1475 1476 IF ( kitinntot == kitinntotmax) THEN 1477 1478 zu_diff(:,:) = 0._wp 1479 zv_diff(:,:) = 0._wp 1480 1481 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1482 1483 zu_diff(ji,jj) = ( ABS ( ( pu(ji-1,jj) - pub_outer(ji-1,jj) ) ) * umask(ji-1,jj,1) & 1484 & + ABS ( ( pu(ji,jj ) - pub_outer(ji,jj) ) ) * umask(ji,jj,1) ) & 1485 & / ( umask(ji-1,jj,1) + umask(ji,jj,1) ) 1486 1487 zv_diff(ji,jj) = ( ABS ( ( pv(ji,jj-1) - pvb_outer(ji,jj-1) ) ) * vmask(ji,jj-1,1) & 1488 & + ABS ( ( pv(ji,jj ) - pvb_outer(ji,jj) ) ) * vmask(ji,jj,1) & 1489 & / ( vmask(ji,jj-1,1) + vmask(ji,jj,1) ) ) 1490 1491 1492 END_2D 1493 1494 CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_diff, 'T', 1., zv_diff , 'T', 1. ) 1495 pvel_diff(:,:) = 0.5_wp * ( zu_diff(:,:) + zv_diff(:,:) ) 1578 1496 1579 1497 ELSE 1580 1498 1581 ! -- Mean residual (N/m^2), zu_res_mean 1582 ! Here we take the residual of the linear system (N/m^2), 1583 ! We define it as in mitgcm: square-root of area-weighted mean square residual 1584 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1585 ! i.e., how close we are to a solution 1586 1587 IF ( lwp ) WRITE(numout,*) ' TEST 1 ' 1588 1589 z1_pglob_area = 1._wp / pglob_area 1590 1591 zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 1592 1593 DO jj = 2, jpj - 1 1594 DO ji = 2, jpi - 1 1595 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1596 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1597 1598 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1599 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1600 1601 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area 1602 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 1603 1604 END DO 1605 END DO 1606 1607 IF ( lwp ) WRITE(numout,*) ' TEST 2 ' 1608 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 1609 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 1610 IF ( lwp ) WRITE(numout,*) ' TEST 3 ' 1611 zvelres = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1612 1613 IF ( lwp ) WRITE(numout,*) ' TEST 4 ' 1614 1615 ! -- Global mean kinetic energy per unit area (J/m2) 1616 zvel2(:,:) = 0._wp 1617 DO jj = 2, jpj - 1 1618 DO ji = 2, jpi - 1 1619 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1620 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1621 zvel2(ji,jj) = zu*zu + zv*zv ! square of ice velocity at T-point 1622 END DO 1623 END DO 1624 1625 IF ( lwp ) WRITE(numout,*) ' TEST 5 ' 1626 1627 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1628 1629 IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 1630 1631 ENDIF ! kiter 1499 pvel_diff(:,:) = 0._wp 1500 1501 ENDIF 1502 1503 !--------------------------------------- 1504 ! 1505 ! --- Mean residual & kinetic energy 1506 ! 1507 !--------------------------------------- 1508 1509 IF ( kitinntot == 1 ) THEN 1510 1511 zu_res_mean = 0._wp 1512 zv_res_mean = 0._wp 1513 zvel_res_mean = 0._wp 1514 zmke = 0._wp 1515 1516 ELSE 1517 1518 ! * Mean residual (N/m2) 1519 ! Here we take the residual of the linear system (N/m2), 1520 ! We define it as in mitgcm: global area-weighted mean of square-root residual 1521 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1522 ! i.e., how close we are to a solution 1523 zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 1524 1525 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1526 1527 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1528 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1529 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1530 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1531 1532 ! zu_res(ji,jj) = pFU(ji,jj) - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) 1533 ! zv_res(ji,jj) = pFV(ji,jj) - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) 1534 1535 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * pat_iu(ji,jj) * e1e2u(ji,jj) * z1_pglob_area 1536 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * pat_iv(ji,jj) * e1e2v(ji,jj) * z1_pglob_area 1537 1538 END_2D 1539 1540 ! Global ice-concentration, grid-cell-area weighted mean 1541 CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res, 'U', 1., zv_res , 'V', 1. ) 1542 1543 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 1544 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 1545 zvel_res_mean = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1546 1547 ! --- Global mean kinetic energy per unit area (J/m2) 1548 zvel2(:,:) = 0._wp 1549 1550 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1632 1551 1633 ! ! ==================== ! 1634 1635 ! time 1636 it = ( kt - nit000 ) * kitermax + kiter 1637 1638 1552 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1553 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1554 zvel2(ji,jj) = zu*zu + zv*zv ! square of ice velocity at T-point 1555 1556 END_2D 1557 1558 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1559 1560 ENDIF ! kitinntot 1561 1562 !--- Spatially-resolved residual at last iteration to send back to main routine (last iteration only) 1563 !--- Calculation @T-point 1564 1565 IF ( kitinntot == kitinntotmax) THEN 1566 1567 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1568 1569 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1570 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1571 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1572 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1573 1574 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) 1575 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) 1576 1577 END_2D 1578 1579 CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res, 'U', 1., zv_res , 'V', 1. ) 1580 1581 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 1582 1583 pvel_res(ji,jj) = 0.25_wp * ( zu_res(ji-1,jj) + zu_res(ji,jj) + zv_res(ji,jj-1) + zv_res(ji,jj) ) 1584 1585 END_2D 1586 CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. ) 1587 1588 ELSE 1589 1590 pvel_res(:,:) = 0._wp 1591 1592 ENDIF 1593 1594 ! ! ==================== ! 1595 1596 it_inn_file = ( kt - nit000 ) * kitinntotmax + kitinntot ! time step in the file 1597 it_out_file = ( kt - nit000 ) * kitoutmax + kitout 1598 1599 ! write variables 1639 1600 IF( lwm ) THEN 1640 ! write variables 1641 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 1642 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 1643 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) ! average of u- and v- residuals 1644 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) ! max U velocity difference, inner iterations 1645 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) ! max V velocity difference, inner iterations 1646 istatus = NF90_PUT_VAR( ncvgid, nvarid_veldif, (/zveldif/), (/it/), (/1/) ) ! max U or V velocity diff between subiterations 1647 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) ! mean kinetic energy 1648 1649 ! 1650 IF ( kiter == kitermax ) THEN 1651 WRITE(numout,*) ' Should plot the spatially dependent residual ' 1652 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) ) ! U-residual, spatially dependent 1653 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) ) ! V-residual, spatially dependent 1601 1602 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures , (/zu_res_mean/), (/it_inn_file/), (/1/) ) ! Residuals of the linear system, area weighted mean 1603 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres , (/zv_res_mean/), (/it_inn_file/), (/1/) ) ! 1604 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvel_res_mean/), (/it_inn_file/), (/1/) ) ! 1605 1606 istatus = NF90_PUT_VAR( ncvgid, nvarid_uerr_max , (/puerr_max/), (/it_inn_file/), (/1/) ) ! Max velocit_inn_filey error, sub-it_inn_fileerates 1607 istatus = NF90_PUT_VAR( ncvgid, nvarid_verr_max , (/pverr_max/), (/it_inn_file/), (/1/) ) ! 1608 istatus = NF90_PUT_VAR( ncvgid, nvarid_velerr_max, (/zvel_err_max/), (/it_inn_file/), (/1/) ) ! 1609 1610 istatus = NF90_PUT_VAR( ncvgid, nvarid_umad , (/zu_mad/) , (/it_inn_file/), (/1/) ) ! velocit_inn_filey MAD, area/sic-weighted, sub-it_inn_fileerates 1611 istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad , (/zv_mad/) , (/it_inn_file/), (/1/) ) ! 1612 istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad , (/zvel_mad/), (/it_inn_file/), (/1/) ) ! 1613 1614 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/kitinntot/), (/1/) ) ! mean kinetic energy 1615 1616 IF ( kitinn == kitinnmax ) THEN ! only print outer mad at the end of inner loop 1617 1618 istatus = NF90_PUT_VAR( ncvgid, nvarid_umad_outer , (/zu_mad_outer/) , (/it_out_file/), (/1/) ) ! velocity MAD, area/sic-weighted, outer-iterates 1619 istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad_outer , (/zv_mad_outer/) , (/it_out_file/), (/1/) ) ! 1620 istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad_outer , (/zvel_mad_outer/), (/it_out_file/), (/1/) ) ! 1621 1654 1622 ENDIF 1655 1623 1656 ! close file 1657 IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax ) istatus = NF90_CLOSE(ncvgid) 1624 IF( kt == nitend - nn_fsbc + 1 .AND. kitinntot == kitinntotmax ) istatus = NF90_CLOSE( ncvgid ) 1658 1625 ENDIF 1659 1626
Note: See TracChangeset
for help on using the changeset viewer.
