- Timestamp:
- 2018-01-04T13:30:03+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r6498 r9176 172 172 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 173 173 fsdepw_b(:,:,1) = 0.0_wp 174 174 !$OPM PARALLEL 175 175 DO jk = 2, jpk 176 !$OMP DO PRIVATE(zcoef) 176 177 DO jj = 1,jpj 177 178 DO ji = 1,jpi … … 189 190 END DO 190 191 END DO 192 !$OMP END DO 191 193 END DO 194 !$OPM PARALLEL 192 195 193 196 ! Before depth and Inverse of the local depth of the water column at u- and v- points … … 214 217 ENDIF 215 218 IF ( ln_vvl_zstar_at_eqtor ) THEN 219 !$OMP PARALLEL DO 216 220 DO jj = 1, jpj 217 221 DO ji = 1, jpi … … 273 277 !!---------------------------------------------------------------------- 274 278 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 275 REAL(wp), POINTER, DIMENSION(:,:) :: zht, z_scale, zwu, zwv, zhdiv279 REAL(wp), DIMENSION(jpi, jpj ) :: zht, z_scale, zwu, zwv, zhdiv 276 280 !! * Arguments 277 281 INTEGER, INTENT( in ) :: kt ! time step … … 285 289 !!---------------------------------------------------------------------- 286 290 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 287 291 ! CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 288 292 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 289 293 … … 306 310 ! z_star coordinate and barotropic z-tilde part ! 307 311 ! ! --------------------------------------------- ! 308 309 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 312 !$OMP PARALLEL 313 !$OMP DO 314 DO jj = 1, jpj 315 z_scale(:,jj) = ( ssha(:,jj) - sshb(:,jj) ) * ssmask(:,jj) / ( ht_0(:,jj) + sshn(:,jj) + 1. - ssmask(:,jj) ) 316 END DO 317 !$OMP END DO 318 !$OMP DO 310 319 DO jk = 1, jpkm1 311 320 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) 312 321 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 313 322 END DO 323 !$OMP END DO 324 !$OMP END PARALLEL 314 325 315 326 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! … … 323 334 zhdiv(:,:) = 0. 324 335 zht(:,:) = 0. 336 !$OMP PARALLEL DO REDUCTION(+:zhdiv, zht) 325 337 DO jk = 1, jpkm1 326 338 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 327 339 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 328 340 END DO 329 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 341 !$OMP DO 342 DO jj = 1, jpj 343 zhdiv(:,jj) = zhdiv(:,jj) / ( zht(:,jj) + 1. - tmask_i(:,jj) ) 344 END DO 345 !$OMP END DO 330 346 331 347 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) … … 333 349 IF( ln_vvl_ztilde ) THEN 334 350 IF( kt .GT. nit000 ) THEN 351 !$OMP PARALLEL DO 335 352 DO jk = 1, jpkm1 336 353 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & … … 347 364 ! ---------------------------------- 348 365 IF( ln_vvl_ztilde ) THEN ! z_tilde case 366 !$OMP PARALLEL DO 349 367 DO jk = 1, jpkm1 350 368 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 351 369 END DO 352 370 ELSE ! layer case 371 !$OMP PARALLEL DO 353 372 DO jk = 1, jpkm1 354 373 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) … … 359 378 ! ------------------ 360 379 IF( ln_vvl_ztilde ) THEN 380 !$OMP PARALLEL DO 361 381 DO jk = 1, jpk 362 382 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) … … 369 389 zwv(:,:) = 0.0_wp 370 390 ! a - first derivative: diffusive fluxes 391 !$OMP PARALLEL 392 !$OMP DO 371 393 DO jk = 1, jpkm1 372 394 DO jj = 1, jpjm1 … … 376 398 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 377 399 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 400 END DO 401 END DO 402 END DO 403 !$OMP END DO 404 !$OMP DO REDUCTION(+:zwu, zwv) 405 DO jk = 1, jpkm1 406 DO jj = 1, jpjm1 407 DO ji = 1, fs_jpim1 ! vector opt. 378 408 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 379 409 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) … … 381 411 END DO 382 412 END DO 413 !$OMP END DO 414 !$OMP END PARALLEL 383 415 ! b - correction for last oceanic u-v points 416 !$OMP PARALLEL DO 384 417 DO jj = 1, jpj 385 418 DO ji = 1, jpi … … 389 422 END DO 390 423 ! c - second derivative: divergence of diffusive fluxes 424 !$OMP PARALLEL DO 391 425 DO jk = 1, jpkm1 392 426 DO jj = 2, jpjm1 … … 413 447 ENDIF 414 448 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 415 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 416 449 !$OMP PARALLEL DO 450 DO jk = 1, jpk 451 tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) 452 ENDDO 417 453 ! Maximum deformation control 418 454 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 419 455 ze3t(:,:,jpk) = 0.0_wp 420 DO jk = 1, jpkm1 421 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 422 END DO 456 !$OMP PARALLEL DO 457 DO jk = 1, jpkm1 458 DO jj = 1, jpj 459 !dir$ IVDEP 460 DO ji = 1, jpi 461 ze3t(ji, jj ,jk) = tilde_e3t_a(ji, jj ,jk) / e3t_0(ji, jj ,jk) * tmask(ji, jj ,jk) * tmask_i(ji, jj) 462 ENDDO 463 ENDDO 464 END DO 465 !$OMP PARALLEL 466 !$OMP DO 467 DO jk = 1, jpkm1 468 ze3t(:, : ,jk) = tilde_e3t_a(:, : ,jk) / e3t_0(:, : ,jk) * tmask_i(:, :) 469 END DO 470 !$OMP END DO 471 !$OMP DO 472 DO jk = 1, jpkm1 473 ze3t(:, : ,jk) = ze3t(:, : ,jk) * tmask(:, : ,jk) 474 END DO 475 !$OMP END DO 476 !$OMP END PARALLEL 477 423 478 z_tmax = MAXVAL( ze3t(:,:,:) ) 424 479 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain … … 448 503 ! - ML - end test 449 504 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 450 tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 451 tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 452 505 !$OMP PARALLEL DO 506 DO jk =1, jpk 507 tilde_e3t_a(:,:,jk) = MIN( tilde_e3t_a(:,:,jk), rn_zdef_max * e3t_0(:,:,jk) ) 508 tilde_e3t_a(:,:,jk) = MAX( tilde_e3t_a(:,:,jk), - rn_zdef_max * e3t_0(:,:,jk) ) 509 ENDDO 453 510 ! 454 511 ! "tilda" change in the after scale factor 455 512 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 513 !$OMP PARALLEL DO 456 514 DO jk = 1, jpkm1 457 515 dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) … … 464 522 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 465 523 zht(:,:) = 0. 524 !$OMP PARALLEL DO REDUCTION(+:zht) 466 525 DO jk = 1, jpkm1 467 526 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 468 527 END DO 469 528 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 529 !$OMP PARALLEL DO 470 530 DO jk = 1, jpkm1 471 531 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) … … 476 536 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 477 537 ! ! ---baroclinic part--------- ! 538 !$OMP PARALLEL DO 478 539 DO jk = 1, jpkm1 479 540 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) … … 491 552 ! 492 553 zht(:,:) = 0.0_wp 554 !$OMP PARALLEL DO REDUCTION(+:zht) 493 555 DO jk = 1, jpkm1 494 556 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) … … 499 561 ! 500 562 zht(:,:) = 0.0_wp 563 !$OMP PARALLEL DO REDUCTION(+:zht) 501 564 DO jk = 1, jpkm1 502 565 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk) … … 507 570 ! 508 571 zht(:,:) = 0.0_wp 572 !$OMP PARALLEL DO REDUCTION(+:zht) 509 573 DO jk = 1, jpkm1 510 574 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk) … … 540 604 hu_a(:,:) = 0._wp ! Ocean depth at U-points 541 605 hv_a(:,:) = 0._wp ! Ocean depth at V-points 606 542 607 DO jk = 1, jpkm1 543 608 hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) … … 545 610 END DO 546 611 ! ! Inverse of the local depth 547 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 548 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 549 550 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 612 !$OMP PARALLEL DO 613 DO jj = 1, jpj 614 hur_a(:,jj) = 1._wp / ( hu_a(:,jj) + 1._wp - umask_i(:,jj) ) * umask_i(:,jj) 615 hvr_a(:,jj) = 1._wp / ( hv_a(:,jj) + 1._wp - vmask_i(:,jj) ) * vmask_i(:,jj) 616 ENDDO 617 618 ! CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 551 619 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 552 620 … … 603 671 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 604 672 ELSE 605 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 673 !$OMP PARALLEL DO 674 DO jk = 1,jpk 675 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 606 676 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 677 ENDDO 607 678 ENDIF 608 679 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) … … 636 707 fsdepw_n(:,:,1) = 0.0_wp 637 708 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 638 709 !$OMP PARALLEL SHARED(jk) 639 710 DO jk = 2, jpk 711 !$OMP DO PRIVATE(zcoef) 640 712 DO jj = 1,jpj 641 713 DO ji = 1,jpi … … 649 721 END DO 650 722 END DO 723 !$OMP END DO 651 724 END DO 725 !$OMP END PARALLEL 652 726 653 727 ! Local depth and Inverse of the local depth of the water column at u- and v- points … … 663 737 ! -------------------------------------------- 664 738 ht(:,:) = 0. 739 !$OMP PARALLEL DO REDUCTION(+:ht) 665 740 DO jk = 1, jpkm1 666 741 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) … … 705 780 ! ! ------------------------------------- ! 706 781 ! horizontal surface weighted interpolation 782 !$OMP PARALLEL DO 707 783 DO jk = 1, jpk 708 784 DO jj = 1, jpjm1 … … 718 794 ! boundary conditions 719 795 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 720 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 796 !$OMP PARALLEL DO 797 DO jk = 1, jpk 798 pe3_out(:,:,jk) = pe3_out(:,:,jk) + e3u_0(:,:,jk) 799 ENDDO 721 800 ! ! ------------------------------------- ! 722 801 CASE( 'V' ) ! interpolation from T-point to V-point ! 723 802 ! ! ------------------------------------- ! 724 803 ! horizontal surface weighted interpolation 804 !$OMP PARALLEL DO 725 805 DO jk = 1, jpk 726 806 DO jj = 1, jpjm1 … … 736 816 ! boundary conditions 737 817 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 738 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 818 !$OMP PARALLEL DO 819 DO jk = 1, jpk 820 pe3_out(:,:,jk) = pe3_out(:,:,jk) + e3v_0(:,:,jk) 821 ENDDO 739 822 ! ! ------------------------------------- ! 740 823 CASE( 'F' ) ! interpolation from U-point to F-point ! 741 824 ! ! ------------------------------------- ! 742 825 ! horizontal surface weighted interpolation 826 !$OMP PARALLEL DO 743 827 DO jk = 1, jpk 744 828 DO jj = 1, jpjm1 … … 754 838 ! boundary conditions 755 839 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 756 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 840 !$OMP PARALLEL DO 841 DO jk = 1, jpk 842 pe3_out(:,:,jk) = pe3_out(:,:,jk) + e3f_0(:,:,jk) 843 ENDDO 757 844 ! ! ------------------------------------- ! 758 845 CASE( 'W' ) ! interpolation from T-point to W-point ! … … 761 848 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 762 849 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 850 !$OMP PARALLEL DO 763 851 DO jk = 2, jpk 764 852 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & … … 771 859 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 772 860 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 861 !$OMP PARALLEL DO 773 862 DO jk = 2, jpk 774 863 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & … … 781 870 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 782 871 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 872 !$OMP PARALLEL DO 783 873 DO jk = 2, jpk 784 874 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & … … 857 947 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 858 948 IF(lwp) write(numout,*) 'neuler is forced to 0' 949 !$OMP PARALLEL DO 859 950 DO jk=1,jpk 860 951 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
Note: See TracChangeset
for help on using the changeset viewer.
