- Timestamp:
- 2021-11-28T18:59:49+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ICE/icedyn_rhg_vp.F90
r14433 r15548 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_xy 60 61 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 62 57 58 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice 59 60 !! * Substitutions 61 # include "do_loop_substitute.h90" 63 62 !!---------------------------------------------------------------------- 64 63 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 86 85 !! 87 86 !! f1(u) = g1(v) 88 !! f2(v) = g2( v)87 !! f2(v) = g2(u) 89 88 !! 90 89 !! The right-hand side (RHS) is explicit … … 139 138 ! 140 139 INTEGER :: ji, ji2, jj, jj2, jn ! dummy loop indices 141 INTEGER :: jter, i_out, i_inn!140 INTEGER :: i_out, i_inn, i_inn_tot ! 142 141 INTEGER :: ji_min, jj_min ! 143 142 INTEGER :: nn_zebra_vp ! number of zebra steps 144 143 145 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp)146 144 ! 147 145 REAL(wp) :: zrhoco ! rho0 * rn_cio … … 150 148 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 151 149 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 152 REAL(wp) :: zd eltat, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars153 REAL(wp) :: zp_del tastar_f!150 REAL(wp) :: zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 151 REAL(wp) :: zp_delstar_f ! 154 152 REAL(wp) :: zu_cV, zv_cU ! 155 153 REAL(wp) :: zfac, zfac1, zfac2, zfac3 … … 158 156 REAL(wp) :: zAA3, zw, ztau, zuerr_max, zverr_max 159 157 ! 160 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice161 158 REAL(wp), DIMENSION(jpi,jpj) :: za_iU , za_iV ! ice fraction on U/V points 162 159 REAL(wp), DIMENSION(jpi,jpj) :: zmU_t, zmV_t ! Acceleration term contribution to RHS 163 160 REAL(wp), DIMENSION(jpi,jpj) :: zmassU_t, zmassV_t ! Mass per unit area divided by time step 164 161 ! 165 REAL(wp), DIMENSION(jpi,jpj) :: zdelta star_t !Delta* at T-points166 REAL(wp), DIMENSION(jpi,jpj) :: zten _i ! Tension167 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 168 165 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zet ! Viscosity pre-factors at T points 169 166 REAL(wp), DIMENSION(jpi,jpj) :: zef ! Viscosity pre-factor at F point … … 193 190 REAL(wp), DIMENSION(jpi,jpj) :: zFU, zFU_prime, zBU_prime ! Rearranged linear system coefficients, U equation 194 191 REAL(wp), DIMENSION(jpi,jpj) :: zFV, zFV_prime, zBV_prime ! Rearranged linear system coefficients, V equation 195 REAL(wp), DIMENSION(jpi,jpj) :: zCU_prime, zCV_prime ! Rearranged linear system coefficients, V equation196 192 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) 197 193 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 198 194 ! 195 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 199 196 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! mask for lots of ice (1), little ice (0) 200 197 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence (1), no ice (0) … … 204 201 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 205 202 !! --- diags 206 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y203 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 207 204 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12, zs12f ! stress tensor components 208 205 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p … … 212 209 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp, zdiag_yatrp ! X/Y-component of area transport (m2/s, SIMIP) 213 210 214 215 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 216 214 217 215 !!---------------------------------------------------------------------------------------------------------------------- 218 ! DEBUG put all forcing terms to zero219 ! air-ice drag220 utau_ice(:,:) = 0._wp221 vtau_ice(:,:) = 0._wp222 ! coriolis223 ff_t(:,:) = 0._wp224 ! ice-ocean drag225 rn_cio = 0._wp226 ! ssh227 ! done line 330 !!! dont forget to act there228 ! END DEBUG229 216 230 217 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' … … 238 225 239 226 ! for diagnostics and convergence tests 240 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 241 DO jj = 1, jpj 242 DO ji = 1, jpi 243 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 244 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 245 END DO 246 END DO 227 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 228 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 229 END_2D 247 230 248 231 IF ( lp_zebra_vp ) THEN; nn_zebra_vp = 2 … … 264 247 IF( nn_rhg_chkcvg /= 0 ) THEN 265 248 266 ! ice area for global mean kinetic energy 267 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) ) 268 251 269 252 ENDIF … … 276 259 277 260 zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 278 zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp; 279 zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp; 280 zrhsu(:,:) = 0._wp; zrhsv(:,:) = 0._wp 281 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 282 264 283 265 !------------------------------------------------------------------------------! … … 289 271 CALL ice_strength ! strength at T points 290 272 291 !------------------------------ 292 ! -- F-mask (code from EVP) 293 !------------------------------ 294 ! MartinV: 295 ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip) 296 ! I am not sure the same recipe applies here 297 298 ! - ocean/land mask 299 DO jj = 1, jpj - 1 300 DO ji = 1, jpi - 1 301 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 302 END DO 303 END DO 304 305 ! Lateral boundary conditions on velocity (modify zfmask) 306 ! Can be computed once for all, at first time step, for all rheologies 307 DO jj = 2, jpj - 1 308 DO ji = 2, jpi - 1 ! vector opt. 309 IF( zfmask(ji,jj) == 0._wp ) THEN 310 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 311 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 312 ENDIF 313 END DO 314 END DO 315 DO jj = 2, jpj - 1 316 IF( zfmask(1,jj) == 0._wp ) THEN 317 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 273 !--------------------------- 274 ! -- F-mask (code from EVP) 275 !--------------------------- 276 IF( kt == nit000 ) THEN 277 ! MartinV: 278 ! In EVP routine, fimask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip) 279 ! I am not sure the same recipe applies here 280 281 ! - ocean/land mask 282 ALLOCATE( fimask(jpi,jpj) ) 283 IF( rn_ishlat == 0._wp ) THEN 284 DO_2D( 0, 0, 0, 0 ) 285 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 286 END_2D 287 ELSE 288 DO_2D( 0, 0, 0, 0 ) 289 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 290 ! Lateral boundary conditions on velocity (modify fimask) 291 IF( fimask(ji,jj) == 0._wp ) THEN 292 fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 293 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 294 ENDIF 295 END_2D 318 296 ENDIF 319 IF( zfmask(jpi,jj) == 0._wp ) THEN 320 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpi - 1,jj,1), umask(jpi,jj-1,1) ) ) 321 ENDIF 322 END DO 323 DO ji = 2, jpi - 1 324 IF( zfmask(ji,1) == 0._wp ) THEN 325 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 326 ENDIF 327 IF( zfmask(ji,jpj) == 0._wp ) THEN 328 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpj - 1,1) ) ) 329 ENDIF 330 END DO 331 332 CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp ) 297 298 CALL lbc_lnk( 'icedyn_rhg_vp', fimask, 'F', 1._wp ) 299 ENDIF 333 300 334 301 !---------------------------------------------------------------------------------------------------------- … … 340 307 ! embedded sea ice: compute representative ice top surface 341 308 ! non-embedded sea ice: use ocean surface for slope calculation 342 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 343 zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 344 345 zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:) ! Snow and ice mass at T-point 346 zmf(:,:) = zmt(:,:) * ff_t(:,:) ! Coriolis factor at T points (m*f) 347 348 DO jj = 2, jpj - 1 349 DO ji = 2, jpi - 1 350 351 ! Ice fraction at U-V points 352 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) 353 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) 354 355 ! Snow and ice mass at U-V points 356 zm1 = zmt(ji,jj) 357 zm2 = zmt(ji+1,jj) 358 zm3 = zmt(ji,jj+1) 359 360 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 361 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 362 363 ! Mass per unit area divided by time step 364 zmassU_t(ji,jj) = zmassU * r1_Dt_ice 365 zmassV_t(ji,jj) = zmassV * r1_Dt_ice 366 367 ! Acceleration term contribution to RHS (depends on velocity at previous time step) 368 zmU_t(ji,jj) = zmassU_t(ji,jj) * u_ice(ji,jj) 369 zmV_t(ji,jj) = zmassV_t(ji,jj) * v_ice(ji,jj) 370 371 ! Ocean currents at U-V points 372 v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 373 u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 374 375 ! Wind stress 376 ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj) 377 ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj) 378 379 ! Force due to sea surface tilt(- m*g*GRAD(ssh)) 380 zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 381 zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 382 383 ! Mask for ice presence (1) or absence (0) 384 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 385 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 386 387 ! Mask for lots of ice (1) or little ice (0) 388 IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 389 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 390 IF ( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 391 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 392 393 ! MV TEST DEBUG 394 IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji+1,jj) <= zmmin ) .AND. & 395 & ( at_i(ji,jj) <= zamin .OR. at_i(ji+1,jj) <= zamin ) ) THEN ; zmsk01x(ji,jj) = 0._wp 396 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 397 398 IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji,jj+1) <= zmmin ) .AND. & 399 & ( at_i(ji,jj) <= zamin .OR. at_i(ji,jj+1) <= zamin ) ) THEN ; zmsk01y(ji,jj) = 0._wp 400 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 401 ! END MV TEST DEBUG 402 403 END DO 404 END DO 405 406 CALL iom_put( 'zmsk00x' , zmsk00x ) ! MV DEBUG 407 CALL iom_put( 'zmsk00y' , zmsk00y ) ! MV DEBUG 408 CALL iom_put( 'zmsk01x' , zmsk01x ) ! MV DEBUG 409 CALL iom_put( 'zmsk01y' , zmsk01y ) ! MV DEBUG 410 CALL iom_put( 'ztaux_ai' , ztaux_ai ) ! MV DEBUG 411 CALL iom_put( 'ztauy_ai' , ztauy_ai ) ! MV DEBUG 412 CALL iom_put( 'zspgU' , zspgU ) ! MV DEBUG 413 CALL iom_put( 'zspgV' , zspgV ) ! MV DEBUG 309 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 310 311 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 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) 314 END_2D 315 316 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 317 318 ! Ice fraction at U-V points 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) 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) 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) 326 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 327 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 328 329 ! Mass per unit area divided by time step 330 zmassU_t(ji,jj) = zmassU * r1_Dt_ice 331 zmassV_t(ji,jj) = zmassV * r1_Dt_ice 332 333 ! Acceleration term contribution to RHS (depends on velocity at previous time step) 334 zmU_t(ji,jj) = zmassU_t(ji,jj) * u_ice(ji,jj) 335 zmV_t(ji,jj) = zmassV_t(ji,jj) * v_ice(ji,jj) 336 337 ! Ocean currents at U-V points 338 v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 339 u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 340 341 ! Wind stress 342 ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj) 343 ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj) 344 345 ! Force due to sea surface tilt(- m*g*GRAD(ssh)) 346 zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 347 zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 348 349 ! Mask for ice presence (1) or absence (0) 350 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 351 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 352 353 ! Mask for lots of ice (1) or little ice (0) 354 IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 355 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 356 IF ( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 357 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 358 359 END_2D 414 360 415 361 !------------------------------------------------------------------------------! … … 422 368 zv_c(:,:) = v_ice(:,:) 423 369 424 jter= 0370 i_inn_tot = 0 425 371 426 372 DO i_out = 1, nn_vp_nout 427 373 428 IF( lwp ) WRITE(numout,*) ' outer loop i_out : ', i_out429 430 374 ! Velocities used in the non linear terms are the average of the past two iterates 431 ! u_it = 0.5 * ( u_{it-1} + u_{it-2} )375 ! u_it = 0.5 * ( u_{it-1} + u_{it-2} ) 432 376 ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009) 433 377 zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) ) … … 441 385 ! In the outer loop, one needs to update all RHS terms 442 386 ! with explicit velocity dependencies (viscosities, coriolis, ocean stress) 443 ! as a function of uc 444 ! 387 ! as a function of "current" velocities (uc, vc) 445 388 446 389 !------------------------------------------ … … 449 392 450 393 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 451 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 452 DO ji = 1, jpi - 1 453 454 ! shear at F points 455 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) & 456 & + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 457 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 458 459 END DO 460 END DO 461 462 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 463 CALL iom_put( 'zds' , zds ) ! MV DEBUG 464 465 IF( lwp ) WRITE(numout,*) ' outer loop 1a i_out : ', i_out 466 467 !DO jj = 2, jpj - 1 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 468 ! DO ji = 2, jpi - 1 ! 469 470 ! MV DEBUG 471 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 472 DO ji = 2, jpi ! 473 ! END MV DEBUG 474 475 ! shear**2 at T points (doc eq. A16) 476 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 477 & + 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) & 478 & ) * 0.25_wp * r1_e1e2t(ji,jj) 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 397 ! shear at F points 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) & 399 & + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 400 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 401 402 END_2D 403 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) 479 413 480 481 482 483 484 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 485 419 486 487 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) &488 &- ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &489 &) * r1_e1e2t(ji,jj)490 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 491 425 492 493 zdeltat= SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )426 ! delta at T points 427 zdeltat(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 494 428 495 496 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 497 431 498 ! P/deltaat T-points499 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) 500 434 501 502 zzt(ji,jj) = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)503 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 504 438 505 END DO 506 END DO 507 508 CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 509 510 CALL iom_put( 'zzt' , zzt ) ! MV DEBUG 511 CALL iom_put( 'zet' , zet ) ! MV DEBUG 512 CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 513 514 IF( lwp ) WRITE(numout,*) ' outer loop 1b i_out : ', i_out 515 516 DO jj = 1, jpj - 1 517 DO ji = 1, jpi - 1 518 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 519 445 ! P/delta* at F points 520 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) ) 521 447 522 448 ! Temporary zef factor at F-point 523 zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) 524 525 END DO 526 END DO 527 528 CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 529 CALL iom_put( 'zef' , zef ) ! MV DEBUG 530 IF( lwp ) WRITE(numout,*) ' outer loop 1c i_out : ', i_out 531 449 zef(ji,jj) = zp_delstar_f * r1_e1e2f(ji,jj) * z1_ecc2 * fimask(ji,jj) * 0.5_wp 450 451 END_2D 452 532 453 !--------------------------------------------------- 533 454 ! -- Ocean-ice drag and Coriolis RHS contributions 534 455 !--------------------------------------------------- 535 456 536 DO jj = 2, jpj - 1 537 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) 538 462 539 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 540 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) 541 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) 542 472 543 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 544 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) ) & 545 & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 546 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) ) & 547 & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 548 549 !--- Ocean-ice drag contributions to RHS 550 ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj) 551 ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj) 552 553 ! --- U-component of Coriolis Force (energy conserving formulation) 554 ! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done 555 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 556 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * zv_c(ji ,jj) + e1v(ji ,jj-1) * zv_c(ji ,jj-1) ) & 557 & + 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) ) ) 558 477 559 ! --- V-component of Coriolis Force (energy conserving formulation) 560 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 561 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * zu_c(ji,jj ) + e2u(ji-1,jj ) * zu_c(ji-1,jj ) ) & 562 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 563 564 END DO 565 END DO 566 567 IF( lwp ) WRITE(numout,*) ' outer loop 1d i_out : ', i_out 568 569 CALL lbc_lnk( 'icedyn_rhg_vp', zCwU , 'U', -1., zCwV, 'V', -1. ) 570 CALL lbc_lnk( 'icedyn_rhg_vp', zCorU, 'U', -1., zCorV, 'V', -1. ) 571 572 CALL iom_put( 'zCwU' , zCwU ) ! MV DEBUG 573 CALL iom_put( 'zCwV' , zCwV ) ! MV DEBUG 574 CALL iom_put( 'zCorU' , zCorU ) ! MV DEBUG 575 CALL iom_put( 'zCorV' , zCorV ) ! MV DEBUG 576 577 IF( lwp ) WRITE(numout,*) ' outer loop 1f i_out : ', i_out 578 579 ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system, 580 ! so there is no need for lbclnk on drag and coriolis 581 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 582 485 !------------------------------------- 583 486 ! -- Internal stress RHS contribution 584 487 !------------------------------------- 585 488 586 ! --- Stress contributions at T-points 587 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 588 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 589 493 590 ! sig1 contribution to RHS of U-equation at T-points 591 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) 592 497 593 ! sig2 contribution to RHS of U-equation at T-points 594 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) 595 596 ! sig1 contribution to RHS of V-equation at T-points 597 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 ) 598 599 ! sig2 contribution to RHS of V-equation at T-points 600 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) 601 602 END DO 603 END DO 604 605 CALL iom_put( 'zs1_rhsu' , zs1_rhsu ) ! MV DEBUG 606 CALL iom_put( 'zs2_rhsu' , zs2_rhsu ) ! MV DEBUG 607 CALL iom_put( 'zs1_rhsv' , zs1_rhsv ) ! MV DEBUG 608 CALL iom_put( 'zs2_rhsv' , zs2_rhsv ) ! MV DEBUG 609 610 ! a priori, no lbclnk, because rhsu is only used in the inner domain 611 612 ! --- Stress contributions at f-points 613 ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing 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 511 ! MV NOTE: I applied fimask on zds, by mimetism on EVP, but without deep understanding of what I was doing 614 512 ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 615 513 616 IF( lwp ) WRITE(numout,*) ' outer loop 2 i_out : ', i_out 617 618 DO jj = 1, jpj - 1 619 DO ji = 1, jpi - 1 620 621 ! sig12 contribution to RHS of U equation at F-points 622 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) 623 624 ! sig12 contribution to RHS of V equation at F-points 625 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) 626 627 END DO 628 END DO 629 630 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 631 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) 632 633 CALL iom_put( 'zs12_rhsu' , zs12_rhsu ) ! MV DEBUG 634 CALL iom_put( 'zs12_rhsv' , zs12_rhsv ) ! MV DEBUG 635 636 ! a priori, no lbclnk, because rhsu are only used in the inner domain 637 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 638 524 ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) 639 525 ! OPT: merge with next loop and use intermediate scalars for zf_rhsu 640 641 DO jj = 2, jpj - 1 642 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 643 528 ! --- U component of internal force contribution to RHS at U points 644 529 zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & … … 650 535 zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * & 651 536 & ( e1v(ji,jj) * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) ) & 652 & + 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) ) & 653 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) ) ) 654 655 END DO 656 END DO 657 658 CALL iom_put( 'zf_rhsu' , zf_rhsu ) ! MV DEBUG 659 CALL iom_put( 'zf_rhsv' , zf_rhsv ) ! MV DEBUG 539 540 END_2D 660 541 661 542 !--------------------------- … … 664 545 ! 665 546 ! OPT: could use intermediate scalars to reduce memory access 666 DO jj = 2, jpj - 1 667 DO ji = 2, jpi - 1 668 669 ! still miss ice ocean stress and acceleration contribution 670 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) 671 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) 672 673 END DO 674 END DO 675 676 CALL lbc_lnk( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V', -1.) 677 CALL lbc_lnk( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V', -1.) 678 CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V', -1.) 679 680 CALL iom_put( 'zmU_t' , zmU_t ) ! MV DEBUG 681 CALL iom_put( 'zmV_t' , zmV_t ) ! MV DEBUG 682 CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu ) ! MV DEBUG 683 CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv ) ! MV DEBUG 684 CALL iom_put( 'zrhsu' , zrhsu ) ! MV DEBUG 685 CALL iom_put( 'zrhsv' , zrhsv ) ! MV DEBUG 686 687 ! inner domain calculations -> no lbclnk 688 689 IF( lwp ) WRITE(numout,*) ' outer loop 4 i_out : ', i_out 690 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 691 554 !---------------------------------------------------------------------------------------! 692 555 ! … … 706 569 ! only zzt and zet are iteration-dependent, other only depend on scale factors 707 570 708 DO ji = 2, jpi - 1 ! internal domain do loop 709 DO jj = 2, jpj - 1 710 711 !------------------------------------- 712 ! -- Internal forces LHS contribution 713 !------------------------------------- 714 ! 715 ! --- U-component 716 ! 717 ! "T" factors (intermediate results) 718 ! 719 zfac = 0.5_wp * r1_e1e2u(ji,jj) 720 zfac1 = zfac * e2u(ji,jj) 721 zfac2 = zfac * r1_e2u(ji,jj) 722 zfac3 = 2._wp * zfac * r1_e1u(ji,jj) 723 724 zt12U = - zfac1 * zzt(ji+1,jj) 725 zt11U = zfac1 * zzt(ji,jj) 726 727 zt22U = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 728 zt21U = zfac2 * zet(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) 729 730 zt122U = - zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) 731 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) 732 594 733 734 735 736 zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U* r1_e2u(ji-1,jj)737 zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj)738 zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U* r1_e2u(ji+1,jj)739 740 zDU(ji,jj) =zt121U * r1_e1u(ji,jj-1)741 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) 742 604 743 744 745 746 747 748 749 zfac1 = zfac * e2v(ji,jj)750 751 752 753 zt12V = - zfac1 * zzt(ji,jj+1)754 zt11V = zfac1 * zzt(ji,jj)755 756 zt22V = zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1)757 zt21V = - zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj)758 759 zt122V = zfac3 * zef(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj)760 zt121V = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)761 762 763 764 765 zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V* r1_e1v(ji,jj-1)766 zBV(ji,jj) = ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) -( zt122V + zt121V ) * r1_e2v(ji,jj)767 zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V* r1_e1v(ji,jj+1)768 769 zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here770 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) 771 633 772 !----------------------------------------------------- 773 ! -- Ocean-ice drag and acceleration LHS contribution 774 !----------------------------------------------------- 775 zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 776 zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 777 778 END DO 779 END DO 780 781 CALL lbc_lnk( 'icedyn_rhg_vp', zAU , 'U', 1., zAV , 'V', 1. ) 782 CALL lbc_lnk( 'icedyn_rhg_vp', zBU , 'U', 1., zBV , 'V', 1. ) 783 CALL lbc_lnk( 'icedyn_rhg_vp', zCU , 'U', 1., zCV , 'V', 1. ) 784 CALL lbc_lnk( 'icedyn_rhg_vp', zDU , 'U', 1., zDV , 'V', 1. ) 785 CALL lbc_lnk( 'icedyn_rhg_vp', zEU , 'U', 1., zEV , 'V', 1. ) 786 787 CALL iom_put( 'zAU' , zAU ) ! MV DEBUG 788 CALL iom_put( 'zBU' , zBU ) ! MV DEBUG 789 CALL iom_put( 'zCU' , zCU ) ! MV DEBUG 790 CALL iom_put( 'zDU' , zDU ) ! MV DEBUG 791 CALL iom_put( 'zEU' , zEU ) ! MV DEBUG 792 CALL iom_put( 'zAV' , zAV ) ! MV DEBUG 793 CALL iom_put( 'zBV' , zBV ) ! MV DEBUG 794 CALL iom_put( 'zCV' , zCV ) ! MV DEBUG 795 CALL iom_put( 'zDV' , zDV ) ! MV DEBUG 796 CALL iom_put( 'zEV' , zEV ) ! MV DEBUG 797 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 798 642 !------------------------------------------------------------------------------! 799 643 ! … … 808 652 DO i_inn = 1, nn_vp_ninn ! inner loop iterations 809 653 810 IF( lwp ) WRITE(numout,*) ' inner loop 1 i_inn : ', i_inn811 812 654 !--- mitgcm computes initial value of residual here... 813 655 814 jter = jter + 1 815 ! l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 816 817 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 818 zv_b(:,:) = v_ice(:,:) 819 820 ! zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 821 ! 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(:,:) 822 661 823 662 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN … … 832 671 ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 833 672 834 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; 835 674 836 675 DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) … … 841 680 ELSE ; jj_min = 3 842 681 ENDIF 843 844 IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min845 682 846 683 DO jj = jj_min, jpj - 1, nn_zebra_vp … … 850 687 !------------------------ 851 688 DO ji = 2, jpi - 1 852 853 ! 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 854 693 ! see Zhang and Hibler, 1997, Appendix B 855 694 zAA3 = 0._wp … … 867 706 END DO 868 707 869 CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. )870 871 708 !--------------- 872 709 ! Forward sweep … … 874 711 DO jj = jj_min, jpj - 1, nn_zebra_vp 875 712 713 zBU_prime(2,jj) = zBU(2,jj) 714 zFU_prime(2,jj) = zFU(2,jj) 715 876 716 DO ji = 3, jpi - 1 877 717 … … 884 724 885 725 END DO 886 887 CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime, 'U', 1., zBU_prime, 'U', 1. ) 888 726 889 727 !----------------------------- 890 728 ! Backward sweep & relaxation … … 894 732 895 733 ! --- Backward sweep 734 896 735 ! last row 897 736 zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 898 737 u_ice(jpi-1,jj) = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 899 738 & * umask(jpi-1,jj,1) 900 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 901 741 zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 902 742 u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) & … … 904 744 END DO 905 745 906 !--- Relaxation 907 ! and velocity masking for little-ice and no-ice cases 746 !--- Relaxation and masking (for low-ice/no-ice cases) 908 747 DO ji = 2, jpi - 1 909 748 … … 917 756 918 757 END DO ! jj 758 759 CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1. ) 919 760 920 761 END DO ! zebra loop … … 932 773 !!! ZH97 explain it is critical for convergence speed 933 774 934 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; 935 776 936 777 DO jn = 1, nn_zebra_vp ! "zebra" loop … … 940 781 ENDIF 941 782 942 IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min943 944 783 DO ji = ji_min, jpi - 1, nn_zebra_vp 945 784 … … 949 788 DO jj = 2, jpj - 1 950 789 951 ! boundary condition substitution (check it is correctly applied !!!)790 ! subdomain boundary condition substitution (check it is correctly applied !!!) 952 791 ! see Zhang and Hibler, 1997, Appendix B 953 792 zAA3 = 0._wp … … 965 804 END DO 966 805 967 CALL lbc_lnk( 'icedyn_rhg_vp', zFV, 'V', 1.)968 969 806 !--------------- 970 807 ! Forward sweep … … 972 809 DO ji = ji_min, jpi - 1, nn_zebra_vp 973 810 974 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 975 815 976 816 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) … … 983 823 END DO 984 824 985 CALL lbc_lnk( 'icedyn_rhg_vp', zFV_prime, 'V', 1., zBV_prime, 'V', 1. )986 987 825 !----------------------------- 988 826 ! Backward sweep & relaxation … … 1003 841 END DO 1004 842 1005 ! --- Relaxation & masking (should it be now or later)843 ! --- Relaxation & masking 1006 844 DO jj = 2, jpj - 1 1007 845 … … 1015 853 1016 854 END DO ! ji 855 856 CALL lbc_lnk( 'icedyn_rhg_vp', v_ice, 'V', -1. ) 1017 857 1018 858 END DO ! zebra loop … … 1020 860 ENDIF ! ll_v_iterate 1021 861 1022 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 1023 863 1024 864 !-------------------------------------------------------------------------------------- … … 1031 871 ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 1032 872 ! then we must optimize the use of the mpp_max, which is prohibitive 1033 zuerr_max = 0._wp873 zuerr_max = 0._wp 1034 874 1035 875 IF ( ll_u_iterate .AND. MOD ( i_inn, nn_vp_chkcvg ) == 0 ) THEN … … 1037 877 ! - Maximum U-velocity difference 1038 878 zuerr(:,:) = 0._wp 1039 DO jj = 2, jpj - 1 1040 DO ji = 2, jpi - 1 1041 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 1042 END DO 1043 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 1044 885 zuerr_max = MAXVAL( zuerr ) 1045 886 CALL mpp_max( 'icedyn_rhg_evp', zuerr_max ) ! max over the global domain - damned! 1046 1047 ! - 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) 1048 889 IF ( i_inn > 1 .AND. zuerr_max > zuerr_max_vp ) THEN 1049 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 " … … 1068 909 ! - Maximum V-velocity difference 1069 910 zverr(:,:) = 0._wp 1070 DO jj = 2, jpj -11071 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 1072 913 zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1) 1073 END DO1074 END DO1075 914 915 END_2D 916 1076 917 zverr_max = MAXVAL( zverr ) 1077 918 CALL mpp_max( 'icedyn_rhg_evp', zverr_max ) ! max over the global domain - damned! … … 1098 939 ! 1099 940 !--------------------------------------------------------------------------------------- 1100 1101 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, & 1102 & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 1103 1104 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 1105 951 1106 952 END DO ! i_inn, end of inner loop … … 1108 954 END DO ! End of outer loop (i_out) ============================================================================================= 1109 955 1110 IF ( lwp ) WRITE(numout,*) ' We are out of outer loop ' 1111 1112 CALL lbc_lnk( 'icedyn_rhg_vp', zFU , 'U', 1., zFV , 'V', 1. ) 1113 CALL lbc_lnk( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zBV_prime , 'V', 1. ) 1114 CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime , 'U', 1., zFV_prime , 'V', 1. ) 1115 CALL lbc_lnk( 'icedyn_rhg_vp', zCU_prime , 'U', 1., zCV_prime , 'V', 1. ) 1116 1117 CALL iom_put( 'zFU' , zFU ) ! MV DEBUG 1118 CALL iom_put( 'zBU_prime' , zBU_prime ) ! MV DEBUG 1119 CALL iom_put( 'zCU_prime' , zCU_prime ) ! MV DEBUG 1120 CALL iom_put( 'zFU_prime' , zFU_prime ) ! MV DEBUG 1121 1122 CALL iom_put( 'zFV' , zFV ) ! MV DEBUG 1123 CALL iom_put( 'zBV_prime' , zBV_prime ) ! MV DEBUG 1124 CALL iom_put( 'zCV_prime' , zCV_prime ) ! MV DEBUG 1125 CALL iom_put( 'zFV_prime' , zFV_prime ) ! MV DEBUG 1126 1127 CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1128 1129 IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg ' 1130 IF( iom_use('uice_dbg' ) ) CALL iom_put( 'uice_dbg' , u_ice ) ! ice velocity u after solver 1131 IF( iom_use('vice_dbg' ) ) CALL iom_put( 'vice_dbg' , v_ice ) ! ice velocity v after solver 1132 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 1133 967 !------------------------------------------------------------------------------! 1134 968 ! 1135 ! --- Convergence diagnostics969 ! --- Recompute delta, shear and div (inputs for mechanical redistribution) 1136 970 ! 1137 971 !------------------------------------------------------------------------------! 1138 1139 IF( nn_rhg_chkcvg /= 0 ) THEN1140 1141 IF( iom_use('uice_cvg') ) THEN1142 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & ! ice velocity difference at last iteration1143 & ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )1144 ENDIF1145 1146 ENDIF1147 1148 ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead1149 DO jj = 2, jpj - 11150 DO ji = 2, jpi - 11151 1152 u_ice(ji,jj) = zmsk00x(ji,jj) &1153 & * ( zmsk01x(ji,jj) * u_oce(ji,jj) * 0.01_wp &1154 + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp )1155 1156 v_ice(ji,jj) = zmsk00y(ji,jj) &1157 & * ( zmsk01y(ji,jj) * v_oce(ji,jj) * 0.01_wp &1158 + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp )1159 1160 END DO1161 END DO1162 1163 CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )1164 1165 IF ( lwp ) WRITE(numout,*) ' Velocity replaced '1166 1167 ! END DEBUG1168 1169 !------------------------------------------------------------------------------!1170 !1171 ! --- Recompute delta, shear and div (inputs for mechanical redistribution)1172 !1173 !------------------------------------------------------------------------------!1174 972 ! 1175 973 ! MV OPT: subroutinize ? 1176 1177 DO jj = 1, jpj - 1 1178 DO ji = 1, jpi - 1 974 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) ! 1->jpj-1; 1->jpi-1 1179 975 1180 976 ! shear at F points 1181 977 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 1182 978 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 1183 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 1184 1185 END DO 1186 END DO 1187 1188 DO jj = 2, jpj - 1 1189 DO ji = 2, jpi - 1 ! 979 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 980 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 1190 984 1191 985 ! tension**2 at T points … … 1195 989 zdt2 = zdt * zdt 1196 990 1197 zten _i(ji,jj)= zdt991 ztens(ji,jj) = zdt 1198 992 1199 993 ! shear**2 at T points (doc eq. A16) … … 1202 996 & ) * 0.25_wp * r1_e1e2t(ji,jj) 1203 997 1204 ! shear at T points 1205 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 ) 1206 1003 1207 1004 ! divergence at T points … … 1209 1006 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 1210 1007 & ) * r1_e1e2t(ji,jj) 1008 1009 zdiv2 = pdivu_i(ji,jj) * pdivu_i(ji,jj) 1211 1010 1212 1011 ! delta at T points 1213 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )1012 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 1214 1013 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 1215 1216 !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 1217 pdelta_i(ji,jj) = zdelta + rn_creepl 1218 1219 END DO 1220 END DO 1221 1222 IF ( lwp ) WRITE(numout,*) ' Deformation recalculated ' 1014 1015 pdelta_i(ji,jj) = zdelta + rn_creepl ! * rswitch 1016 1017 END_2D 1223 1018 1224 1019 CALL lbc_lnk( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) … … 1237 1032 ! 1238 1033 ! ---- Sea ice stresses at T-points 1239 IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1240 1241 DO jj = 2, jpj - 1 1242 DO ji = 2, jpi - 1 1243 zp_deltastar_t(ji,jj) = strength(ji,jj) / pdelta_i(ji,jj) 1244 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) 1245 1042 zs1(ji,jj) = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 1246 zs2(ji,jj) = zfac * z1_ecc2 * zten _i(ji,jj)1247 zs12(ji,jj) = zfac * z1_ecc2 * pshear_i(ji,jj)1248 END DO 1249 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 1250 1047 1251 1048 CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. ) … … 1256 1053 IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1257 1054 1258 DO jj = 1, jpj - 1 1259 DO ji = 1, jpi - 1 1055 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) ! 1->jpj-1; 1->jpi-1 1260 1056 1261 1057 ! P/delta* at F points 1262 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) ) 1263 1059 1264 1060 ! s12 at F-points 1265 zs12f(ji,jj) = zp_del tastar_f * z1_ecc2 * zds(ji,jj)1061 zs12f(ji,jj) = zp_delstar_f * z1_ecc2 * zds(ji,jj) 1266 1062 1267 END DO 1268 END DO 1063 END_2D 1269 1064 1270 1065 CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. ) 1271 1066 1272 1067 ENDIF 1273 1274 IF ( lwp ) WRITE(numout,*) ' zs12f recalculated '1275 1068 1276 1069 ! … … 1286 1079 1287 1080 !--- Recalculate oceanic stress at last inner iteration 1288 DO jj = 2, jpj - 1 1289 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 1290 1082 1291 1083 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) … … 1303 1095 ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) ) 1304 1096 1305 END DO 1306 END DO 1097 END_2D 1307 1098 1308 1099 ! 1309 1100 CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, & 1310 ! & 1101 ! & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 1311 1102 ! 1312 1103 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) … … 1323 1114 ! --- Divergence, shear and strength --- ! 1324 1115 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 1325 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 1326 1117 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 1327 1118 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 1328 1119 1329 IF ( lwp ) WRITE(numout,*) 'Some terms recalculated '1330 1331 1120 ! --- Stress tensor invariants (SIMIP diags) --- ! 1332 1121 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN … … 1340 1129 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 1341 1130 ! 1342 DO jj = 2, jpj - 1 1343 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 1344 1132 ! Stress invariants 1345 zsig_I(ji,jj) = zs1(ji,jj) * 0.5_wp ! 1st invariant, aka average normal stress aka negative pressure 1346 zsig_II(ji,jj) = SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) ) ! 2nd invariant, aka maximum shear stress 1347 END DO 1348 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 1349 1136 1350 1137 CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.) … … 1356 1143 1357 1144 ENDIF 1358 1359 IF ( lwp ) WRITE(numout,*) 'SIMIP work done'1360 1145 1361 1146 ! --- Normalized stress tensor principal components --- ! … … 1370 1155 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 1371 1156 ! 1372 DO jj = 2, jpj -11373 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 1374 1159 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 1375 1160 ! and **deformations** at current iterates 1376 1161 ! following Lemieux & Dupont (2020) 1377 zfac = zp_deltastar_t(ji,jj) 1378 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 1379 zsig1 = 0._wp !!! FUCKING DEBUG TEST !!! 1380 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 1381 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 1382 1166 1383 1167 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 1384 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st invariant1385 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 1386 1170 1387 1171 ! Normalized principal stresses (used to display the ellipse) … … 1389 1173 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 1390 1174 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 1391 END DO 1392 END DO 1393 IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done' 1175 1176 END_2D 1394 1177 ! 1395 1178 CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 1396 !1397 IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll '1398 1179 ! 1399 1180 CALL iom_put( 'sig1_pnorm' , zsig1_p ) … … 1401 1182 1402 1183 DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 1403 1404 IF ( lwp ) WRITE(numout,*) ' So what ??? '1405 1184 1406 1185 ENDIF … … 1411 1190 1412 1191 ! --- Recalculate Coriolis stress at last inner iteration 1413 DO jj = 2, jpj - 1 1414 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 1415 1193 ! --- U-component 1416 1194 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & … … 1420 1198 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 1421 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) ) ) 1422 END DO 1423 END DO 1200 END_2D 1424 1201 ! 1425 1202 CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & … … 1436 1213 1437 1214 ! Recalculate internal forces (divergence of stress tensor) at last inner iteration 1438 DO jj = 2, jpj -11439 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 1440 1217 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 1441 1218 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & … … 1444 1221 & ) * 2._wp * r1_e1u(ji,jj) & 1445 1222 & ) * r1_e1e2u(ji,jj) 1223 1446 1224 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 1447 1225 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & … … 1450 1228 & ) * 2._wp * r1_e2v(ji,jj) & 1451 1229 & ) * r1_e1e2v(ji,jj) 1452 END DO 1453 END DO1230 1231 END_2D 1454 1232 1455 1233 CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. ) … … 1467 1245 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 1468 1246 ! 1469 DO jj = 2, jpj - 1 1470 DO ji = 2, jpi - 1 1471 ! 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 1472 1249 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 1473 1250 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) … … 1482 1259 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 1483 1260 1484 END DO 1485 END DO 1486 1261 END_2D 1262 1487 1263 CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 1488 1264 & zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & … … 1501 1277 ENDIF 1502 1278 1503 DEALLOCATE( zmsk00, zmsk15 )1504 1505 1279 END SUBROUTINE ice_dyn_rhg_vp 1506 1280 1507 1281 1508 1509 SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, & 1510 & prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 1511 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 !! 1512 1289 !!---------------------------------------------------------------------- 1513 1290 !! *** ROUTINE rhg_cvg_vp *** … … 1524 1301 !! 1525 1302 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 1303 !! 1526 1304 !!---------------------------------------------------------------------- 1527 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1528 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pmt ! now velocity and mass per unit area 1529 REAL(wp), INTENT(in) :: puerr_max, pverr_max ! absolute mean velocity difference 1530 REAL(wp), INTENT(in) :: pglob_area ! global ice area 1531 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients 1532 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV 1533 !! 1534 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 1535 1321 INTEGER :: ji, jj ! dummy loop indices 1536 REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 1537 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 1538 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 1539 1331 1540 1332 CHARACTER(len=20) :: clname 1541 1333 !!---------------------------------------------------------------------- 1542 1334 1335 1543 1336 IF( lwp ) THEN 1337 1544 1338 WRITE(numout,*) 1545 1339 WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 1546 1340 WRITE(numout,*) '~~~~~~~~~~~' 1547 WRITE(numout,*) ' kiter = : ', kiter 1548 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 1549 1350 ENDIF 1550 1351 1352 z1_pglob_area = 1._wp / pglob_area ! inverse global ice area 1353 1551 1354 ! create file 1552 IF( kt == nit000 .AND. kit er== 1 ) THEN1355 IF( kt == nit000 .AND. kitinntot == 1 ) THEN 1553 1356 ! 1554 1357 IF( lwm ) THEN … … 1562 1365 istatus = NF90_DEF_DIM( ncvgid, 'y' , jpj, iy_dim ) 1563 1366 1564 ! i suggest vel_dif instead 1565 istatus = NF90_DEF_VAR( ncvgid, 'u_res' , NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1566 istatus = NF90_DEF_VAR( ncvgid, 'v_res' , NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1567 istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1568 istatus = NF90_DEF_VAR( ncvgid, 'u_dif' , NF90_DOUBLE , (/ idtime /), nvarid_udif ) 1569 istatus = NF90_DEF_VAR( ncvgid, 'v_dif' , NF90_DOUBLE , (/ idtime /), nvarid_vdif ) 1570 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 1571 1383 istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke ) 1572 1384 1573 istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy)1574 istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy)1575 1576 1385 istatus = NF90_ENDDEF(ncvgid) 1577 1386 … … 1580 1389 ENDIF 1581 1390 1582 IF ( lwp ) WRITE(numout,*) ' File created ' 1583 1584 ! --- Max absolute velocity difference with previous iterate (zveldif) 1585 zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code 1586 ! if puerrmask and pverrmax are masked at 15% (TEST) 1587 1588 ! --- Mean residual and kinetic energy 1589 IF ( kiter == 1 ) THEN 1590 1591 zu_res_mean = 0._wp 1592 zv_res_mean = 0._wp 1593 zvelres = 0._wp 1594 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(:,:) ) 1595 1496 1596 1497 ELSE 1597 1498 1598 ! -- Mean residual (N/m^2), zu_res_mean 1599 ! Here we take the residual of the linear system (N/m^2), 1600 ! We define it as in mitgcm: square-root of area-weighted mean square residual 1601 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1602 ! i.e., how close we are to a solution 1603 1604 IF ( lwp ) WRITE(numout,*) ' TEST 1 ' 1605 1606 z1_pglob_area = 1._wp / pglob_area 1607 1608 zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 1609 1610 DO jj = 2, jpj - 1 1611 DO ji = 2, jpi - 1 1612 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1613 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1614 1615 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1616 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1617 1618 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area 1619 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 1620 1621 END DO 1622 END DO 1623 1624 IF ( lwp ) WRITE(numout,*) ' TEST 2 ' 1625 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 1626 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 1627 IF ( lwp ) WRITE(numout,*) ' TEST 3 ' 1628 zvelres = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1629 1630 IF ( lwp ) WRITE(numout,*) ' TEST 4 ' 1631 1632 ! -- Global mean kinetic energy per unit area (J/m2) 1633 zvel2(:,:) = 0._wp 1634 DO jj = 2, jpj - 1 1635 DO ji = 2, jpi - 1 1636 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1637 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1638 zvel2(ji,jj) = zu*zu + zv*zv ! square of ice velocity at T-point 1639 END DO 1640 END DO 1641 1642 IF ( lwp ) WRITE(numout,*) ' TEST 5 ' 1643 1644 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1645 1646 IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 1647 1648 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 1649 1551 1650 ! ! ==================== ! 1651 1652 ! time 1653 it = ( kt - 1 ) * kitermax + kiter 1654 1655 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 1656 1600 IF( lwm ) THEN 1657 ! write variables 1658 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 1659 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 1660 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) ! average of u- and v- residuals 1661 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) ! max U velocity difference, inner iterations 1662 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) ! max V velocity difference, inner iterations 1663 istatus = NF90_PUT_VAR( ncvgid, nvarid_veldif, (/zveldif/), (/it/), (/1/) ) ! max U or V velocity diff between subiterations 1664 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) ! mean kinetic energy 1665 1666 ! 1667 IF ( kiter == kitermax ) THEN 1668 WRITE(numout,*) ' Should plot the spatially dependent residual ' 1669 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) ) ! U-residual, spatially dependent 1670 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 1671 1622 ENDIF 1672 1623 1673 ! close file 1674 IF( kt == nitend ) istatus = NF90_CLOSE( ncvgid ) 1624 IF( kt == nitend - nn_fsbc + 1 .AND. kitinntot == kitinntotmax ) istatus = NF90_CLOSE( ncvgid ) 1675 1625 ENDIF 1676 1626
Note: See TracChangeset
for help on using the changeset viewer.