- Timestamp:
- 2020-11-27T17:26:33+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynvor.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900
- Property svn:externals
-
NEMO/branches/2020/tickets_icb_1900/src/OCE/DYN/dynvor.F90
r13237 r13899 80 80 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2u_2 ! = di(e2u)/2 used in T-point metric term calculation 81 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - - 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2 u)/(2*e1e2f) used in F-point metric term calculation83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1 v)/(2*e1e2f) - - - -82 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2v)/(2*e1e2f) used in F-point metric term calculation 83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1u)/(2*e1e2f) - - - - 84 84 85 85 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 … … 217 217 INTEGER :: ji, jj, jk ! dummy loop indices 218 218 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 219 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace220 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zwz ! 3D workspace219 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace 220 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 221 221 !!---------------------------------------------------------------------- 222 222 ! … … 231 231 CASE ( np_RVO ) !* relative vorticity 232 232 DO jk = 1, jpkm1 ! Horizontal slab 233 DO_2D _10_10233 DO_2D( 1, 0, 1, 0 ) 234 234 zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 235 235 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 236 236 END_2D 237 237 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 238 DO_2D _10_10238 DO_2D( 1, 0, 1, 0 ) 239 239 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 240 240 END_2D … … 246 246 CASE ( np_CRV ) !* Coriolis + relative vorticity 247 247 DO jk = 1, jpkm1 ! Horizontal slab 248 DO_2D _10_10248 DO_2D( 1, 0, 1, 0 ) ! relative vorticity 249 249 zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 250 250 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 251 251 END_2D 252 252 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 253 DO_2D _10_10253 DO_2D( 1, 0, 1, 0 ) 254 254 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 255 255 END_2D … … 269 269 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm) 270 270 CASE ( np_RVO ) !* relative vorticity 271 DO_2D _01_01271 DO_2D( 0, 1, 0, 1 ) 272 272 zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 273 273 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & … … 275 275 END_2D 276 276 CASE ( np_MET ) !* metric term 277 DO_2D _01_01277 DO_2D( 0, 1, 0, 1 ) 278 278 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 279 279 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) & … … 281 281 END_2D 282 282 CASE ( np_CRV ) !* Coriolis + relative vorticity 283 DO_2D _01_01283 DO_2D( 0, 1, 0, 1 ) 284 284 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 285 285 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) & … … 287 287 END_2D 288 288 CASE ( np_CME ) !* Coriolis + metric 289 DO_2D _01_01289 DO_2D( 0, 1, 0, 1 ) 290 290 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 291 291 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & … … 298 298 ! 299 299 ! !== compute and add the vorticity term trend =! 300 DO_2D _00_00300 DO_2D( 0, 0, 0, 0 ) 301 301 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 302 302 & * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) & … … 358 358 zwz(:,:) = ff_f(:,:) 359 359 CASE ( np_RVO ) !* relative vorticity 360 DO_2D _10_10360 DO_2D( 1, 0, 1, 0 ) 361 361 zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 362 362 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 363 363 END_2D 364 364 CASE ( np_MET ) !* metric term 365 DO_2D _10_10365 DO_2D( 1, 0, 1, 0 ) 366 366 zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 367 367 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 368 368 END_2D 369 369 CASE ( np_CRV ) !* Coriolis + relative vorticity 370 DO_2D _10_10370 DO_2D( 1, 0, 1, 0 ) 371 371 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 372 372 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 373 373 END_2D 374 374 CASE ( np_CME ) !* Coriolis + metric 375 DO_2D _10_10375 DO_2D( 1, 0, 1, 0 ) 376 376 zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 377 377 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 382 382 ! 383 383 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 384 DO_2D _10_10384 DO_2D( 1, 0, 1, 0 ) 385 385 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 386 386 END_2D … … 396 396 ENDIF 397 397 ! !== compute and add the vorticity term trend =! 398 DO_2D _00_00398 DO_2D( 0, 0, 0, 0 ) 399 399 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 400 400 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) … … 454 454 zwz(:,:) = ff_f(:,:) 455 455 CASE ( np_RVO ) !* relative vorticity 456 DO_2D _10_10456 DO_2D( 1, 0, 1, 0 ) 457 457 zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 458 458 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 459 459 END_2D 460 460 CASE ( np_MET ) !* metric term 461 DO_2D _10_10461 DO_2D( 1, 0, 1, 0 ) 462 462 zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 463 463 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 464 464 END_2D 465 465 CASE ( np_CRV ) !* Coriolis + relative vorticity 466 DO_2D _10_10466 DO_2D( 1, 0, 1, 0 ) 467 467 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 468 468 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 469 469 END_2D 470 470 CASE ( np_CME ) !* Coriolis + metric 471 DO_2D _10_10471 DO_2D( 1, 0, 1, 0 ) 472 472 zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 473 473 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 478 478 ! 479 479 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 480 DO_2D _10_10480 DO_2D( 1, 0, 1, 0 ) 481 481 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 482 482 END_2D … … 492 492 ENDIF 493 493 ! !== compute and add the vorticity term trend =! 494 DO_2D _00_00494 DO_2D( 0, 0, 0, 0 ) 495 495 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 496 496 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) … … 533 533 REAL(wp) :: zua, zva ! local scalars 534 534 REAL(wp) :: zmsk, ze3f ! local scalars 535 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f536 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse537 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zwz535 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f 536 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 537 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined 538 538 !!---------------------------------------------------------------------- 539 539 ! … … 550 550 SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point 551 551 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 552 DO_2D _10_10552 DO_2D( 1, 0, 1, 0 ) 553 553 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 554 554 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & … … 560 560 END_2D 561 561 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 562 DO_2D _10_10562 DO_2D( 1, 0, 1, 0 ) 563 563 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 564 564 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & … … 575 575 SELECT CASE( kvor ) !== vorticity considered ==! 576 576 CASE ( np_COR ) !* Coriolis (planetary vorticity) 577 DO_2D _10_10577 DO_2D( 1, 0, 1, 0 ) 578 578 zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj) 579 579 END_2D 580 580 CASE ( np_RVO ) !* relative vorticity 581 DO_2D _10_10581 DO_2D( 1, 0, 1, 0 ) 582 582 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 583 583 & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 584 584 END_2D 585 585 CASE ( np_MET ) !* metric term 586 DO_2D _10_10586 DO_2D( 1, 0, 1, 0 ) 587 587 zwz(ji,jj,jk) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 588 588 & - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 589 589 END_2D 590 590 CASE ( np_CRV ) !* Coriolis + relative vorticity 591 DO_2D _10_10591 DO_2D( 1, 0, 1, 0 ) 592 592 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 593 593 & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & … … 595 595 END_2D 596 596 CASE ( np_CME ) !* Coriolis + metric 597 DO_2D _10_10597 DO_2D( 1, 0, 1, 0 ) 598 598 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 599 599 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) … … 604 604 ! 605 605 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 606 DO_2D _10_10606 DO_2D( 1, 0, 1, 0 ) 607 607 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 608 608 END_2D … … 635 635 END DO 636 636 END DO 637 DO_2D _00_00637 DO_2D( 0, 0, 0, 0 ) 638 638 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 639 639 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) … … 677 677 REAL(wp) :: zua, zva ! local scalars 678 678 REAL(wp) :: zmsk, z1_e3t ! local scalars 679 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy680 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse681 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zwz679 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy 680 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 681 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined 682 682 !!---------------------------------------------------------------------- 683 683 ! … … 695 695 SELECT CASE( kvor ) !== vorticity considered ==! 696 696 CASE ( np_COR ) !* Coriolis (planetary vorticity) 697 DO_2D _10_10697 DO_2D( 1, 0, 1, 0 ) 698 698 zwz(ji,jj,jk) = ff_f(ji,jj) 699 699 END_2D 700 700 CASE ( np_RVO ) !* relative vorticity 701 DO_2D _10_10701 DO_2D( 1, 0, 1, 0 ) 702 702 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 703 703 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & … … 705 705 END_2D 706 706 CASE ( np_MET ) !* metric term 707 DO_2D _10_10707 DO_2D( 1, 0, 1, 0 ) 708 708 zwz(ji,jj,jk) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 709 709 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 710 710 END_2D 711 711 CASE ( np_CRV ) !* Coriolis + relative vorticity 712 DO_2D _10_10712 DO_2D( 1, 0, 1, 0 ) 713 713 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 714 714 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & … … 716 716 END_2D 717 717 CASE ( np_CME ) !* Coriolis + metric 718 DO_2D _10_10718 DO_2D( 1, 0, 1, 0 ) 719 719 zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 720 720 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 725 725 ! 726 726 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 727 DO_2D _10_10727 DO_2D( 1, 0, 1, 0 ) 728 728 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 729 729 END_2D … … 758 758 END DO 759 759 END DO 760 DO_2D _00_00760 DO_2D( 0, 0, 0, 0 ) 761 761 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 762 762 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) … … 818 818 IF(lwp) WRITE(numout,*) ' change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 819 819 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 820 DO_3D _10_10(1, jpk )820 DO_3D( 1, 0, 1, 0, 1, jpk ) 821 821 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 822 822 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp … … 857 857 CASE( np_ENT ) !* T-point metric term : pre-compute di(e2u)/2 and dj(e1v)/2 858 858 ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) ) 859 DO_2D _00_00859 DO_2D( 0, 0, 0, 0 ) 860 860 di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj ) ) * 0.5_wp 861 861 dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp … … 865 865 CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f) 866 866 ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) ) 867 DO_2D _10_10867 DO_2D( 1, 0, 1, 0 ) 868 868 di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj ) - e2v(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 869 869 dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj)
Note: See TracChangeset
for help on using the changeset viewer.
