Changeset 13982 for NEMO/trunk/src/OCE/TRA
- Timestamp:
- 2020-12-02T11:57:05+01:00 (3 years ago)
- Location:
- NEMO/trunk/src/OCE/TRA
- Files:
-
- 23 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/eosbn2.F90
r13497 r13982 39 39 !!---------------------------------------------------------------------- 40 40 USE dom_oce ! ocean space and time domain 41 USE domutl, ONLY : is_tile 41 42 USE phycst ! physical constants 42 43 USE stopar ! Stochastic T/S fluctuations … … 189 190 190 191 SUBROUTINE eos_insitu( pts, prd, pdep ) 192 !! 193 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 194 ! ! 2 : salinity [psu] 195 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 196 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 197 !! 198 CALL eos_insitu_t( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) ) 199 END SUBROUTINE eos_insitu 200 201 SUBROUTINE eos_insitu_t( pts, ktts, prd, ktrd, pdep, ktdep ) 191 202 !!---------------------------------------------------------------------- 192 203 !! *** ROUTINE eos_insitu *** … … 222 233 !! TEOS-10 Manual, 2010 223 234 !!---------------------------------------------------------------------- 224 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 225 ! ! 2 : salinity [psu] 226 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 227 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 235 INTEGER , INTENT(in ) :: ktts, ktrd, ktdep 236 REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 237 ! ! 2 : salinity [psu] 238 REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-] 239 REAL(wp), DIMENSION(A2D_T(ktdep),JPK ), INTENT(in ) :: pdep ! depth [m] 228 240 ! 229 241 INTEGER :: ji, jj, jk ! dummy loop indices … … 238 250 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 239 251 ! 240 DO_3D( 1, 1, 1, 1, 1, jpkm1 )252 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 241 253 ! 242 254 zh = pdep(ji,jj,jk) * r1_Z0 ! depth … … 274 286 CASE( np_seos ) !== simplified EOS ==! 275 287 ! 276 DO_3D( 1, 1, 1, 1, 1, jpkm1 )288 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 277 289 zt = pts (ji,jj,jk,jp_tem) - 10._wp 278 290 zs = pts (ji,jj,jk,jp_sal) - 35._wp … … 293 305 IF( ln_timing ) CALL timing_stop('eos-insitu') 294 306 ! 295 END SUBROUTINE eos_insitu 307 END SUBROUTINE eos_insitu_t 296 308 297 309 298 310 SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 311 !! 312 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 313 ! ! 2 : salinity [psu] 314 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 315 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prhop ! potential density (surface referenced) 316 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m] 317 !! 318 CALL eos_insitu_pot_t( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) ) 319 END SUBROUTINE eos_insitu_pot 320 321 322 SUBROUTINE eos_insitu_pot_t( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep ) 299 323 !!---------------------------------------------------------------------- 300 324 !! *** ROUTINE eos_insitu_pot *** … … 309 333 !! 310 334 !!---------------------------------------------------------------------- 311 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 312 ! ! 2 : salinity [psu] 313 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 314 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced) 315 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 335 INTEGER , INTENT(in ) :: ktts, ktrd, ktrhop, ktdep 336 REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 337 ! ! 2 : salinity [psu] 338 REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK ), INTENT( out) :: prd ! in situ density [-] 339 REAL(wp), DIMENSION(A2D_T(ktrhop),JPK ), INTENT( out) :: prhop ! potential density (surface referenced) 340 REAL(wp), DIMENSION(A2D_T(ktdep) ,JPK ), INTENT(in ) :: pdep ! depth [m] 316 341 ! 317 342 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices … … 338 363 END DO 339 364 ! 340 DO_3D( 1, 1, 1, 1, 1, jpkm1 )365 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 341 366 ! 342 367 ! compute density (2*nn_sto_eos) times: … … 388 413 ! Non-stochastic equation of state 389 414 ELSE 390 DO_3D( 1, 1, 1, 1, 1, jpkm1 )415 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 391 416 ! 392 417 zh = pdep(ji,jj,jk) * r1_Z0 ! depth … … 426 451 CASE( np_seos ) !== simplified EOS ==! 427 452 ! 428 DO_3D( 1, 1, 1, 1, 1, jpkm1 )453 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 429 454 zt = pts (ji,jj,jk,jp_tem) - 10._wp 430 455 zs = pts (ji,jj,jk,jp_sal) - 35._wp … … 444 469 END SELECT 445 470 ! 446 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 471 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', & 472 & tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 447 473 ! 448 474 IF( ln_timing ) CALL timing_stop('eos-pot') 449 475 ! 450 END SUBROUTINE eos_insitu_pot 476 END SUBROUTINE eos_insitu_pot_t 451 477 452 478 453 479 SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 480 !! 481 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 482 ! ! 2 : salinity [psu] 483 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m] 484 REAL(wp), DIMENSION(:,:) , INTENT( out) :: prd ! in situ density 485 !! 486 CALL eos_insitu_2d_t( pts, is_tile(pts), pdep, is_tile(pdep), prd, is_tile(prd) ) 487 END SUBROUTINE eos_insitu_2d 488 489 490 SUBROUTINE eos_insitu_2d_t( pts, ktts, pdep, ktdep, prd, ktrd ) 454 491 !!---------------------------------------------------------------------- 455 492 !! *** ROUTINE eos_insitu_2d *** … … 462 499 !! 463 500 !!---------------------------------------------------------------------- 464 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 465 ! ! 2 : salinity [psu] 466 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 467 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 501 INTEGER , INTENT(in ) :: ktts, ktdep, ktrd 502 REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 503 ! ! 2 : salinity [psu] 504 REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m] 505 REAL(wp), DIMENSION(A2D_T(ktrd) ), INTENT( out) :: prd ! in situ density 468 506 ! 469 507 INTEGER :: ji, jj, jk ! dummy loop indices … … 480 518 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 481 519 ! 482 DO_2D( 1, 1, 1, 1)520 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 483 521 ! 484 522 zh = pdep(ji,jj) * r1_Z0 ! depth … … 515 553 CASE( np_seos ) !== simplified EOS ==! 516 554 ! 517 DO_2D( 1, 1, 1, 1)555 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 518 556 ! 519 557 zt = pts (ji,jj,jp_tem) - 10._wp … … 535 573 IF( ln_timing ) CALL timing_stop('eos2d') 536 574 ! 537 END SUBROUTINE eos_insitu_2d 575 END SUBROUTINE eos_insitu_2d_t 538 576 539 577 540 578 SUBROUTINE rab_3d( pts, pab, Kmm ) 579 !! 580 INTEGER , INTENT(in ) :: Kmm ! time level index 581 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity 582 REAL(wp), DIMENSION(:,:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio 583 !! 584 CALL rab_3d_t( pts, is_tile(pts), pab, is_tile(pab), Kmm ) 585 END SUBROUTINE rab_3d 586 587 588 SUBROUTINE rab_3d_t( pts, ktts, pab, ktab, Kmm ) 541 589 !!---------------------------------------------------------------------- 542 590 !! *** ROUTINE rab_3d *** … … 548 596 !! ** Action : - pab : thermal/haline expansion ratio at T-points 549 597 !!---------------------------------------------------------------------- 550 INTEGER , INTENT(in ) :: Kmm ! time level index 551 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 552 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab ! thermal/haline expansion ratio 598 INTEGER , INTENT(in ) :: Kmm ! time level index 599 INTEGER , INTENT(in ) :: ktts, ktab 600 REAL(wp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in ) :: pts ! pot. temperature & salinity 601 REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio 553 602 ! 554 603 INTEGER :: ji, jj, jk ! dummy loop indices … … 563 612 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 564 613 ! 565 DO_3D( 1, 1, 1, 1, 1, jpkm1 )614 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 566 615 ! 567 616 zh = gdept(ji,jj,jk,Kmm) * r1_Z0 ! depth … … 616 665 CASE( np_seos ) !== simplified EOS ==! 617 666 ! 618 DO_3D( 1, 1, 1, 1, 1, jpkm1 )667 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 619 668 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 620 669 zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0) … … 641 690 IF( ln_timing ) CALL timing_stop('rab_3d') 642 691 ! 643 END SUBROUTINE rab_3d 692 END SUBROUTINE rab_3d_t 644 693 645 694 646 695 SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 696 !! 697 INTEGER , INTENT(in ) :: Kmm ! time level index 698 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity 699 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pdep ! depth [m] 700 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pab ! thermal/haline expansion ratio 701 !! 702 CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm) 703 END SUBROUTINE rab_2d 704 705 706 SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm ) 647 707 !!---------------------------------------------------------------------- 648 708 !! *** ROUTINE rab_2d *** … … 652 712 !! ** Action : - pab : thermal/haline expansion ratio at T-points 653 713 !!---------------------------------------------------------------------- 654 INTEGER , INTENT(in ) :: Kmm ! time level index 655 REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: pts ! pot. temperature & salinity 656 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 657 REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio 714 INTEGER , INTENT(in ) :: Kmm ! time level index 715 INTEGER , INTENT(in ) :: ktts, ktdep, ktab 716 REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! pot. temperature & salinity 717 REAL(wp), DIMENSION(A2D_T(ktdep) ), INTENT(in ) :: pdep ! depth [m] 718 REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT( out) :: pab ! thermal/haline expansion ratio 658 719 ! 659 720 INTEGER :: ji, jj, jk ! dummy loop indices … … 670 731 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 671 732 ! 672 DO_2D( 1, 1, 1, 1)733 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 673 734 ! 674 735 zh = pdep(ji,jj) * r1_Z0 ! depth … … 723 784 CASE( np_seos ) !== simplified EOS ==! 724 785 ! 725 DO_2D( 1, 1, 1, 1)786 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 726 787 ! 727 788 zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) … … 748 809 IF( ln_timing ) CALL timing_stop('rab_2d') 749 810 ! 750 END SUBROUTINE rab_2d 811 END SUBROUTINE rab_2d_t 751 812 752 813 … … 849 910 850 911 SUBROUTINE bn2( pts, pab, pn2, Kmm ) 912 !! 913 INTEGER , INTENT(in ) :: Kmm ! time level index 914 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu] 915 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1] 916 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2] 917 !! 918 CALL bn2_t( pts, pab, is_tile(pab), pn2, is_tile(pn2), Kmm ) 919 END SUBROUTINE bn2 920 921 922 SUBROUTINE bn2_t( pts, pab, ktab, pn2, ktn2, Kmm ) 851 923 !!---------------------------------------------------------------------- 852 924 !! *** ROUTINE bn2 *** … … 862 934 !! 863 935 !!---------------------------------------------------------------------- 864 INTEGER , INTENT(in ) :: Kmm ! time level index 865 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu] 866 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1] 867 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2] 936 INTEGER , INTENT(in ) :: Kmm ! time level index 937 INTEGER , INTENT(in ) :: ktab, ktn2 938 REAL(wp), DIMENSION(jpi,jpj, jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu] 939 REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1] 940 REAL(wp), DIMENSION(A2D_T(ktn2),JPK ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2] 868 941 ! 869 942 INTEGER :: ji, jj, jk ! dummy loop indices … … 873 946 IF( ln_timing ) CALL timing_start('bn2') 874 947 ! 875 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90948 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 ) ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90 876 949 zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) & 877 950 & / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) … … 889 962 IF( ln_timing ) CALL timing_stop('bn2') 890 963 ! 891 END SUBROUTINE bn2 964 END SUBROUTINE bn2_t 892 965 893 966 … … 949 1022 950 1023 951 SUBROUTINE eos_fzp_2d( psal, ptf, pdep ) 1024 SUBROUTINE eos_fzp_2d( psal, ptf, pdep ) 1025 !! 1026 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 1027 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m] 1028 REAL(wp), DIMENSION(:,:) , INTENT(out ) :: ptf ! freezing temperature [Celsius] 1029 !! 1030 CALL eos_fzp_2d_t( psal, ptf, is_tile(ptf), pdep ) 1031 END SUBROUTINE eos_fzp_2d 1032 1033 1034 SUBROUTINE eos_fzp_2d_t( psal, ptf, kttf, pdep ) 952 1035 !!---------------------------------------------------------------------- 953 1036 !! *** ROUTINE eos_fzp *** … … 961 1044 !! Reference : UNESCO tech. papers in the marine science no. 28. 1978 962 1045 !!---------------------------------------------------------------------- 963 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 964 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m] 965 REAL(wp), DIMENSION(jpi,jpj), INTENT(out ) :: ptf ! freezing temperature [Celsius] 1046 INTEGER , INTENT(in ) :: kttf 1047 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: psal ! salinity [psu] 1048 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ), OPTIONAL :: pdep ! depth [m] 1049 REAL(wp), DIMENSION(A2D_T(kttf)), INTENT(out ) :: ptf ! freezing temperature [Celsius] 966 1050 ! 967 1051 INTEGER :: ji, jj ! dummy loop indices … … 996 1080 END SELECT 997 1081 ! 998 END SUBROUTINE eos_fzp_2d 1082 END SUBROUTINE eos_fzp_2d_t 999 1083 1000 1084 -
NEMO/trunk/src/OCE/TRA/traadv.F90
r13237 r13982 18 18 USE oce ! ocean dynamics and active tracers 19 19 USE dom_oce ! ocean space and time domain 20 ! TEMP: [tiling] This change not necessary after extended haloes development 21 USE domain, ONLY : dom_tile 20 22 USE domvvl ! variable vertical scale factors 21 23 USE sbcwave ! wave module … … 23 25 USE traadv_cen ! centered scheme (tra_adv_cen routine) 24 26 USE traadv_fct ! FCT scheme (tra_adv_fct routine) 27 USE traadv_fct_lf ! FCT scheme (tra_adv_fct routine - loop fusion version) 25 28 USE traadv_mus ! MUSCL scheme (tra_adv_mus routine) 29 USE traadv_mus_lf ! MUSCL scheme (tra_adv_mus routine - loop fusion version) 26 30 USE traadv_ubs ! UBS scheme (tra_adv_ubs routine) 27 31 USE traadv_qck ! QUICKEST scheme (tra_adv_qck routine) … … 65 69 INTEGER, PARAMETER :: np_UBS = 4 ! 3rd order Upstream Biased Scheme 66 70 INTEGER, PARAMETER :: np_QCK = 5 ! QUICK scheme 67 71 72 !! * Substitutions 73 # include "do_loop_substitute.h90" 68 74 # include "domzgr_substitute.h90" 69 75 !!---------------------------------------------------------------------- … … 86 92 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 87 93 ! 88 INTEGER :: jk ! dummy loop index 89 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zuu, zvv, zww ! 3D workspace 90 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds 94 INTEGER :: ji, jj, jk ! dummy loop index 95 ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) if using XIOS (subdomain support) 96 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zuu, zvv, zww ! 3D workspace 97 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds 98 ! TEMP: [tiling] This change not necessary after extra haloes development 99 LOGICAL :: lskip 91 100 !!---------------------------------------------------------------------- 92 101 ! 93 102 IF( ln_timing ) CALL timing_start('tra_adv') 94 103 ! 95 ! !== effective transport ==! 96 zuu(:,:,jpk) = 0._wp 97 zvv(:,:,jpk) = 0._wp 98 zww(:,:,jpk) = 0._wp 99 IF( ln_wave .AND. ln_sdw ) THEN 100 DO jk = 1, jpkm1 ! eulerian transport + Stokes Drift 101 zuu(:,:,jk) = & 102 & e2u (:,:) * e3u(:,:,jk,Kmm) * ( uu(:,:,jk,Kmm) + usd(:,:,jk) ) 103 zvv(:,:,jk) = & 104 & e1v (:,:) * e3v(:,:,jk,Kmm) * ( vv(:,:,jk,Kmm) + vsd(:,:,jk) ) 105 zww(:,:,jk) = & 106 & e1e2t(:,:) * ( ww(:,:,jk) + wsd(:,:,jk) ) 107 END DO 108 ELSE 109 DO jk = 1, jpkm1 110 zuu(:,:,jk) = e2u (:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm) ! eulerian transport only 111 zvv(:,:,jk) = e1v (:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm) 112 zww(:,:,jk) = e1e2t(:,:) * ww(:,:,jk) 113 END DO 114 ENDIF 115 ! 116 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! add z-tilde and/or vvl corrections 117 zuu(:,:,:) = zuu(:,:,:) + un_td(:,:,:) 118 zvv(:,:,:) = zvv(:,:,:) + vn_td(:,:,:) 119 ENDIF 120 ! 121 zuu(:,:,jpk) = 0._wp ! no transport trough the bottom 122 zvv(:,:,jpk) = 0._wp 123 zww(:,:,jpk) = 0._wp 124 ! 125 IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad ) & 126 & CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm, Krhs ) ! add the eiv transport (if necessary) 127 ! 128 IF( ln_mle ) CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm ) ! add the mle transport (if necessary) 129 ! 130 CALL iom_put( "uocetr_eff", zuu ) ! output effective transport 131 CALL iom_put( "vocetr_eff", zvv ) 132 CALL iom_put( "wocetr_eff", zww ) 133 ! 134 !!gm ??? 135 CALL dia_ptr( kt, Kmm, zvv ) ! diagnose the effective MSF 136 !!gm ??? 137 ! 138 139 IF( l_trdtra ) THEN !* Save ta and sa trends 140 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 141 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 142 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 143 ENDIF 144 ! 145 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 146 ! 147 CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order 148 CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) 149 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 150 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 151 CASE ( np_MUS ) ! MUSCL 152 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 153 CASE ( np_UBS ) ! UBS 154 CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v ) 155 CASE ( np_QCK ) ! QUICKEST 156 CALL tra_adv_qck ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 157 ! 158 END SELECT 159 ! 160 IF( l_trdtra ) THEN ! save the advective trends for further diagnostics 161 DO jk = 1, jpkm1 162 ztrdt(:,:,jk) = pts(:,:,jk,jp_tem,Krhs) - ztrdt(:,:,jk) 163 ztrds(:,:,jk) = pts(:,:,jk,jp_sal,Krhs) - ztrds(:,:,jk) 164 END DO 165 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_totad, ztrdt ) 166 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_totad, ztrds ) 167 DEALLOCATE( ztrdt, ztrds ) 104 lskip = .FALSE. 105 106 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 107 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 108 ALLOCATE( zuu(jpi,jpj,jpk), zvv(jpi,jpj,jpk), zww(jpi,jpj,jpk) ) 109 ENDIF 110 111 ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia) 112 IF( nadv /= np_CEN .OR. (nadv == np_CEN .AND. nn_cen_h == 4) .OR. ln_ldfeiv_dia ) THEN 113 IF( ln_tile ) THEN 114 IF( ntile == 1 ) THEN 115 CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 116 ELSE 117 lskip = .TRUE. 118 ENDIF 119 ENDIF 120 ENDIF 121 IF( .NOT. lskip ) THEN 122 ! !== effective transport ==! 123 IF( ln_wave .AND. ln_sdw ) THEN 124 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 125 zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * ( uu(ji,jj,jk,Kmm) + usd(ji,jj,jk) ) 126 zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * ( vv(ji,jj,jk,Kmm) + vsd(ji,jj,jk) ) 127 zww(ji,jj,jk) = e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) ) 128 END_3D 129 ELSE 130 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 131 zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) ! eulerian transport only 132 zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) 133 zww(ji,jj,jk) = e1e2t(ji,jj) * ww(ji,jj,jk) 134 END_3D 135 ENDIF 136 ! 137 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! add z-tilde and/or vvl corrections 138 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 139 zuu(ji,jj,jk) = zuu(ji,jj,jk) + un_td(ji,jj,jk) 140 zvv(ji,jj,jk) = zvv(ji,jj,jk) + vn_td(ji,jj,jk) 141 END_3D 142 ENDIF 143 ! 144 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 145 zuu(ji,jj,jpk) = 0._wp ! no transport trough the bottom 146 zvv(ji,jj,jpk) = 0._wp 147 zww(ji,jj,jpk) = 0._wp 148 END_2D 149 ! 150 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 151 IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad ) & 152 & CALL ldf_eiv_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), & 153 & 'TRA', Kmm, Krhs ) ! add the eiv transport (if necessary) 154 ! 155 IF( ln_mle ) CALL tra_mle_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), & 156 & 'TRA', Kmm ) ! add the mle transport (if necessary) 157 ! 158 ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support) 159 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 160 CALL iom_put( "uocetr_eff", zuu ) ! output effective transport 161 CALL iom_put( "vocetr_eff", zvv ) 162 CALL iom_put( "wocetr_eff", zww ) 163 ENDIF 164 ! 165 !!gm ??? 166 ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support) 167 CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) ) ! diagnose the effective MSF 168 !!gm ??? 169 ! 170 171 IF( l_trdtra ) THEN !* Save ta and sa trends 172 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 173 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 174 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 175 ENDIF 176 ! 177 ! NOTE: [tiling-comms-merge] These lbc_lnk calls are still needed (pts in the zco case because zps_hde is not called in step, zuu/zvv/zww in all cases, I think because DO loop bounds need to be updated in DYN as done in TRA) 178 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 179 ! 180 CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order 181 IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kmm), 'T', 1. ) 182 CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) 183 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 184 IF (nn_hls.EQ.2) THEN 185 CALL lbc_lnk_multi( 'traadv', pts(:,:,:,:,Kbb), 'T', 1., pts(:,:,:,:,Kmm), 'T', 1.) 186 CALL lbc_lnk_multi( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1., zww(:,:,:), 'W', 1.) 187 #if defined key_loop_fusion 188 CALL tra_adv_fct_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 189 #else 190 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 191 #endif 192 ELSE 193 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 194 END IF 195 CASE ( np_MUS ) ! MUSCL 196 ! NOTE: [tiling-comms-merge] I added this lbc_lnk as it did not validate against the trunk when using ln_zco 197 IF (nn_hls.EQ.2) THEN 198 CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.) 199 #if defined key_loop_fusion 200 CALL tra_adv_mus_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 201 #else 202 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 203 #endif 204 ELSE 205 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 206 END IF 207 CASE ( np_UBS ) ! UBS 208 IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.) 209 CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v ) 210 CASE ( np_QCK ) ! QUICKEST 211 IF (nn_hls.EQ.2) THEN 212 CALL lbc_lnk_multi( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1.) 213 CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.) 214 END IF 215 CALL tra_adv_qck ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 216 ! 217 END SELECT 218 ! 219 IF( l_trdtra ) THEN ! save the advective trends for further diagnostics 220 DO jk = 1, jpkm1 221 ztrdt(:,:,jk) = pts(:,:,jk,jp_tem,Krhs) - ztrdt(:,:,jk) 222 ztrds(:,:,jk) = pts(:,:,jk,jp_sal,Krhs) - ztrds(:,:,jk) 223 END DO 224 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_totad, ztrdt ) 225 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_totad, ztrds ) 226 DEALLOCATE( ztrdt, ztrds ) 227 ENDIF 228 229 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia) 230 IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) 231 168 232 ENDIF 169 233 ! ! print mean trends (used for debugging) 170 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' adv - Ta: ', mask1=tmask, 234 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' adv - Ta: ', mask1=tmask, & 171 235 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 236 237 ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support) 238 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 239 DEALLOCATE( zuu, zvv, zww ) 240 ENDIF 172 241 ! 173 242 IF( ln_timing ) CALL timing_stop( 'tra_adv' ) -
NEMO/trunk/src/OCE/TRA/traadv_cen.F90
r13497 r13982 71 71 INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme) 72 72 INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme) 73 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 73 74 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 74 75 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 78 79 REAL(wp) :: zC2t_u, zC4t_u ! local scalars 79 80 REAL(wp) :: zC2t_v, zC4t_v ! - - 80 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwx, zwy, zwz, ztu, ztv, ztw81 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwx, zwy, zwz, ztu, ztv, ztw 81 82 !!---------------------------------------------------------------------- 82 83 ! 83 IF( kt == kit000 ) THEN 84 IF(lwp) WRITE(numout,*) 85 IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v 86 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ ' 84 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 85 IF( kt == kit000 ) THEN 86 IF(lwp) WRITE(numout,*) 87 IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v 88 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ ' 89 ENDIF 90 ! ! set local switches 91 l_trd = .FALSE. 92 l_hst = .FALSE. 93 l_ptr = .FALSE. 94 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 95 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 96 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 97 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 87 98 ENDIF 88 ! ! set local switches89 l_trd = .FALSE.90 l_hst = .FALSE.91 l_ptr = .FALSE.92 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.93 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.94 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &95 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.96 99 ! 97 100 ! … … 112 115 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 113 116 ztv(:,:,jpk) = 0._wp 114 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! masked gradient117 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) ! masked gradient 115 118 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 116 119 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 117 120 END_3D 118 CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond.121 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. 119 122 ! 120 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes123 DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 ) ! Horizontal advective fluxes 121 124 zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! C2 interpolation of T at u- & v-points (x2) 122 125 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) … … 128 131 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v 129 132 END_3D 130 CALL lbc_lnk_multi( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. )133 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. ) 131 134 ! 132 135 CASE DEFAULT … … 155 158 END_2D 156 159 ELSE ! no ice-shelf cavities (only ocean surface) 157 zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm) 160 DO_2D( 1, 1, 1, 1 ) 161 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm) 162 END_2D 158 163 ENDIF 159 164 ENDIF … … 171 176 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) ) 172 177 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) ) 173 END 174 ! ! "Poleward" heat and salt transports 178 ENDIF 179 ! ! "Poleward" heat and salt transports 175 180 IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) ) 176 181 ! ! heat and salt transport -
NEMO/trunk/src/OCE/TRA/traadv_fct.F90
r13497 r13982 34 34 PUBLIC tra_adv_fct ! called by traadv.F90 35 35 PUBLIC interp_4th_cpt ! called by traadv_cen.F90 36 PUBLIC tridia_solver ! called by traadv_fct_lf.F90 37 PUBLIC nonosc ! called by traadv_fct_lf.F90 - key_agrif 36 38 37 39 LOGICAL :: l_trd ! flag to compute trends … … 79 81 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 80 82 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 83 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 84 ! NOTE: [tiling-comms-merge] These were changed to INTENT(inout) but they are not modified, so it is reverted 81 85 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 82 86 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 83 87 ! 84 INTEGER :: ji, jj, jk, jn ! dummy loop indices 88 INTEGER :: ji, jj, jk, jn ! dummy loop indices 85 89 REAL(wp) :: ztra ! local scalar 86 90 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - - 87 91 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - - 88 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw92 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 89 93 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry 90 94 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup … … 92 96 !!---------------------------------------------------------------------- 93 97 ! 94 IF( kt == kit000 ) THEN 95 IF(lwp) WRITE(numout,*) 96 IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype 97 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 98 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 99 IF( kt == kit000 ) THEN 100 IF(lwp) WRITE(numout,*) 101 IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype 102 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 103 ENDIF 104 ! NOTE: [tiling-comms-merge] Bug fix- move array zeroing out of this IF block 105 ! 106 l_trd = .FALSE. ! set local switches 107 l_hst = .FALSE. 108 l_ptr = .FALSE. 109 ll_zAimp = .FALSE. 110 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 111 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 112 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 113 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 114 ! 98 115 ENDIF 116 99 117 !! -- init to 0 100 118 zwi(:,:,:) = 0._wp … … 108 126 ztw(:,:,:) = 0._wp 109 127 ! 110 l_trd = .FALSE. ! set local switches111 l_hst = .FALSE.112 l_ptr = .FALSE.113 ll_zAimp = .FALSE.114 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.115 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.116 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &117 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.118 !119 128 IF( l_trd .OR. l_hst ) THEN 120 ALLOCATE( ztrdx( jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )129 ALLOCATE( ztrdx(A2D(nn_hls),jpk), ztrdy(A2D(nn_hls),jpk), ztrdz(A2D(nn_hls),jpk) ) 121 130 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 122 131 ENDIF 123 132 ! 124 IF( l_ptr ) THEN 125 ALLOCATE( zptry( jpi,jpj,jpk) )133 IF( l_ptr ) THEN 134 ALLOCATE( zptry(A2D(nn_hls),jpk) ) 126 135 zptry(:,:,:) = 0._wp 127 136 ENDIF 128 ! ! surface & bottom value : flux set to zero one for all129 zwz(:,:, 1 ) = 0._wp130 zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp131 !132 zwi(:,:,:) = 0._wp133 137 ! 134 138 ! If adaptive vertical advection, check if it is needed on this PE at this time 135 139 IF( ln_zad_Aimp ) THEN 136 IF( MAXVAL( ABS( wi( :,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.140 IF( MAXVAL( ABS( wi(A2D(nn_hls),:) ) ) > 0._wp ) ll_zAimp = .TRUE. 137 141 END IF 138 142 ! If active adaptive vertical advection, build tridiagonal matrix 139 143 IF( ll_zAimp ) THEN 140 ALLOCATE(zwdia( jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))141 DO_3D( 0, 0, 0, 0, 1, jpkm1 )144 ALLOCATE(zwdia(A2D(nn_hls),jpk), zwinf(A2D(nn_hls),jpk), zwsup(A2D(nn_hls),jpk)) 145 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 142 146 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 143 147 & / e3t(ji,jj,jk,Krhs) … … 151 155 ! !== upstream advection with initial mass fluxes & intermediate update ==! 152 156 ! !* upstream tracer flux in the i and j direction 153 DO_3D( 1, 0, 1, 0, 1, jpkm1 )157 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 154 158 ! upstream scheme 155 159 zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) … … 178 182 ENDIF 179 183 ! 180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme184 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 181 185 ! ! total intermediate advective trends 182 186 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & … … 194 198 ! 195 199 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 196 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask)200 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 197 201 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 198 202 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 206 210 ! 207 211 END IF 208 ! 212 ! 209 213 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 210 214 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) … … 218 222 ! 219 223 CASE( 2 ) !- 2nd order centered 220 DO_3D( 1, 0, 1, 0, 1, jpkm1 )224 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 221 225 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk) 222 226 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk) … … 238 242 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 239 243 ! 240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! Horizontal advective fluxes244 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 241 245 zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points 242 246 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 243 ! ! C4 minus upstream advective fluxes 247 ! ! C4 minus upstream advective fluxes 244 248 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 245 249 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 246 250 END_3D 251 IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'traadv_fct', zwx, 'U', -1.0_wp, zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 247 252 ! 248 253 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 249 254 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 250 255 ztv(:,:,jpk) = 0._wp 251 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! 1st derivative (gradient)256 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) ! 1st derivative (gradient) 252 257 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 253 258 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 254 259 END_3D 255 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 260 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 261 ! 262 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 256 263 ! 257 264 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes … … 265 272 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 266 273 END_3D 274 IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 267 275 ! 268 276 END SELECT … … 271 279 ! 272 280 CASE( 2 ) !- 2nd order centered 273 DO_3D( 0, 0, 0, 0, 2, jpkm1 )281 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 274 282 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) & 275 283 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) … … 278 286 CASE( 4 ) !- 4th order COMPACT 279 287 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 280 DO_3D( 0, 0, 0, 0, 2, jpkm1 )288 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 281 289 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 282 290 END_3D … … 286 294 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 287 295 ENDIF 288 ! 296 ! 297 IF (nn_hls.EQ.1) THEN 298 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 299 ELSE 300 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 301 END IF 302 ! 303 IF (nn_hls.EQ.1) THEN 304 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 305 ELSE 306 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 307 END IF 308 ! 289 309 IF ( ll_zAimp ) THEN 290 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme310 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 291 311 ! ! total intermediate advective trends 292 312 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 293 313 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 294 314 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 295 ztw(ji,jj,jk) 315 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 296 316 END_3D 297 317 ! 298 318 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 299 319 ! 300 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask)320 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 301 321 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 302 322 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 303 zwz(ji,jj,jk) = 323 zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 304 324 END_3D 305 325 END IF 306 !307 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'W', 1.0_wp )308 326 ! 309 327 ! !== monotonicity algorithm ==! … … 334 352 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 335 353 END_3D 336 END IF 337 ! 354 END IF 355 ! NOTE: [tiling-comms-merge] I tested this 356 ! NOT TESTED - NEED l_trd OR l_hst TRUE 338 357 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport 339 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes 358 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes 340 359 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes 341 360 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! … … 350 369 ! 351 370 ENDIF 371 ! NOTE: [tiling-comms-merge] I tested this 372 ! NOT TESTED - NEED l_ptr TRUE 352 373 IF( l_ptr ) THEN ! "Poleward" transports 353 374 zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes … … 360 381 DEALLOCATE( zwdia, zwinf, zwsup ) 361 382 ENDIF 362 IF( l_trd .OR. l_hst ) THEN 383 IF( l_trd .OR. l_hst ) THEN 363 384 DEALLOCATE( ztrdx, ztrdy, ztrdz ) 364 385 ENDIF … … 383 404 !! in-space based differencing for fluid 384 405 !!---------------------------------------------------------------------- 385 INTEGER , INTENT(in ) :: Kmm ! time level index 386 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 387 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 388 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 406 INTEGER , INTENT(in ) :: Kmm ! time level index 407 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 408 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pbef ! before field 409 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(in ) :: paft ! after field 410 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 389 411 ! 390 412 INTEGER :: ji, jj, jk ! dummy loop indices … … 392 414 REAL(dp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 393 415 REAL(dp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 394 REAL(dp), DIMENSION( jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo416 REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo 395 417 !!---------------------------------------------------------------------- 396 418 ! … … 402 424 ! -------------------- 403 425 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 404 zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), & 405 & paft * tmask - zbig * ( 1._wp - tmask ) ) 406 zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), & 407 & paft * tmask + zbig * ( 1._wp - tmask ) ) 426 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 427 zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ), & 428 & paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 429 zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ), & 430 & paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 431 END_3D 408 432 409 433 DO jk = 1, jpkm1 410 434 ikm1 = MAX(jk-1,1) 411 DO_2D( 0, 0, 0, 0)435 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 412 436 413 437 ! search maximum in neighbourhood … … 439 463 END_2D 440 464 END DO 441 CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign)465 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign) 442 466 443 467 ! 3. monotonic flux in the i & j direction (paa & pbb) 444 468 ! ---------------------------------------- 445 DO_3D( 0, 0, 0, 0, 1, jpkm1 )469 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 446 470 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 447 471 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) … … 461 485 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 462 486 END_3D 463 CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1.0_wp , pbb, 'V', -1.0_wp ) ! lateral boundary condition (changed sign)464 487 ! 465 488 END SUBROUTINE nonosc … … 537 560 !!---------------------------------------------------------------------- 538 561 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point 539 REAL(wp),DIMENSION( jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point562 REAL(wp),DIMENSION(A2D(nn_hls) ,jpk), INTENT( out) :: pt_out ! field interpolated at w-point 540 563 ! 541 564 INTEGER :: ji, jj, jk ! dummy loop integers 542 565 INTEGER :: ikt, ikb ! local integers 543 REAL(wp),DIMENSION( jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt566 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwd, zwi, zws, zwrm, zwt 544 567 !!---------------------------------------------------------------------- 545 568 ! 546 569 ! !== build the three diagonal matrix & the RHS ==! 547 570 ! 548 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! interior (from jk=3 to jpk-1)571 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 549 572 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 550 573 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal … … 565 588 END IF 566 589 ! 567 DO_2D( 0, 0, 0, 0) ! 2nd order centered at top & bottom590 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2nd order centered at top & bottom 568 591 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 569 592 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point … … 582 605 ! !== tridiagonal solver ==! 583 606 ! 584 DO_2D( 0, 0, 0, 0) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1607 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 585 608 zwt(ji,jj,2) = zwd(ji,jj,2) 586 609 END_2D 587 DO_3D( 0, 0, 0, 0, 3, jpkm1 )610 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 588 611 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 589 612 END_3D 590 613 ! 591 DO_2D( 0, 0, 0, 0) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1614 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 592 615 pt_out(ji,jj,2) = zwrm(ji,jj,2) 593 616 END_2D 594 DO_3D( 0, 0, 0, 0, 3, jpkm1 )617 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 595 618 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 596 619 END_3D 597 620 598 DO_2D( 0, 0, 0, 0) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk621 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 599 622 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 600 623 END_2D 601 DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )624 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 602 625 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 603 626 END_3D … … 626 649 !! The 3d array zwt is used as a work space array. 627 650 !!---------------------------------------------------------------------- 628 REAL(wp),DIMENSION( :,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix629 REAL(wp),DIMENSION( :,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side630 REAL(wp),DIMENSION( :,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev)631 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level632 ! ! =0 pt at t-level651 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix 652 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pRHS ! Right-Hand-Side 653 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: pt_out !!gm field at level=F(klev) 654 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 655 ! ! =0 pt at t-level 633 656 INTEGER :: ji, jj, jk ! dummy loop integers 634 657 INTEGER :: kstart ! local indices 635 REAL(wp),DIMENSION( jpi,jpj,jpk) :: zwt ! 3D work array658 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwt ! 3D work array 636 659 !!---------------------------------------------------------------------- 637 660 ! 638 661 kstart = 1 + klev 639 662 ! 640 DO_2D( 0, 0, 0, 0) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1663 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 641 664 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 642 665 END_2D 643 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )666 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 644 667 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 645 668 END_3D 646 669 ! 647 DO_2D( 0, 0, 0, 0) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1670 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 648 671 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 649 672 END_2D 650 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )673 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 651 674 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 652 675 END_3D 653 676 654 DO_2D( 0, 0, 0, 0) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk677 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 655 678 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 656 679 END_2D 657 DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 )680 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, kstart, -1 ) 658 681 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 659 682 END_3D -
NEMO/trunk/src/OCE/TRA/traadv_mus.F90
r13497 r13982 81 81 LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl 82 82 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 83 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 83 84 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 84 85 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 88 89 REAL(wp) :: zu, z0u, zzwx, zw , zalpha ! local scalars 89 90 REAL(wp) :: zv, z0v, zzwy, z0w ! - - 90 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwx, zslpx ! 3D workspace91 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwy, zslpy ! - -91 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwx, zslpx ! 3D workspace 92 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwy, zslpy ! - - 92 93 !!---------------------------------------------------------------------- 93 94 ! 94 IF( kt == kit000 ) THEN 95 IF(lwp) WRITE(numout,*) 96 IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype 97 IF(lwp) WRITE(numout,*) ' : mixed up-stream ', ld_msc_ups 98 IF(lwp) WRITE(numout,*) '~~~~~~~' 99 IF(lwp) WRITE(numout,*) 100 ! 101 ! Upstream / MUSCL scheme indicator 102 ! 103 ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr ) 104 xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed 105 ! 106 IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked) 107 ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 108 upsmsk(:,:) = 0._wp ! not upstream by default 109 ! 110 DO jk = 1, jpkm1 111 xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed 112 & - MAX ( rnfmsk(:,:) * rnfmsk_z(jk), & ! =>0 near runoff mouths (& closed sea outflows) 113 & upsmsk(:,:) ) * tmask(:,:,jk) ! =>0 in some user defined area 114 END DO 115 ENDIF 116 ! 117 ENDIF 118 ! 119 l_trd = .FALSE. 120 l_hst = .FALSE. 121 l_ptr = .FALSE. 122 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 123 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 124 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 125 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 95 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 96 IF( kt == kit000 ) THEN 97 IF(lwp) WRITE(numout,*) 98 IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype 99 IF(lwp) WRITE(numout,*) ' : mixed up-stream ', ld_msc_ups 100 IF(lwp) WRITE(numout,*) '~~~~~~~' 101 IF(lwp) WRITE(numout,*) 102 ! 103 ! Upstream / MUSCL scheme indicator 104 ! 105 ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr ) 106 xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed 107 ! 108 IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked) 109 ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 110 upsmsk(:,:) = 0._wp ! not upstream by default 111 ! 112 DO jk = 1, jpkm1 113 xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed 114 & - MAX ( rnfmsk(:,:) * rnfmsk_z(jk), & ! =>0 near runoff mouths (& closed sea outflows) 115 & upsmsk(:,:) ) * tmask(:,:,jk) ! =>0 in some user defined area 116 END DO 117 ENDIF 118 ! 119 ENDIF 120 ! 121 l_trd = .FALSE. 122 l_hst = .FALSE. 123 l_ptr = .FALSE. 124 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 125 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 126 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 127 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 128 ENDIF 126 129 ! 127 130 DO jn = 1, kjpt !== loop over the tracers ==! … … 132 135 zwx(:,:,jpk) = 0._wp ! bottom values 133 136 zwy(:,:,jpk) = 0._wp 134 DO_3D( 1, 0, 1, 0, 1, jpkm1 )137 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 135 138 zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 136 139 zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 137 140 END_3D 138 141 ! lateral boundary conditions (changed sign) 139 CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )142 IF ( nn_hls.EQ.1 ) CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) 140 143 ! !-- Slopes of tracer 141 144 zslpx(:,:,jpk) = 0._wp ! bottom values 142 145 zslpy(:,:,jpk) = 0._wp 143 DO_3D( 0, 1, 0, 1, 1, jpkm1 )146 DO_3D( nn_hls-1, 1, nn_hls-1, 1, 1, jpkm1 ) 144 147 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) & 145 148 & * ( 0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) ) … … 148 151 END_3D 149 152 ! 150 DO_3D( 0, 1, 0, 1, 1, jpkm1 ) !-- Slopes limitation153 DO_3D( nn_hls-1, 1, nn_hls-1, 1, 1, jpkm1 ) !-- Slopes limitation 151 154 zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), & 152 155 & 2.*ABS( zwx (ji-1,jj,jk) ), & … … 157 160 END_3D 158 161 ! 159 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- MUSCL horizontal advective fluxes162 DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 ) !-- MUSCL horizontal advective fluxes 160 163 ! MUSCL fluxes 161 164 z0u = SIGN( 0.5_wp, pU(ji,jj,jk) ) … … 173 176 zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 174 177 END_3D 175 CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! lateral boundary conditions (changed sign)178 IF ( nn_hls.EQ.1 ) CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! lateral boundary conditions (changed sign) 176 179 ! 177 180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend … … 195 198 zwx(:,:, 1 ) = 0._wp ! surface & bottom boundary conditions 196 199 zwx(:,:,jpk) = 0._wp 197 DO jk = 2, jpkm1! interior values198 zwx( :,:,jk) = tmask(:,:,jk) * ( pt(:,:,jk-1,jn,Kbb) - pt(:,:,jk,jn,Kbb) )199 END DO200 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! interior values 201 zwx(ji,jj,jk) = tmask(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 202 END_3D 200 203 ! !-- Slopes of tracer 201 204 zslpx(:,:,1) = 0._wp ! surface values … … 223 226 END_2D 224 227 ELSE ! no cavities: only at the ocean surface 225 zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb) 228 DO_2D( 1, 1, 1, 1 ) 229 zwx(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 230 END_2D 226 231 ENDIF 227 232 ENDIF -
NEMO/trunk/src/OCE/TRA/traadv_qck.F90
r13497 r13982 91 91 INTEGER , INTENT(in ) :: kjpt ! number of tracers 92 92 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 93 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 94 ! NOTE: [tiling-comms-merge] These were changed to INTENT(inout) but they are not modified, so it is reverted 93 95 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume transport components 94 96 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 95 97 !!---------------------------------------------------------------------- 96 98 ! 97 IF( kt == kit000 ) THEN 98 IF(lwp) WRITE(numout,*) 99 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype 100 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 101 IF(lwp) WRITE(numout,*) 99 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 100 IF( kt == kit000 ) THEN 101 IF(lwp) WRITE(numout,*) 102 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype 103 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 104 IF(lwp) WRITE(numout,*) 105 ENDIF 106 ! 107 l_trd = .FALSE. 108 l_ptr = .FALSE. 109 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 110 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 102 111 ENDIF 103 !104 l_trd = .FALSE.105 l_ptr = .FALSE.106 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.107 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.108 !109 112 ! 110 113 ! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme … … 127 130 INTEGER , INTENT(in ) :: kjpt ! number of tracers 128 131 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 132 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 129 133 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU ! i-velocity components 130 134 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation … … 132 136 INTEGER :: ji, jj, jk, jn ! dummy loop indices 133 137 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 134 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwx, zfu, zfc, zfd138 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwx, zfu, zfc, zfd 135 139 !---------------------------------------------------------------------- 136 140 ! … … 142 146 ! 143 147 !!gm why not using a SHIFT instruction... 144 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !--- Computation of the ustream and downstream value of the tracer and the mask148 DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 ) !--- Computation of the ustream and downstream value of the tracer and the mask 145 149 zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb) ! Upstream in the x-direction for the tracer 146 150 zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb) ! Downstream in the x-direction for the tracer 147 151 END_3D 148 CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions152 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 149 153 150 154 ! 151 155 ! Horizontal advective fluxes 152 156 ! --------------------------- 153 DO_3D( 0, 0, 0, 0, 1, jpkm1 )157 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 154 158 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 155 159 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 156 160 END_3D 157 161 ! 158 DO_3D( 0, 0, 0, 0, 1, jpkm1 )162 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 159 163 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 160 164 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) … … 164 168 END_3D 165 169 !--- Lateral boundary conditions 166 CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwx(:,:,:), 'T', 1.0_wp )170 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwx(:,:,:), 'T', 1.0_wp ) 167 171 168 172 !--- QUICKEST scheme … … 170 174 ! 171 175 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 172 DO_3D( 0, 0, 0, 0, 1, jpkm1 )176 DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 ) 173 177 zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 174 178 END_3D 175 CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions179 IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 176 180 177 181 ! 178 182 ! Tracer flux on the x-direction 179 DO jk = 1, jpkm1 180 ! 181 DO_2D( 0, 0, 0, 0 ) 182 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 183 !--- If the second ustream point is a land point 184 !--- the flux is computed by the 1st order UPWIND scheme 185 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 186 zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 187 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk) 188 END_2D 189 END DO 190 ! 191 CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 183 DO_3D( 0, 0, 1, 0, 1, jpkm1 ) 184 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 185 !--- If the second ustream point is a land point 186 !--- the flux is computed by the 1st order UPWIND scheme 187 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 188 zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 189 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk) 190 END_3D 192 191 ! 193 192 ! Computation of the trend … … 216 215 INTEGER , INTENT(in ) :: kjpt ! number of tracers 217 216 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 217 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 218 218 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pV ! j-velocity components 219 219 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation … … 221 221 INTEGER :: ji, jj, jk, jn ! dummy loop indices 222 222 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 223 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwy, zfu, zfc, zfd ! 3D workspace223 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwy, zfu, zfc, zfd ! 3D workspace 224 224 !---------------------------------------------------------------------- 225 225 ! … … 233 233 ! 234 234 !--- Computation of the ustream and downstream value of the tracer and the mask 235 DO_2D( 0, 0, 0, 0 )235 DO_2D( nn_hls-1, nn_hls-1, 0, 0 ) 236 236 ! Upstream in the x-direction for the tracer 237 237 zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) … … 240 240 END_2D 241 241 END DO 242 CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 243 242 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 244 243 245 244 ! … … 247 246 ! --------------------------- 248 247 ! 249 DO_3D( 0, 0, 0, 0, 1, jpkm1 )248 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 250 249 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 251 250 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 252 251 END_3D 253 252 ! 254 DO_3D( 0, 0, 0, 0, 1, jpkm1 )253 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 255 254 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 256 255 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) … … 261 260 262 261 !--- Lateral boundary conditions 263 CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )262 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp ) 264 263 265 264 !--- QUICKEST scheme … … 267 266 ! 268 267 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 269 DO_3D( 0, 0, 0, 0, 1, jpkm1 )268 DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 ) 270 269 zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 271 270 END_3D 272 CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) !--- Lateral boundary conditions271 IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) !--- Lateral boundary conditions 273 272 ! 274 273 ! Tracer flux on the x-direction 275 DO jk = 1, jpkm1 276 ! 277 DO_2D( 0, 0, 0, 0 ) 278 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 279 !--- If the second ustream point is a land point 280 !--- the flux is computed by the 1st order UPWIND scheme 281 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 282 zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 283 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk) 284 END_2D 285 END DO 286 ! 287 CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 274 DO_3D( 1, 0, 0, 0, 1, jpkm1 ) 275 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 276 !--- If the second ustream point is a land point 277 !--- the flux is computed by the 1st order UPWIND scheme 278 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 279 zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 280 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk) 281 END_3D 288 282 ! 289 283 ! Computation of the trend … … 313 307 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 314 308 INTEGER , INTENT(in ) :: kjpt ! number of tracers 315 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pW ! vertical velocity 309 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 310 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pW ! vertical velocity 316 311 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation 317 312 ! 318 313 INTEGER :: ji, jj, jk, jn ! dummy loop indices 319 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwz ! 3D workspace314 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwz ! 3D workspace 320 315 !!---------------------------------------------------------------------- 321 316 ! … … 332 327 IF( ln_linssh ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask) 333 328 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 334 DO_2D( 1, 1, 1, 1)329 DO_2D( 0, 0, 0, 0 ) 335 330 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) ! linear free surface 336 331 END_2D 337 332 ELSE ! no ocean cavities (only ocean surface) 338 zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm) 333 DO_2D( 0, 0, 0, 0 ) 334 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm) 335 END_2D 339 336 ENDIF 340 337 ENDIF … … 359 356 !! ** Method : 360 357 !!---------------------------------------------------------------------- 361 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point362 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point363 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)364 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux358 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfu ! second upwind point 359 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfd ! first douwning point 360 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point) 361 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux 365 362 !! 366 363 INTEGER :: ji, jj, jk ! dummy loop indices … … 369 366 !---------------------------------------------------------------------- 370 367 ! 371 DO_3D( 1, 1, 1, 1, 1, jpkm1 )368 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 372 369 zc = puc(ji,jj,jk) ! Courant number 373 370 zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) -
NEMO/trunk/src/OCE/TRA/traadv_ubs.F90
r13497 r13982 92 92 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 93 93 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 94 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 94 95 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume transport components 95 96 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 99 100 REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - - 100 101 REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - - 101 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu, ztv, zltu, zltv, zti, ztw ! 3D workspace 102 !!---------------------------------------------------------------------- 103 ! 104 IF( kt == kit000 ) THEN 105 IF(lwp) WRITE(numout,*) 106 IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype 107 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 102 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztu, ztv, zltu, zltv, zti, ztw ! 3D workspace 103 !!---------------------------------------------------------------------- 104 ! 105 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 106 IF( kt == kit000 ) THEN 107 IF(lwp) WRITE(numout,*) 108 IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype 109 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 110 ENDIF 111 ! 112 l_trd = .FALSE. 113 l_hst = .FALSE. 114 l_ptr = .FALSE. 115 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 116 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 117 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 118 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 108 119 ENDIF 109 !110 l_trd = .FALSE.111 l_hst = .FALSE.112 l_ptr = .FALSE.113 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.114 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.115 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &116 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.117 120 ! 118 121 ztw (:,:, 1 ) = 0._wp ! surface & bottom value : set to zero for all tracers 119 122 zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp 120 123 ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp 121 !122 124 ! ! =========== 123 125 DO jn = 1, kjpt ! tracer loop … … 125 127 ! 126 128 DO jk = 1, jpkm1 !== horizontal laplacian of before tracer ==! 127 DO_2D( 1, 0, 1, 0) ! First derivative (masked gradient)129 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! First derivative (masked gradient) 128 130 zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk) 129 131 zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk) … … 131 133 ztv(ji,jj,jk) = zeev * ( pt(ji ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 132 134 END_2D 133 DO_2D( 0, 0, 0, 0) ! Second derivative (divergence)135 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! Second derivative (divergence) 134 136 zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) ) 135 137 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef … … 138 140 ! 139 141 END DO 140 CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp ) ; CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn)142 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 141 143 ! 142 144 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== Horizontal advective fluxes ==! (UBS) … … 153 155 END_3D 154 156 ! 155 zltu(:,:,:) = pt(:,:,:,jn,Krhs) ! store the initial trends before its update 157 DO_3D( 1, 1, 1, 1, 1, jpk ) 158 zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) ! store the initial trends before its update 159 END_3D 156 160 ! 157 161 DO jk = 1, jpkm1 !== add the horizontal advective trend ==! … … 165 169 END DO 166 170 ! 167 zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:) ! Horizontal advective trend used in vertical 2nd order FCT case 168 ! ! and/or in trend diagnostic (l_trd=T) 169 ! 171 DO_3D( 1, 1, 1, 1, 1, jpk ) 172 zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltu(ji,jj,jk) ! Horizontal advective trend used in vertical 2nd order FCT case 173 END_3D ! and/or in trend diagnostic (l_trd=T) 174 ! 170 175 IF( l_trd ) THEN ! trend diagnostics 171 176 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztu, pU, pt(:,:,:,jn,Kmm) ) … … 185 190 CASE( 2 ) ! 2nd order FCT 186 191 ! 187 IF( l_trd ) zltv(:,:,:) = pt(:,:,:,jn,Krhs) ! store pt(:,:,:,:,Krhs) if trend diag. 192 IF( l_trd ) THEN 193 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 194 zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) ! store pt(:,:,:,:,Krhs) if trend diag. 195 END_3D 196 ENDIF 188 197 ! 189 198 ! !* upstream advection with initial mass fluxes & intermediate update ==! … … 199 208 END_2D 200 209 ELSE ! no cavities: only at the ocean surface 201 ztw(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb) 210 DO_2D( 1, 1, 1, 1 ) 211 ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 212 END_2D 202 213 ENDIF 203 214 ENDIF … … 209 220 zti(ji,jj,jk) = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 210 221 END_3D 211 CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1.0_wp ) ! Lateral boundary conditions on zti, zsi (unchanged sign)212 222 ! 213 223 ! !* anti-diffusive flux : high order minus low order … … 226 236 ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 227 237 END_3D 228 IF( ln_linssh ) ztw(:,:, 1 ) = pW(:,:,1) * pt(:,:,1,jn,Kmm) !!gm ISF & 4th COMPACT doesn't work 238 IF( ln_linssh ) THEN 239 DO_2D( 1, 1, 1, 1 ) 240 ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm) !!gm ISF & 4th COMPACT doesn't work 241 END_2D 242 ENDIF 229 243 ! 230 244 END SELECT … … 262 276 !! in-space based differencing for fluid 263 277 !!---------------------------------------------------------------------- 264 INTEGER , INTENT(in ) 265 REAL(wp), INTENT(in ) 266 REAL(wp), DIMENSION 267 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field268 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pcc ! monotonic flux in the k direction278 INTEGER , INTENT(in ) :: Kmm ! time level index 279 REAL(wp), INTENT(in ) :: p2dt ! tracer time-step 280 REAL(wp), DIMENSION(jpi,jpj,jpk) :: pbef ! before field 281 REAL(wp), INTENT(inout), DIMENSION(A2D(nn_hls) ,jpk) :: paft ! after field 282 REAL(wp), INTENT(inout), DIMENSION(A2D(nn_hls) ,jpk) :: pcc ! monotonic flux in the k direction 269 283 ! 270 284 INTEGER :: ji, jj, jk ! dummy loop indices 271 285 INTEGER :: ikm1 ! local integer 272 286 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 273 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zbetup, zbetdo! 3D workspace287 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo ! 3D workspace 274 288 !!---------------------------------------------------------------------- 275 289 ! … … 281 295 ! -------------------- 282 296 ! ! large negative value (-zbig) inside land 283 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 284 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 297 DO_3D( 0, 0, 0, 0, 1, jpk ) 298 pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) ) 299 paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) ) 300 END_3D 285 301 ! 286 302 DO jk = 1, jpkm1 ! search maximum in neighbourhood … … 293 309 END DO 294 310 ! ! large positive value (+zbig) inside land 295 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 296 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 311 DO_3D( 0, 0, 0, 0, 1, jpk ) 312 pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) ) 313 paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) ) 314 END_3D 297 315 ! 298 316 DO jk = 1, jpkm1 ! search minimum in neighbourhood … … 305 323 END DO 306 324 ! ! restore masked values to zero 307 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 308 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 325 DO_3D( 0, 0, 0, 0, 1, jpk ) 326 pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) 327 paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) 328 END_3D 309 329 ! 310 330 ! Positive and negative part of fluxes and beta terms -
NEMO/trunk/src/OCE/TRA/traatf.F90
r13295 r13982 156 156 ENDIF 157 157 ! 158 CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kbb) , 'T', 1.0_wp, & 159 & pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp, & 160 & pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp ) 161 ! 158 CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kmm) , 'T', 1.0_wp, pts(:,:,:,jp_sal,Kmm) , 'T', 1.0_wp ) 159 162 160 ENDIF 163 161 ! … … 210 208 DO jn = 1, kjpt 211 209 ! 212 DO_3D( 0, 0, 0, 0, 1, jpkm1 )210 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 213 211 ztn = pt(ji,jj,jk,jn,Kmm) 214 212 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers … … 275 273 zfact2 = zfact1 * r1_rho0 276 274 DO jn = 1, kjpt 277 DO_3D( 0, 0, 0, 0, 1, jpkm1 )275 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 278 276 ze3t_b = e3t(ji,jj,jk,Kbb) 279 277 ze3t_n = e3t(ji,jj,jk,Kmm) -
NEMO/trunk/src/OCE/TRA/traatf_qco.F90
r13295 r13982 149 149 ENDIF 150 150 ! 151 CALL lbc_lnk_multi( 'traatfqco', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., & 152 & pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., & 153 & pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 154 ! 151 CALL lbc_lnk_multi( 'traatfqco', pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1. ) 152 155 153 ENDIF 156 154 ! … … 203 201 DO jn = 1, kjpt 204 202 ! 205 DO_3D( 0, 0, 0, 0, 1, jpkm1 )203 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 206 204 ztn = pt(ji,jj,jk,jn,Kmm) 207 205 ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb) ! time laplacian on tracers … … 268 266 zfact2 = zfact1 * r1_rho0 269 267 DO jn = 1, kjpt 270 DO_3D( 0, 0, 0, 0, 1, jpkm1 )268 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 271 269 ze3t_b = e3t(ji,jj,jk,Kbb) 272 270 ze3t_n = e3t(ji,jj,jk,Kmm) -
NEMO/trunk/src/OCE/TRA/trabbc.F90
r13295 r13982 80 80 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 81 81 ! 82 INTEGER :: ji, jj ! dummy loop indices82 INTEGER :: ji, jj, jk ! dummy loop indices 83 83 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt ! 3D workspace 84 84 !!---------------------------------------------------------------------- … … 86 86 IF( ln_timing ) CALL timing_start('tra_bbc') 87 87 ! 88 IF( l_trdtra ) THEN! Save the input temperature trend88 IF( l_trdtra ) THEN ! Save the input temperature trend 89 89 ALLOCATE( ztrdt(jpi,jpj,jpk) ) 90 90 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) … … 96 96 END_2D 97 97 ! 98 CALL lbc_lnk( 'trabbc', pts(:,:,:,jp_tem,Krhs) , 'T', 1.0_wp )99 !100 98 IF( l_trdtra ) THEN ! Send the trend for diagnostics 101 99 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) … … 104 102 ENDIF 105 103 ! 106 CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 104 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 105 CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 106 ENDIF 107 107 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 108 108 ! -
NEMO/trunk/src/OCE/TRA/trabbl.F90
r13532 r13982 106 106 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 107 107 ! 108 INTEGER :: ji, jj, jk ! Dummy loop indices 108 109 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds 109 110 !!---------------------------------------------------------------------- … … 112 113 ! 113 114 IF( l_trdtra ) THEN !* Save the T-S input trends 114 ALLOCATE( ztrdt(jpi,jpj,jpk) 115 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 115 116 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 116 117 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) … … 125 126 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf - Ta: ', mask1=tmask, & 126 127 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 127 ! lateral boundary conditions ; just need for outputs128 CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1.0_wp , ahv_bbl, 'V', 1.0_wp )129 CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef130 CALL iom_put( "ahv_bbl", ahv_bbl ) ! bbl diffusive flux j-coef128 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 129 CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef 130 CALL iom_put( "ahv_bbl", ahv_bbl ) ! bbl diffusive flux j-coef 131 ENDIF 131 132 ! 132 133 ENDIF … … 136 137 CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm ) 137 138 IF(sn_cfctl%l_prtctl) & 138 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv - Ta: ', mask1=tmask, 139 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv - Ta: ', mask1=tmask, & 139 140 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 140 ! lateral boundary conditions ; just need for outputs 141 CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp ) 142 CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport 143 CALL iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport 141 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 142 ! lateral boundary conditions ; just need for outputs 143 ! NOTE: [tiling-comms-merge] The diagnostic results change along the north fold if this is removed 144 CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp ) 145 CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport 146 CALL iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport 147 ENDIF 144 148 ! 145 149 ENDIF … … 187 191 INTEGER :: ik ! local integers 188 192 REAL(wp) :: zbtr ! local scalars 189 REAL(wp), DIMENSION( jpi,jpj) :: zptb ! workspace193 REAL(wp), DIMENSION(A2D(nn_hls)) :: zptb ! workspace 190 194 !!---------------------------------------------------------------------- 191 195 ! … … 235 239 INTEGER :: iis , iid , ijs , ijd ! local integers 236 240 INTEGER :: ikus, ikud, ikvs, ikvd ! - - 241 INTEGER :: isi, isj ! - - 237 242 REAL(wp) :: zbtr, ztra ! local scalars 238 243 REAL(wp) :: zu_bbl, zv_bbl ! - - 239 244 !!---------------------------------------------------------------------- 240 245 ! 246 IF( ntsi == Nis0 ) THEN ; isi = 1 ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 247 IF( ntsj == Njs0 ) THEN ; isj = 1 ; ELSE ; isj = 0 ; ENDIF 241 248 ! ! =========== 242 249 DO jn = 1, kjpt ! tracer loop 243 250 ! ! =========== 244 DO jj = 1, jpjm1 245 DO ji = 1, jpim1 ! CAUTION start from i=1 to update i=2 when cyclic east-west 246 IF( utr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection 247 ! down-slope i/k-indices (deep) & up-slope i/k indices (shelf) 248 iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 249 ikud = mbku_d(ji,jj) ; ikus = mbku(ji,jj) 250 zu_bbl = ABS( utr_bbl(ji,jj) ) 251 ! 252 ! ! up -slope T-point (shelf bottom point) 253 zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm) 254 ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr 255 pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra 256 ! 257 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 258 zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm) 259 ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr 260 pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra 261 END DO 262 ! 263 zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm) 264 ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr 265 pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra 266 ENDIF 267 ! 268 IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero j-direction bbl advection 269 ! down-slope j/k-indices (deep) & up-slope j/k indices (shelf) 270 ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 271 ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj) 272 zv_bbl = ABS( vtr_bbl(ji,jj) ) 273 ! 274 ! up -slope T-point (shelf bottom point) 275 zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm) 276 ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr 277 pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra 278 ! 279 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 280 zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm) 281 ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr 282 pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn) + ztra 283 END DO 284 ! ! down-slope T-point (deep bottom point) 285 zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm) 286 ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr 287 pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra 288 ENDIF 289 END DO 251 ! NOTE: [tiling-comms-merge] Bug fix- correct order of indices 252 DO_2D( isj, 0, isi, 0 ) ! CAUTION start from i=1 to update i=2 when cyclic east-west 253 IF( utr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection 254 ! down-slope i/k-indices (deep) & up-slope i/k indices (shelf) 255 iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 256 ikud = mbku_d(ji,jj) ; ikus = mbku(ji,jj) 257 zu_bbl = ABS( utr_bbl(ji,jj) ) 258 ! 259 ! ! up -slope T-point (shelf bottom point) 260 zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm) 261 ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr 262 pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra 263 ! 264 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 265 zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm) 266 ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr 267 pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra 268 END DO 269 ! 270 zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm) 271 ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr 272 pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra 273 ENDIF 290 274 ! 291 END DO 292 ! ! =========== 293 END DO ! end tracer 294 ! ! =========== 275 IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero j-direction bbl advection 276 ! down-slope j/k-indices (deep) & up-slope j/k indices (shelf) 277 ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 278 ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj) 279 zv_bbl = ABS( vtr_bbl(ji,jj) ) 280 ! 281 ! up -slope T-point (shelf bottom point) 282 zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm) 283 ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr 284 pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra 285 ! 286 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 287 zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm) 288 ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr 289 pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn) + ztra 290 END DO 291 ! ! down-slope T-point (deep bottom point) 292 zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm) 293 ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr 294 pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra 295 ENDIF 296 END_2D 297 ! ! =========== 298 END DO ! end tracer 299 ! ! =========== 295 300 END SUBROUTINE tra_bbl_adv 296 301 … … 333 338 REAL(wp) :: za, zb, zgdrho ! local scalars 334 339 REAL(wp) :: zsign, zsigna, zgbbl ! - - 335 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts, zab ! 3D workspace 336 REAL(wp), DIMENSION(jpi,jpj) :: zub, zvb, zdep ! 2D workspace 337 !!---------------------------------------------------------------------- 338 ! 339 IF( kt == kit000 ) THEN 340 IF(lwp) WRITE(numout,*) 341 IF(lwp) WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype 342 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 340 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zts, zab ! 3D workspace 341 REAL(wp), DIMENSION(A2D(nn_hls)) :: zub, zvb, zdep ! 2D workspace 342 !!---------------------------------------------------------------------- 343 ! 344 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 345 IF( kt == kit000 ) THEN 346 IF(lwp) WRITE(numout,*) 347 IF(lwp) WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype 348 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 349 ENDIF 343 350 ENDIF 344 351 ! !* bottom variables (T, S, alpha, beta, depth, velocity) -
NEMO/trunk/src/OCE/TRA/tradmp.F90
r13295 r13982 95 95 ! 96 96 INTEGER :: ji, jj, jk, jn ! dummy loop indices 97 REAL(wp), DIMENSION( jpi,jpj,jpk,jpts) :: zts_dta97 REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta 98 98 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdts 99 99 !!---------------------------------------------------------------------- … … 102 102 ! 103 103 IF( l_trdtra ) THEN !* Save ta and sa trends 104 ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 105 ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) 104 ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 105 ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) 106 106 ENDIF 107 107 ! !== input T-S data at kt ==! … … 144 144 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) ) 145 145 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) ) 146 DEALLOCATE( ztrdts ) 146 DEALLOCATE( ztrdts ) 147 147 ENDIF 148 148 ! ! Control print -
NEMO/trunk/src/OCE/TRA/traisf.F90
r13295 r13982 11 11 !!---------------------------------------------------------------------- 12 12 USE isf_oce ! Ice shelf variables 13 USE par_oce , ONLY : nijtile, ntile, ntsi, ntei, ntsj, ntej 13 14 USE dom_oce ! ocean space domain variables 14 15 USE isfutils, ONLY : debug ! debug option … … 46 47 IF( ln_timing ) CALL timing_start('tra_isf') 47 48 ! 48 IF( kt == nit000 ) THEN 49 IF(lwp) WRITE(numout,*) 50 IF(lwp) WRITE(numout,*) 'tra_isf : Ice shelf heat fluxes' 51 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 49 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 50 IF( kt == nit000 ) THEN 51 IF(lwp) WRITE(numout,*) 52 IF(lwp) WRITE(numout,*) 'tra_isf : Ice shelf heat fluxes' 53 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 54 ENDIF 52 55 ENDIF 53 56 ! … … 76 79 ! 77 80 IF ( ln_isfdebug ) THEN 78 CALL debug('tra_isf: pts(:,:,:,:,Krhs) T', pts(:,:,:,1,Krhs)) 79 CALL debug('tra_isf: pts(:,:,:,:,Krhs) S', pts(:,:,:,2,Krhs)) 81 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 82 CALL debug('tra_isf: pts(:,:,:,:,Krhs) T', pts(:,:,:,1,Krhs)) 83 CALL debug('tra_isf: pts(:,:,:,:,Krhs) S', pts(:,:,:,2,Krhs)) 84 ENDIF 80 85 END IF 81 86 ! … … 101 106 INTEGER :: ji,jj,jk ! loop index 102 107 INTEGER :: ikt, ikb ! top and bottom level of the tbl 103 REAL(wp), DIMENSION( jpi,jpj):: ztc ! total ice shelf tracer trend108 REAL(wp), DIMENSION(A2D(nn_hls)) :: ztc ! total ice shelf tracer trend 104 109 !!---------------------------------------------------------------------- 105 110 ! 106 111 ! compute 2d total trend due to isf 107 ztc(:,:) = 0.5_wp * ( ptsc(:,:,jp_tem) + ptsc_b(:,:,jp_tem) ) / phtbl(:,:) 112 DO_2D( 0, 0, 0, 0 ) 113 ztc(ji,jj) = 0.5_wp * ( ptsc(ji,jj,jp_tem) + ptsc_b(ji,jj,jp_tem) ) / phtbl(ji,jj) 114 END_2D 108 115 ! 109 116 ! update pts(:,:,:,:,Krhs) 110 DO_2D( 1, 1, 1, 1)117 DO_2D( 0, 0, 0, 0 ) 111 118 ! 112 119 ikt = ktop(ji,jj) … … 137 144 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: ptsc 138 145 !!---------------------------------------------------------------------- 139 INTEGER :: j k146 INTEGER :: ji, jj, jk 140 147 !!---------------------------------------------------------------------- 141 148 ! 142 DO jk = 1,jpk 143 ptsa(:,:,jk,jp_tem) = & 144 & ptsa(:,:,jk,jp_tem) + ptsc(:,:,jk,jp_tem) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) 145 ptsa(:,:,jk,jp_sal) = & 146 & ptsa(:,:,jk,jp_sal) + ptsc(:,:,jk,jp_sal) * r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) 147 END DO 149 DO_3D( 0, 0, 0, 0, 1, jpk ) 150 ptsa(ji,jj,jk,jp_tem) = ptsa(ji,jj,jk,jp_tem) + ptsc(ji,jj,jk,jp_tem) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 151 ptsa(ji,jj,jk,jp_sal) = ptsa(ji,jj,jk,jp_sal) + ptsc(ji,jj,jk,jp_sal) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 152 END_3D 148 153 ! 149 154 END SUBROUTINE tra_isf_cpl -
NEMO/trunk/src/OCE/TRA/traldf.F90
r12377 r13982 17 17 USE oce ! ocean dynamics and tracers 18 18 USE dom_oce ! ocean space and time domain 19 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 20 USE domain, ONLY : dom_tile 19 21 USE phycst ! physical constants 20 22 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. … … 37 39 PUBLIC tra_ldf ! called by step.F90 38 40 PUBLIC tra_ldf_init ! called by nemogcm.F90 39 41 40 42 !!---------------------------------------------------------------------- 41 43 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 56 58 !! 57 59 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds 60 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 61 LOGICAL :: lskip 58 62 !!---------------------------------------------------------------------- 59 63 ! 60 64 IF( ln_timing ) CALL timing_start('tra_ldf') 61 65 ! 66 lskip = .FALSE. 67 62 68 IF( l_trdtra ) THEN !* Save ta and sa trends 63 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 64 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 69 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 70 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 65 71 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 66 72 ENDIF 67 ! 68 SELECT CASE ( nldf_tra ) !* compute lateral mixing trend and add it to the general trend 69 CASE ( np_lap ) ! laplacian: iso-level operator 70 CALL tra_ldf_lap ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 71 CASE ( np_lap_i ) ! laplacian: standard iso-neutral operator (Madec) 72 CALL tra_ldf_iso ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 73 CASE ( np_lap_it ) ! laplacian: triad iso-neutral operator (griffies) 74 CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 75 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: iso-level & iso-neutral operators 76 CALL tra_ldf_blp ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, nldf_tra ) 77 END SELECT 78 ! 79 IF( l_trdtra ) THEN !* save the horizontal diffusive trends for further diagnostics 80 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 81 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 82 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ldf, ztrdt ) 83 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds ) 84 DEALLOCATE( ztrdt, ztrds ) 73 74 ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 75 IF( nldf_tra == np_blp .OR. nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) THEN 76 IF( ln_tile ) THEN 77 IF( ntile == 1 ) THEN 78 CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 79 ELSE 80 lskip = .TRUE. 81 ENDIF 82 ENDIF 83 ENDIF 84 IF( .NOT. lskip ) THEN 85 ! 86 SELECT CASE ( nldf_tra ) !* compute lateral mixing trend and add it to the general trend 87 CASE ( np_lap ) ! laplacian: iso-level operator 88 CALL tra_ldf_lap ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 89 CASE ( np_lap_i ) ! laplacian: standard iso-neutral operator (Madec) 90 CALL tra_ldf_iso ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 91 CASE ( np_lap_it ) ! laplacian: triad iso-neutral operator (griffies) 92 CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 93 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: iso-level & iso-neutral operators 94 ! NOTE: [tiling-comms-merge] This lbc_lnk is still needed in the zco case, because zps_hde is not called in step 95 IF(nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.) 96 CALL tra_ldf_blp ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, nldf_tra ) 97 END SELECT 98 ! 99 IF( l_trdtra ) THEN !* save the horizontal diffusive trends for further diagnostics 100 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 101 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 102 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ldf, ztrdt ) 103 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds ) 104 DEALLOCATE( ztrdt, ztrds ) 105 ENDIF 106 107 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 108 IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) 85 109 ENDIF 86 110 ! !* print mean trends (used for debugging) 87 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' ldf - Ta: ', mask1=tmask, 111 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' ldf - Ta: ', mask1=tmask, & 88 112 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 89 113 ! -
NEMO/trunk/src/OCE/TRA/traldf_iso.F90
r13497 r13982 19 19 USE oce ! ocean dynamics and active tracers 20 20 USE dom_oce ! ocean space and time domain 21 USE domutl, ONLY : is_tile 21 22 USE trc_oce ! share passive tracers/Ocean variables 22 23 USE zdf_oce ! ocean vertical physics … … 49 50 CONTAINS 50 51 51 SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv, & 52 & pgu , pgv , pgui, pgvi, & 53 & pt , pt2 , pt_rhs , kjpt , kpass ) 52 SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv, & 53 & pgu , pgv , pgui, pgvi, & 54 & pt, pt2, pt_rhs, kjpt, kpass ) 55 !! 56 INTEGER , INTENT(in ) :: kt ! ocean time-step index 57 INTEGER , INTENT(in ) :: kit000 ! first time step index 58 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 59 INTEGER , INTENT(in ) :: kjpt ! number of tracers 60 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 61 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 62 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 63 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 64 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 65 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 66 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 67 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pt_rhs ! tracer trend 68 !! 69 CALL tra_ldf_iso_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu), & 70 & pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui), & 71 & pt, is_tile(pt), pt2, is_tile(pt2), pt_rhs, is_tile(pt_rhs), kjpt, kpass ) 72 END SUBROUTINE tra_ldf_iso 73 74 75 SUBROUTINE tra_ldf_iso_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah, & 76 & pgu , pgv , ktg , pgui, pgvi, ktgi, & 77 & pt, ktt, pt2, ktt2, pt_rhs, ktt_rhs, kjpt, kpass ) 54 78 !!---------------------------------------------------------------------- 55 79 !! *** ROUTINE tra_ldf_iso *** … … 92 116 !! ** Action : Update pt_rhs arrays with the before rotated diffusion 93 117 !!---------------------------------------------------------------------- 94 INTEGER , INTENT(in ) :: kt ! ocean time-step index 95 INTEGER , INTENT(in ) :: kit000 ! first time step index 96 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 97 INTEGER , INTENT(in ) :: kjpt ! number of tracers 98 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 99 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 100 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 101 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 102 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 103 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 104 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 105 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend 118 INTEGER , INTENT(in ) :: kt ! ocean time-step index 119 INTEGER , INTENT(in ) :: kit000 ! first time step index 120 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 121 INTEGER , INTENT(in ) :: kjpt ! number of tracers 122 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 123 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 124 INTEGER , INTENT(in ) :: ktah, ktg, ktgi, ktt, ktt2, ktt_rhs 125 REAL(wp), DIMENSION(A2D_T(ktah) ,JPK) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 126 REAL(wp), DIMENSION(A2D_T(ktg) ,KJPT), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 127 REAL(wp), DIMENSION(A2D_T(ktgi) ,KJPT), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 128 REAL(wp), DIMENSION(A2D_T(ktt) ,JPK,KJPT), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 129 REAL(wp), DIMENSION(A2D_T(ktt2) ,JPK,KJPT), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 130 REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend 106 131 ! 107 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 111 136 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - 112 137 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 113 REAL(wp), DIMENSION( jpi,jpj) :: zdkt, zdk1t, z2d114 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw138 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdkt, zdk1t, z2d 139 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdit, zdjt, zftu, zftv, ztfw 115 140 !!---------------------------------------------------------------------- 116 141 ! 117 142 IF( kpass == 1 .AND. kt == kit000 ) THEN 118 IF(lwp) WRITE(numout,*) 119 IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 120 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 121 ! 122 akz (:,:,:) = 0._wp 123 ah_wslp2(:,:,:) = 0._wp 143 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 144 IF(lwp) WRITE(numout,*) 145 IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 146 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 147 ENDIF 148 ! 149 DO_3D( 0, 0, 0, 0, 1, jpk ) 150 akz (ji,jj,jk) = 0._wp 151 ah_wslp2(ji,jj,jk) = 0._wp 152 END_3D 124 153 ENDIF 125 ! 126 l_hst = .FALSE. 127 l_ptr = .FALSE. 128 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE. 129 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 130 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 154 ! 155 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 156 l_hst = .FALSE. 157 l_ptr = .FALSE. 158 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE. 159 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 160 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 161 ENDIF 131 162 ! 132 163 ! … … 167 198 ! 168 199 IF( ln_traldf_blp ) THEN ! bilaplacian operator 169 DO_3D( 1, 0, 1, 0, 2, jpkm1 )200 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 170 201 akz(ji,jj,jk) = 16._wp & 171 202 & * ah_wslp2 (ji,jj,jk) & … … 175 206 END_3D 176 207 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 177 DO_3D( 1, 0, 1, 0, 2, jpkm1 )208 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 178 209 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 179 210 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) … … 183 214 ! 184 215 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 185 akz(:,:,:) = ah_wslp2(:,:,:) 216 DO_3D( 0, 0, 0, 0, 1, jpk ) 217 akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 218 END_3D 186 219 ENDIF 187 220 ENDIF … … 195 228 !!---------------------------------------------------------------------- 196 229 !!gm : bug.... why (x,:,:)? (1,jpj,:) and (jpi,1,:) should be sufficient.... 197 zdit ( 1,:,:) = 0._wp ; zdit (jpi,:,:) = 0._wp198 zdjt ( 1,:,:) = 0._wp ; zdjt (jpi,:,:) = 0._wp230 zdit (ntsi-nn_hls,:,:) = 0._wp ; zdit (ntei+nn_hls,:,:) = 0._wp 231 zdjt (ntsi-nn_hls,:,:) = 0._wp ; zdjt (ntei+nn_hls,:,:) = 0._wp 199 232 !!end 200 233 … … 223 256 DO jk = 1, jpkm1 ! Horizontal slab 224 257 ! 225 ! !== Vertical tracer gradient 226 zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1) ! level jk+1 227 ! 228 IF( jk == 1 ) THEN ; zdkt(:,:) = zdk1t(:,:) ! surface: zdkt(jk=1)=zdkt(jk=2) 229 ELSE ; zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk) 230 ENDIF 258 DO_2D( 1, 1, 1, 1 ) 259 ! !== Vertical tracer gradient 260 zdk1t(ji,jj) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) ! level jk+1 261 ! 262 IF( jk == 1 ) THEN ; zdkt(ji,jj) = zdk1t(ji,jj) ! surface: zdkt(jk=1)=zdkt(jk=2) 263 ELSE ; zdkt(ji,jj) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 264 ENDIF 265 END_2D 266 ! 231 267 DO_2D( 1, 0, 1, 0 ) !== Horizontal fluxes 232 268 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) … … 330 366 END DO ! end tracer loop 331 367 ! 332 END SUBROUTINE tra_ldf_iso 368 END SUBROUTINE tra_ldf_iso_t 333 369 334 370 !!============================================================================== -
NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90
r13497 r13982 13 13 USE oce ! ocean dynamics and active tracers 14 14 USE dom_oce ! ocean space and time domain 15 USE domutl, ONLY : is_tile 15 16 USE ldftra ! lateral physics: eddy diffusivity 16 17 USE traldf_iso ! iso-neutral lateral diffusion (standard operator) (tra_ldf_iso routine) … … 46 47 CONTAINS 47 48 48 SUBROUTINE tra_ldf_lap( kt, Kmm, kit000, cdtype, pahu, pahv , & 49 & pgu , pgv , pgui, pgvi, & 50 & pt , pt_rhs, kjpt, kpass ) 49 SUBROUTINE tra_ldf_lap( kt, Kmm, kit000, cdtype, pahu, pahv, & 50 & pgu , pgv , pgui, pgvi, & 51 & pt, pt_rhs, kjpt, kpass ) 52 !! 53 INTEGER , INTENT(in ) :: kt ! ocean time-step index 54 INTEGER , INTENT(in ) :: kit000 ! first time step index 55 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 56 INTEGER , INTENT(in ) :: kjpt ! number of tracers 57 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 58 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 59 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 60 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 61 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 62 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt ! before tracer fields 63 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pt_rhs ! tracer trend 64 !! 65 CALL tra_ldf_lap_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu), & 66 & pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui), & 67 & pt, is_tile(pt), pt_rhs, is_tile(pt_rhs), kjpt, kpass ) 68 END SUBROUTINE tra_ldf_lap 69 70 71 SUBROUTINE tra_ldf_lap_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah, & 72 & pgu , pgv , ktg , pgui, pgvi, ktgi, & 73 & pt, ktt, pt_rhs, ktt_rhs, kjpt, kpass ) 51 74 !!---------------------------------------------------------------------- 52 75 !! *** ROUTINE tra_ldf_lap *** … … 72 95 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 73 96 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 74 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 75 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 76 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 77 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before tracer fields 78 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend 79 ! 80 INTEGER :: ji, jj, jk, jn ! dummy loop indices 81 REAL(wp) :: zsign ! local scalars 82 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu, ztv, zaheeu, zaheev 83 !!---------------------------------------------------------------------- 84 ! 85 IF( kt == nit000 .AND. lwp ) THEN 86 WRITE(numout,*) 87 WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass 88 WRITE(numout,*) '~~~~~~~~~~~ ' 89 ENDIF 90 ! 91 l_hst = .FALSE. 92 l_ptr = .FALSE. 93 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE. 94 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 95 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 97 INTEGER , INTENT(in ) :: ktah, ktg, ktgi, ktt, ktt_rhs 98 REAL(wp), DIMENSION(A2D_T(ktah), JPK) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 99 REAL(wp), DIMENSION(A2D_T(ktg), KJPT), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 100 REAL(wp), DIMENSION(A2D_T(ktgi), KJPT), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 101 REAL(wp), DIMENSION(A2D_T(ktt), JPK,KJPT), INTENT(in ) :: pt ! before tracer fields 102 REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend 103 ! 104 INTEGER :: ji, jj, jk, jn ! dummy loop indices 105 INTEGER :: isi, iei, isj, iej ! local integers 106 REAL(wp) :: zsign ! local scalars 107 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztu, ztv, zaheeu, zaheev 108 !!---------------------------------------------------------------------- 109 ! 110 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 111 IF( kt == nit000 .AND. lwp ) THEN 112 WRITE(numout,*) 113 WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass 114 WRITE(numout,*) '~~~~~~~~~~~ ' 115 ENDIF 116 ! 117 l_hst = .FALSE. 118 l_ptr = .FALSE. 119 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) ) l_ptr = .TRUE. 120 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 121 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 122 ENDIF 96 123 ! 97 124 ! !== Initialization of metric arrays used for all tracers ==! … … 99 126 ELSE ; zsign = -1._wp 100 127 ENDIF 101 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 128 129 IF( ntsi == Nis0 ) THEN ; isi = nn_hls - 1 ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 130 IF( ntsj == Njs0 ) THEN ; isj = nn_hls - 1 ; ELSE ; isj = 0 ; ENDIF 131 IF( ntei == Nie0 ) THEN ; iei = nn_hls - 1 ; ELSE ; iei = 0 ; ENDIF 132 IF( ntej == Nje0 ) THEN ; iej = nn_hls - 1 ; ELSE ; iej = 0 ; ENDIF 133 134 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) !== First derivative (gradient) ==! 102 135 zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) !!gm * umask(ji,jj,jk) pah masked! 103 136 zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) !!gm * vmask(ji,jj,jk) … … 108 141 ! ! =========== ! 109 142 ! 110 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== First derivative (gradient) ==!143 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) !== First derivative (gradient) ==! 111 144 ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) 112 145 ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) 113 146 END_3D 114 147 IF( ln_zps ) THEN ! set gradient at bottom/top ocean level 115 DO_2D( 1, 0, 1, 0) ! bottom148 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! bottom 116 149 ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn) 117 150 ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn) 118 151 END_2D 119 152 IF( ln_isfcav ) THEN ! top in ocean cavities only 120 DO_2D( 1, 0, 1, 0)153 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 121 154 IF( miku(ji,jj) > 1 ) ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 122 155 IF( mikv(ji,jj) > 1 ) ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) … … 125 158 ENDIF 126 159 ! 127 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Second derivative (divergence) added to the general tracer trends ==! 160 ! NOTE: [tiling-comms-merge] Bounds changed to avoid repeating this calculation for overlapping rows when using tiling 161 DO_3D( isj, iej, isi, iei, 1, jpkm1 ) !== Second derivative (divergence) added to the general tracer trends ==! 128 162 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) & 129 163 & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & … … 142 176 ! ! ================== 143 177 ! 144 END SUBROUTINE tra_ldf_lap 178 END SUBROUTINE tra_ldf_lap_t 145 179 146 180 … … 173 207 ! 174 208 INTEGER :: ji, jj, jk, jn ! dummy loop indices 175 REAL(wp), DIMENSION( jpi,jpj,jpk,kjpt) :: zlap ! laplacian at t-point176 REAL(wp), DIMENSION( jpi,jpj, kjpt) :: zglu, zglv ! bottom GRADh of the laplacian (u- and v-points)177 REAL(wp), DIMENSION( jpi,jpj, kjpt) :: zgui, zgvi ! top GRADh of the laplacian (u- and v-points)209 REAL(wp), DIMENSION(A2D(nn_hls),jpk,kjpt) :: zlap ! laplacian at t-point 210 REAL(wp), DIMENSION(A2D(nn_hls), kjpt) :: zglu, zglv ! bottom GRADh of the laplacian (u- and v-points) 211 REAL(wp), DIMENSION(A2D(nn_hls), kjpt) :: zgui, zgvi ! top GRADh of the laplacian (u- and v-points) 178 212 !!--------------------------------------------------------------------- 179 213 ! 180 IF( kt == kit000 .AND. lwp ) THEN 181 WRITE(numout,*) 182 SELECT CASE ( kldf ) 183 CASE ( np_blp ) ; WRITE(numout,*) 'tra_ldf_blp : iso-level bilaplacian operator on ', cdtype 184 CASE ( np_blp_i ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)' 185 CASE ( np_blp_it ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)' 186 END SELECT 187 WRITE(numout,*) '~~~~~~~~~~~' 214 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 215 IF( kt == kit000 .AND. lwp ) THEN 216 WRITE(numout,*) 217 SELECT CASE ( kldf ) 218 CASE ( np_blp ) ; WRITE(numout,*) 'tra_ldf_blp : iso-level bilaplacian operator on ', cdtype 219 CASE ( np_blp_i ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)' 220 CASE ( np_blp_it ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)' 221 END SELECT 222 WRITE(numout,*) '~~~~~~~~~~~' 223 ENDIF 188 224 ENDIF 189 225 … … 200 236 END SELECT 201 237 ! 238 ! NOTE: [tiling-comms-merge] Needed for both nn_hls as tra_ldf_iso and tra_ldf_triad have not yet been adjusted to work with nn_hls = 2. In the zps case the lbc_lnk in zps_hde handles this, but in the zco case zlap always needs this lbc_lnk. I did try adjusting the bounds in tra_ldf_iso and tra_ldf_triad so this lbc_lnk was only needed for nn_hls = 1, but this was not correct and I did not have time to figure out why 202 239 CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp ) ! Lateral boundary conditions (unchanged sign) 203 240 ! ! Partial top/bottom cell: GRADh( zlap ) -
NEMO/trunk/src/OCE/TRA/traldf_triad.F90
r13497 r13982 13 13 USE oce ! ocean dynamics and active tracers 14 14 USE dom_oce ! ocean space and time domain 15 ! TEMP: [tiling] This change not necessary if XIOS has subdomain support 16 USE domain, ONLY : dom_tile 17 USE domutl, ONLY : is_tile 15 18 USE phycst ! physical constants 16 19 USE trc_oce ! share passive tracers/Ocean variables … … 33 36 PUBLIC tra_ldf_triad ! routine called by traldf.F90 34 37 35 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels36 37 38 LOGICAL :: l_ptr ! flag to compute poleward transport 38 39 LOGICAL :: l_hst ! flag to compute heat transport … … 49 50 CONTAINS 50 51 51 SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, & 52 & pgu , pgv , pgui, pgvi , & 53 & pt , pt2, pt_rhs, kjpt, kpass ) 52 SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, & 53 & pgu , pgv , pgui, pgvi, & 54 & pt, pt2, pt_rhs, kjpt, kpass ) 55 !! 56 INTEGER , INTENT(in ) :: kt ! ocean time-step index 57 INTEGER , INTENT(in ) :: kit000 ! first time step index 58 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 59 INTEGER , INTENT(in ) :: kjpt ! number of tracers 60 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 61 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 62 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 63 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 64 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 65 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 66 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 67 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pt_rhs ! tracer trend 68 !! 69 CALL tra_ldf_triad_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu), & 70 & pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui), & 71 & pt, is_tile(pt), pt2, is_tile(pt2), pt_rhs, is_tile(pt_rhs), kjpt, kpass ) 72 END SUBROUTINE tra_ldf_triad 73 74 75 SUBROUTINE tra_ldf_triad_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah, & 76 & pgu , pgv , ktg , pgui, pgvi, ktgi, & 77 & pt, ktt, pt2, ktt2, pt_rhs, ktt_rhs, kjpt, kpass ) 54 78 !!---------------------------------------------------------------------- 55 79 !! *** ROUTINE tra_ldf_triad *** … … 77 101 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 78 102 INTEGER , INTENT(in) :: Kmm ! ocean time level indices 79 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 80 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 81 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 84 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend 103 INTEGER , INTENT(in ) :: ktah, ktg, ktgi, ktt, ktt2, ktt_rhs 104 REAL(wp), DIMENSION(A2D_T(ktah), JPK) , INTENT(in ) :: pahu, pahv ! eddy diffusivity at u- and v-points [m2/s] 105 REAL(wp), DIMENSION(A2D_T(ktg), KJPT), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 106 REAL(wp), DIMENSION(A2D_T(ktgi), KJPT), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 107 REAL(wp), DIMENSION(A2D_T(ktt), JPK,KJPT), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2) 108 REAL(wp), DIMENSION(A2D_T(ktt2), JPK,KJPT), INTENT(in ) :: pt2 ! tracer (only used in kpass=2) 109 REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend 85 110 ! 86 111 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 94 119 REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 95 120 REAL(wp) :: zah, zah_slp, zaei_slp 96 REAL(wp), DIMENSION(jpi,jpj ) :: z2d ! 2D workspace 97 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw ! 3D - 121 REAL(wp), DIMENSION(A2D(nn_hls),0:1) :: zdkt3d ! vertical tracer gradient at 2 levels 122 REAL(wp), DIMENSION(A2D(nn_hls) ) :: z2d ! 2D workspace 123 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk) :: zdit, zdjt, zftu, zftv, ztfw ! 3D - 124 ! TEMP: [tiling] This can be A2D(nn_hls) if XIOS has subdomain support 125 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 98 126 !!---------------------------------------------------------------------- 99 127 ! 100 IF( .NOT.ALLOCATED(zdkt3d) ) THEN 101 ALLOCATE( zdkt3d(jpi,jpj,0:1) , STAT=ierr ) 102 CALL mpp_sum ( 'traldf_triad', ierr ) 103 IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_triad: unable to allocate arrays') 104 ENDIF 105 ! 106 IF( kpass == 1 .AND. kt == kit000 ) THEN 107 IF(lwp) WRITE(numout,*) 108 IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype 109 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 110 ENDIF 111 ! 112 l_hst = .FALSE. 113 l_ptr = .FALSE. 114 IF( cdtype == 'TRA' ) THEN 115 IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') ) l_ptr = .TRUE. 116 IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 117 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) l_hst = .TRUE. 128 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 129 IF( kpass == 1 .AND. kt == kit000 ) THEN 130 IF(lwp) WRITE(numout,*) 131 IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype 132 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 133 ENDIF 134 ! 135 l_hst = .FALSE. 136 l_ptr = .FALSE. 137 IF( cdtype == 'TRA' ) THEN 138 IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') ) l_ptr = .TRUE. 139 IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 140 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) l_hst = .TRUE. 141 ENDIF 118 142 ENDIF 119 143 ! … … 128 152 IF( kpass == 1 ) THEN !== first pass only and whatever the tracer is ==! 129 153 ! 130 akz (:,:,:) = 0._wp 131 ah_wslp2(:,:,:) = 0._wp 132 IF( ln_ldfeiv_dia ) THEN 133 zpsi_uw(:,:,:) = 0._wp 134 zpsi_vw(:,:,:) = 0._wp 135 ENDIF 154 DO_3D( 0, 0, 0, 0, 1, jpk ) 155 akz (ji,jj,jk) = 0._wp 156 ah_wslp2(ji,jj,jk) = 0._wp 157 END_3D 136 158 ! 137 159 DO ip = 0, 1 ! i-k triads 138 160 DO kp = 0, 1 139 DO_3D( 1, 0, 1, 0, 1, jpkm1 )140 ze3wr = 1._wp / e3w(ji +ip,jj,jk+kp,Kmm)141 zbu = e1e2u(ji ,jj) * e3u(ji,jj,jk,Kmm)142 zah = 0.25_wp * pahu(ji ,jj,jk)143 zslope_skew = triadi_g(ji +ip,jj,jk,1-ip,kp)161 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 162 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 163 zbu = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm) 164 zah = 0.25_wp * pahu(ji-ip,jj,jk) 165 zslope_skew = triadi_g(ji,jj,jk,1-ip,kp) 144 166 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 145 zslope2 = zslope_skew + ( gdept(ji +1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)167 zslope2 = zslope_skew + ( gdept(ji-ip+1,jj,jk,Kmm) - gdept(ji-ip,jj,jk,Kmm) ) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 146 168 zslope2 = zslope2 *zslope2 147 ah_wslp2(ji +ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2148 akz (ji +ip,jj,jk+kp) = akz (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj) &149 & * r1_e1u(ji ,jj) * umask(ji,jj,jk+kp)169 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 170 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + zah * r1_e1u(ji-ip,jj) & 171 & * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 150 172 ! 151 IF( ln_ldfeiv_dia ) zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &152 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew153 173 END_3D 154 174 END DO … … 157 177 DO jp = 0, 1 ! j-k triads 158 178 DO kp = 0, 1 159 DO_3D( 1, 0, 1, 0, 1, jpkm1 )160 ze3wr = 1.0_wp / e3w(ji,jj +jp,jk+kp,Kmm)161 zbv = e1e2v(ji,jj ) * e3v(ji,jj,jk,Kmm)162 zah = 0.25_wp * pahv(ji,jj ,jk)163 zslope_skew = triadj_g(ji,jj +jp,jk,1-jp,kp)179 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 180 ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 181 zbv = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm) 182 zah = 0.25_wp * pahv(ji,jj-jp,jk) 183 zslope_skew = triadj_g(ji,jj,jk,1-jp,kp) 164 184 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 165 185 ! (do this by *adding* gradient of depth) 166 zslope2 = zslope_skew + ( gdept(ji,jj +1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)186 zslope2 = zslope_skew + ( gdept(ji,jj-jp+1,jk,Kmm) - gdept(ji,jj-jp,jk,Kmm) ) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 167 187 zslope2 = zslope2 * zslope2 168 ah_wslp2(ji,jj +jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2169 akz (ji,jj +jp,jk+kp) = akz (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj) &170 & * r1_e2v(ji,jj ) * vmask(ji,jj,jk+kp)188 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 189 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + zah * r1_e2v(ji,jj-jp) & 190 & * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 171 191 ! 172 IF( ln_ldfeiv_dia ) zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) &173 & + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew174 192 END_3D 175 193 END DO … … 179 197 ! 180 198 IF( ln_traldf_blp ) THEN ! bilaplacian operator 181 DO_3D( 1, 0, 1, 0, 2, jpkm1 )199 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 182 200 akz(ji,jj,jk) = 16._wp & 183 201 & * ah_wslp2 (ji,jj,jk) & … … 187 205 END_3D 188 206 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 189 DO_3D( 1, 0, 1, 0, 2, jpkm1 )207 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 190 208 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 191 209 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) … … 195 213 ! 196 214 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 197 akz(:,:,:) = ah_wslp2(:,:,:) 198 ENDIF 199 ! 200 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 215 DO_3D( 0, 0, 0, 0, 1, jpk ) 216 akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 217 END_3D 218 ENDIF 219 ! 220 ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support 221 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 222 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 223 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 224 225 zpsi_uw(:,:,:) = 0._wp 226 zpsi_vw(:,:,:) = 0._wp 227 228 DO jp = 0, 1 229 DO kp = 0, 1 230 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 231 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 232 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp) 233 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 234 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp) 235 END_3D 236 END DO 237 END DO 238 CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 239 240 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile ) 241 ENDIF 242 ENDIF 201 243 ! 202 244 ENDIF !== end 1st pass only ==! … … 234 276 DO jk = 1, jpkm1 235 277 ! !== Vertical tracer gradient at level jk and jk+1 236 zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 278 DO_2D( 1, 1, 1, 1 ) 279 zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 280 END_2D 237 281 ! 238 282 ! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 239 283 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 240 ELSE ; zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk) 284 ELSE 285 DO_2D( 1, 1, 1, 1 ) 286 zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 287 END_2D 241 288 ENDIF 242 289 ! … … 380 427 END DO ! end tracer loop 381 428 ! ! =============== 382 END SUBROUTINE tra_ldf_triad 429 END SUBROUTINE tra_ldf_triad_t 383 430 384 431 !!============================================================================== -
NEMO/trunk/src/OCE/TRA/tramle.F90
r13497 r13982 79 79 !! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 80 80 !!---------------------------------------------------------------------- 81 INTEGER 82 INTEGER 83 INTEGER 84 CHARACTER(len=3) 85 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components86 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components87 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport81 INTEGER , INTENT(in ) :: kt ! ocean time-step index 82 INTEGER , INTENT(in ) :: kit000 ! first time step index 83 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 84 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 85 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 86 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pv ! out: same 3 transport components 87 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 88 88 ! 89 89 INTEGER :: ji, jj, jk ! dummy loop indices … … 91 91 REAL(wp) :: zcuw, zmuw, zc ! local scalar 92 92 REAL(wp) :: zcvw, zmvw ! - - 93 INTEGER , DIMENSION(jpi,jpj) :: inml_mle 94 REAL(wp), DIMENSION(jpi,jpj) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 95 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 93 INTEGER , DIMENSION(A2D(nn_hls)) :: inml_mle 94 REAL(wp), DIMENSION(A2D(nn_hls)) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_MH 95 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpsi_uw, zpsi_vw 96 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 97 REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: zLf_NH 98 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zpsiu_mle, zpsiv_mle 96 99 !!---------------------------------------------------------------------- 97 100 ! 98 101 ! !== MLD used for MLE ==! 99 102 ! ! compute from the 10m density to deal with the diurnal cycle 100 inml_mle(:,:) = mbkt(:,:) + 1 ! init. to number of ocean w-level (T-level + 1) 103 DO_2D( 1, 1, 1, 1 ) 104 inml_mle(ji,jj) = mbkt(ji,jj) + 1 ! init. to number of ocean w-level (T-level + 1) 105 END_2D 101 106 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 102 107 DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 (10m) … … 135 140 END SELECT 136 141 ! ! convert density into buoyancy 137 zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 142 DO_2D( 1, 1, 1, 1 ) 143 zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) ) 144 END_2D 138 145 ! 139 146 ! … … 206 213 END DO 207 214 215 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 208 216 IF( cdtype == 'TRA') THEN !== outputs ==! 209 ! 210 zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:) ! Lf = N H / f 211 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 217 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 218 ALLOCATE( zLf_NH(jpi,jpj), zpsiu_mle(jpi,jpj,jpk), zpsiv_mle(jpi,jpj,jpk) ) 219 zpsiu_mle(:,:,:) = 0._wp ; zpsiv_mle(:,:,:) = 0._wp 220 ENDIF 221 ! 222 DO_2D( 0, 0, 0, 0 ) 223 zLf_NH(ji,jj) = SQRT( rb_c * zmld(ji,jj) ) * r1_ft(ji,jj) ! Lf = N H / f 224 END_2D 212 225 ! 213 226 ! divide by cross distance to give streamfunction with dimensions m^2/s 214 DO jk = 1, ikmax+1 215 zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 216 zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 217 END DO 218 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 219 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 227 DO_3D( 0, 0, 0, 0, 1, ikmax+1 ) 228 zpsiu_mle(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj) 229 zpsiv_mle(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj) 230 END_3D 231 232 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 233 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 234 CALL iom_put( "psiu_mle", zpsiu_mle ) ! i-mle streamfunction 235 CALL iom_put( "psiv_mle", zpsiv_mle ) ! j-mle streamfunction 236 DEALLOCATE( zLf_NH, zpsiu_mle, zpsiv_mle ) 237 ENDIF 220 238 ENDIF 221 239 ! … … 283 301 IF( ierr /= 0 ) CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' ) 284 302 z1_t2 = 1._wp / ( rn_time * rn_time ) 285 DO_2D( 0, 1, 0, 1) ! "coriolis+ time^-1" at u- & v-points303 DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! "coriolis+ time^-1" at u- & v-points 286 304 zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp 287 305 zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp … … 289 307 rfv(ji,jj) = SQRT( zfv * zfv + z1_t2 ) 290 308 END_2D 291 CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1.0_wp , rfv, 'V', 1.0_wp )309 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1.0_wp , rfv, 'V', 1.0_wp ) 292 310 ! 293 311 ELSEIF( nn_mle == 1 ) THEN ! MLE array allocation & initialisation -
NEMO/trunk/src/OCE/TRA/tranpc.F90
r13497 r13982 17 17 USE oce ! ocean dynamics and active tracers 18 18 USE dom_oce ! ocean space and time domain 19 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed) 20 USE domain, ONLY : dom_tile 19 21 USE phycst ! physical constants 20 22 USE zdf_oce ! ocean vertical physics … … 32 34 33 35 PUBLIC tra_npc ! routine called by step.F90 36 37 INTEGER :: nnpcc ! number of statically instable water column 34 38 35 39 !! * Substitutions … … 64 68 ! 65 69 INTEGER :: ji, jj, jk ! dummy loop indices 66 INTEGER :: inpcc ! number of statically instable water column67 70 INTEGER :: jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low ! local integers 68 71 LOGICAL :: l_bottom_reached, l_column_treated … … 70 73 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt 71 74 REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp ! acceptance criteria for neutrality (N2==0) 72 REAL(wp), DIMENSION( jpk ) :: zvn2! vertical profile of N2 at 1 given point...73 REAL(wp), DIMENSION( jpk,jpts) :: zvts, zvab! vertical profile of T & S , and alpha & betaat 1 given point74 REAL(wp), DIMENSION( jpi,jpj,jpk ) :: zn2 ! N^275 REAL(wp), DIMENSION( jpi,jpj,jpk,jpts) :: zab! alpha and beta75 REAL(wp), DIMENSION( jpk ) :: zvn2 ! vertical profile of N2 at 1 given point... 76 REAL(wp), DIMENSION( jpk,jpts) :: zvts, zvab ! vertical profile of T & S , and alpha & betaat 1 given point 77 REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zn2 ! N^2 78 REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zab ! alpha and beta 76 79 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds ! 3D workspace 77 80 ! 78 81 LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is 79 82 INTEGER :: ilc1, jlc1, klc1, nncpu ! actually happening in a water column at point "ilc1, jlc1" 83 INTEGER :: isi, isj, iei, iej 80 84 LOGICAL :: lp_monitor_point = .FALSE. ! in CPU domain "nncpu" 81 85 !!---------------------------------------------------------------------- … … 87 91 IF( l_trdtra ) THEN !* Save initial after fields 88 92 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 89 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 93 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 90 94 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 91 95 ENDIF … … 101 105 CALL bn2 ( pts(:,:,:,:,Kaa), zab, zn2, Kmm ) ! after Brunt-Vaisala (given on W-points) 102 106 ! 103 inpcc = 0 104 ! 105 DO_2D( 0, 0, 0, 0 ) ! interior column only 107 IF( ntile == 0 .OR. ntile == 1 ) nnpcc = 0 ! Do only on the first tile 108 ! 109 IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 110 IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 111 IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 112 IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 113 ! 114 ! NOTE: [tiling-comms-merge] Bounds changed to avoid repeating this calculation for overlapping rows when using tiling 115 DO_2D( isj, iej, isi, iei ) ! interior column only 106 116 ! 107 117 IF( tmask(ji,jj,2) == 1 ) THEN ! At least 2 ocean points … … 160 170 ENDIF 161 171 ! 162 IF( jiter == 1 ) inpcc = inpcc + 1172 IF( jiter == 1 ) nnpcc = nnpcc + 1 163 173 ! 164 174 IF( lp_monitor_point ) WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer … … 310 320 ENDIF 311 321 ! 312 CALL lbc_lnk_multi( 'tranpc', pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp )313 !314 IF( lwp .AND. l_LB_debug ) THEN315 WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc316 WRITE(numout,*)322 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 323 IF( lwp .AND. l_LB_debug ) THEN 324 WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', nnpcc 325 WRITE(numout,*) 326 ENDIF 317 327 ENDIF 318 328 ! -
NEMO/trunk/src/OCE/TRA/traqsr.F90
r13970 r13982 22 22 USE phycst ! physical constants 23 23 USE dom_oce ! ocean space and time domain 24 USE domain, ONLY : dom_tile 24 25 USE sbc_oce ! surface boundary condition: ocean 25 26 USE trc_oce ! share SMS/Ocean variables … … 107 108 ! 108 109 INTEGER :: ji, jj, jk ! dummy loop indices 109 INTEGER :: irgb 110 INTEGER :: irgb, isi, iei, isj, iej ! local integers 110 111 REAL(wp) :: zchl, zcoef, z1_2 ! local scalars 111 112 REAL(wp) :: zc0 , zc1 , zc2 , zc3 ! - - … … 120 121 IF( ln_timing ) CALL timing_start('tra_qsr') 121 122 ! 122 IF( kt == nit000 ) THEN 123 IF(lwp) WRITE(numout,*) 124 IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation' 125 IF(lwp) WRITE(numout,*) '~~~~~~~' 123 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 124 IF( kt == nit000 ) THEN 125 IF(lwp) WRITE(numout,*) 126 IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation' 127 IF(lwp) WRITE(numout,*) '~~~~~~~' 128 ENDIF 126 129 ENDIF 127 130 ! 128 131 IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend 129 ALLOCATE( ztrdt(jpi,jpj,jpk) ) 132 ALLOCATE( ztrdt(jpi,jpj,jpk) ) 130 133 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 131 134 ENDIF … … 134 137 ! ! before qsr induced heat content ! 135 138 ! !-----------------------------------! 139 ! NOTE: [tiling-comms-merge] Many DO loop bounds changed (probably more than necessary) to avoid changing results when using tiling. Some bounds were also adjusted to account for those changed in tra_atf 140 IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 141 IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 142 IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 143 IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 144 136 145 IF( kt == nit000 ) THEN !== 1st time step ==! 137 146 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 .AND. .NOT.l_1st_euler ) THEN ! read in restart 138 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file'139 147 z1_2 = 0.5_wp 140 CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux 148 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 149 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file' 150 CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux 151 ENDIF 141 152 ELSE ! No restart or restart not found: Euler forward time stepping 142 153 z1_2 = 1._wp 143 qsr_hc_b(:,:,:) = 0._wp 154 DO_3D( isj, iej, isi, iei, 1, jpk ) 155 qsr_hc_b(ji,jj,jk) = 0._wp 156 END_3D 144 157 ENDIF 145 158 ELSE !== Swap of qsr heat content ==! 146 159 z1_2 = 0.5_wp 147 qsr_hc_b(:,:,:) = qsr_hc(:,:,:) 160 DO_3D( isj, iej, isi, iei, 1, jpk ) 161 qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk) 162 END_3D 148 163 ENDIF 149 164 ! … … 154 169 CASE( np_BIO ) !== bio-model fluxes ==! 155 170 ! 156 DO jk = 1, nksr157 qsr_hc( :,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )158 END DO171 DO_3D( isj, iej, isi, iei, 1, nksr ) 172 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) 173 END_3D 159 174 ! 160 175 CASE( np_RGB , np_RGBc ) !== R-G-B fluxes ==! 161 176 ! 162 ALLOCATE( ze0 ( jpi,jpj) , ze1 (jpi,jpj) , &163 & ze2 ( jpi,jpj) , ze3 (jpi,jpj) , &164 & ztmp3d( jpi,jpj,nksr + 1) )177 ALLOCATE( ze0 (A2D(nn_hls)) , ze1 (A2D(nn_hls)) , & 178 & ze2 (A2D(nn_hls)) , ze3 (A2D(nn_hls)) , & 179 & ztmp3d(A2D(nn_hls),nksr + 1) ) 165 180 ! 166 181 IF( nqsr == np_RGBc ) THEN !* Variable Chlorophyll 167 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 182 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only for the full domain 183 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Use full domain 184 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 185 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) ! Revert to tile domain 186 ENDIF 168 187 ! 169 188 ! Separation in R-G-B depending on the surface Chl … … 172 191 ! most expensive calculations) 173 192 ! 174 DO_2D( 0, 0, 0, 0)193 DO_2D( isj, iej, isi, iei ) 175 194 ! zlogc = log(zchl) 176 195 zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) … … 191 210 192 211 ! 193 DO_3D( 0, 0, 0, 0, 1, nksr + 1 )212 DO_3D( isj, iej, isi, iei, 1, nksr + 1 ) 194 213 ! zchl = ALOG( ze0(ji,jj) ) 195 214 zlogc = ze0(ji,jj) … … 216 235 zlui = 41 + 20.*LOG10(zchl) + 1.e-15 217 236 DO jk = 1, nksr + 1 218 ztmp3d(:,:,jk) = zlui 237 ztmp3d(:,:,jk) = zlui 219 238 END DO 220 239 ENDIF 221 240 ! 222 241 zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B 223 DO_2D( 0, 0, 0, 0)242 DO_2D( isj, iej, isi, iei ) 224 243 ze0(ji,jj) = rn_abs * qsr(ji,jj) 225 244 ze1(ji,jj) = zcoef * qsr(ji,jj) … … 232 251 ! 233 252 ! !* interior equi-partition in R-G-B depending on vertical profile of Chl 234 DO_3D( 0, 0, 0, 0, 2, nksr + 1 )253 DO_3D( isj, iej, isi, iei, 2, nksr + 1 ) 235 254 ze3t = e3t(ji,jj,jk-1,Kmm) 236 255 irgb = NINT( ztmp3d(ji,jj,jk) ) … … 246 265 END_3D 247 266 ! 248 DO_3D( 0, 0, 0, 0, 1, nksr ) !* now qsr induced heat content267 DO_3D( isj, iej, isi, iei, 1, nksr ) !* now qsr induced heat content 249 268 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 250 269 END_3D … … 256 275 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 257 276 zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 258 DO_3D( 0, 0, 0, 0, 1, nksr ) ! solar heat absorbed at T-point in the top 400m277 DO_3D( isj, iej, isi, iei, 1, nksr ) !* now qsr induced heat content 259 278 zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r ) 260 279 zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) … … 274 293 ! 275 294 ! sea-ice: store the 1st ocean level attenuation coefficient 276 DO_2D( 0, 0, 0, 0)295 DO_2D( isj, iej, isi, iei ) 277 296 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 278 297 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 279 298 ENDIF 280 299 END_2D 281 CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp ) 282 ! 283 IF( iom_use('qsr3d') ) THEN ! output the shortwave Radiation distribution 284 ALLOCATE( zetot(jpi,jpj,jpk) ) 285 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 286 DO jk = nksr, 1, -1 287 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 288 END DO 289 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation 290 DEALLOCATE( zetot ) 291 ENDIF 292 ! 293 IF( lrst_oce ) THEN ! write in the ocean restart file 294 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 295 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) 300 ! 301 ! TEMP: [tiling] This change not necessary and working array can use A2D(nn_hls) if using XIOS (subdomain support) 302 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 303 IF( iom_use('qsr3d') ) THEN ! output the shortwave Radiation distribution 304 ALLOCATE( zetot(jpi,jpj,jpk) ) 305 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 306 DO jk = nksr, 1, -1 307 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 308 END DO 309 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation 310 DEALLOCATE( zetot ) 311 ENDIF 312 ENDIF 313 ! 314 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 315 IF( lrst_oce ) THEN ! write in the ocean restart file 316 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 317 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) 318 ENDIF 296 319 ENDIF 297 320 ! … … 299 322 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 300 323 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 301 DEALLOCATE( ztrdt ) 324 DEALLOCATE( ztrdt ) 302 325 ENDIF 303 326 ! ! print mean trends (used for debugging) -
NEMO/trunk/src/OCE/TRA/trasbc.F90
r13970 r13982 76 76 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 77 77 ! 78 INTEGER :: ji, jj, jk, jn ! dummy loop indices79 INTEGER :: ikt, ikb 80 REAL(wp) :: zfact, z1_e3t, zdep, ztim ! local scalar78 INTEGER :: ji, jj, jk, jn ! dummy loop indices 79 INTEGER :: ikt, ikb, isi, iei, isj, iej ! local integers 80 REAL(wp) :: zfact, z1_e3t, zdep, ztim ! local scalar 81 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds 82 82 !!---------------------------------------------------------------------- … … 84 84 IF( ln_timing ) CALL timing_start('tra_sbc') 85 85 ! 86 IF( kt == nit000 ) THEN 87 IF(lwp) WRITE(numout,*) 88 IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition' 89 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 86 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 87 IF( kt == nit000 ) THEN 88 IF(lwp) WRITE(numout,*) 89 IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition' 90 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 91 ENDIF 90 92 ENDIF 91 93 ! 92 94 IF( l_trdtra ) THEN !* Save ta and sa trends 93 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )95 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 94 96 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 95 97 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 96 98 ENDIF 97 99 ! 100 ! NOTE: [tiling-comms-merge] Many DO loop bounds changed to avoid changing results when using tiling. Some bounds were also adjusted to account for those changed in tra_atf 101 IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 102 IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 103 IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 104 IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 105 98 106 !!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist) 99 107 IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration 100 qns(:,:) = qns(:,:) + qsr(:,:) ! total heat flux in qns 101 qsr(:,:) = 0._wp ! qsr set to zero 108 DO_2D( isj, iej, isi, iei ) 109 qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns 110 qsr(ji,jj) = 0._wp ! qsr set to zero 111 END_2D 102 112 ENDIF 103 113 … … 109 119 IF( ln_rstart .AND. & ! Restart: read in restart file 110 120 & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 111 IF(lwp) WRITE(numout,*) ' nit000-1 sbc tracer content field read in the restart file'112 121 zfact = 0.5_wp 113 sbc_tsc(:,:,:) = 0._wp 114 CALL iom_get( numror, jpdom_auto, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) ) ! before heat content sbc trend 115 CALL iom_get( numror, jpdom_auto, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) ) ! before salt content sbc trend 122 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 123 IF(lwp) WRITE(numout,*) ' nit000-1 sbc tracer content field read in the restart file' 124 sbc_tsc(:,:,:) = 0._wp 125 CALL iom_get( numror, jpdom_auto, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) ) ! before heat content sbc trend 126 CALL iom_get( numror, jpdom_auto, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) ) ! before salt content sbc trend 127 ENDIF 116 128 ELSE ! No restart or restart not found: Euler forward time stepping 117 129 zfact = 1._wp 118 sbc_tsc(:,:,:) = 0._wp 119 sbc_tsc_b(:,:,:) = 0._wp 130 DO_2D( isj, iej, isi, iei ) 131 sbc_tsc(ji,jj,:) = 0._wp 132 sbc_tsc_b(ji,jj,:) = 0._wp 133 END_2D 120 134 ENDIF 121 135 ELSE !* other time-steps: swap of forcing fields 122 136 zfact = 0.5_wp 123 sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:) 137 DO_2D( isj, iej, isi, iei ) 138 sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:) 139 END_2D 124 140 ENDIF 125 141 ! !== Now sbc tracer content fields ==! 126 DO_2D( 0, 1, 0, 0)142 DO_2D( isj, iej, isi, iei ) 127 143 sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux 128 144 sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting 129 145 END_2D 130 146 IF( ln_linssh ) THEN !* linear free surface 131 DO_2D( 0, 1, 0, 0) !==>> add concentration/dilution effect due to constant volume cell147 DO_2D( isj, iej, isi, iei ) !==>> add concentration/dilution effect due to constant volume cell 132 148 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 133 149 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm) 134 150 END_2D !==>> output c./d. term 135 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 136 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) ) 151 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 152 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 153 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) ) 154 ENDIF 137 155 ENDIF 138 156 ! 139 157 DO jn = 1, jpts !== update tracer trend ==! 140 DO_2D( 0, 1, 0, 0 ) 158 ! NOTE: [tiling-comms-merge] This looped over nn_hls, which changes the results when using tiling 159 DO_2D( 0, 0, 0, 0 ) 141 160 pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) & 142 161 & / e3t(ji,jj,1,Kmm) … … 144 163 END DO 145 164 ! 146 IF( lrst_oce ) THEN !== write sbc_tsc in the ocean restart file ==! 147 CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) ) 148 CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) ) 165 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 166 IF( lrst_oce ) THEN !== write sbc_tsc in the ocean restart file ==! 167 CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) ) 168 CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) ) 169 ENDIF 149 170 ENDIF 150 171 ! … … 155 176 IF( ln_rnf ) THEN ! input of heat and salt due to river runoff 156 177 zfact = 0.5_wp 157 DO_2D( 0, 1, 0, 0 )178 DO_2D( 0, 0, 0, 0 ) 158 179 IF( rnf(ji,jj) /= 0._wp ) THEN 159 180 zdep = zfact / h_rnf(ji,jj) … … 168 189 ENDIF 169 190 170 IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) ) ! runoff term on sst 171 IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) ) ! runoff term on sss 191 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 192 IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) ) ! runoff term on sst 193 IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) ) ! runoff term on sss 194 ENDIF 172 195 173 196 #if defined key_asminc … … 180 203 ! 181 204 IF( ln_linssh ) THEN 182 DO_2D( 0, 1, 0, 0 )205 DO_2D( 0, 0, 0, 0 ) 183 206 ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm) 184 207 pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim … … 186 209 END_2D 187 210 ELSE 188 DO_2D( 0, 1, 0, 0 )211 DO_2D( 0, 0, 0, 0 ) 189 212 ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) ) 190 213 pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim … … 202 225 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt ) 203 226 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds ) 204 DEALLOCATE( ztrdt , ztrds ) 227 DEALLOCATE( ztrdt , ztrds ) 205 228 ENDIF 206 229 ! -
NEMO/trunk/src/OCE/TRA/trazdf.F90
r13497 r13982 13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 15 USE dom_oce ! ocean space and time domain variables 16 16 USE domvvl ! variable volume 17 17 USE phycst ! physical constant … … 55 55 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 56 56 ! 57 INTEGER :: j k ! Dummy loop indices57 INTEGER :: ji, jj, jk ! Dummy loop indices 58 58 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace 59 59 !!--------------------------------------------------------------------- … … 62 62 ! 63 63 IF( kt == nit000 ) THEN 64 IF(lwp)WRITE(numout,*) 65 IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S' 66 IF(lwp)WRITE(numout,*) '~~~~~~~ ' 64 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 65 IF(lwp)WRITE(numout,*) 66 IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S' 67 IF(lwp)WRITE(numout,*) '~~~~~~~ ' 68 ENDIF 67 69 ENDIF 68 70 ! 69 71 IF( l_trdtra ) THEN !* Save ta and sa trends 70 ALLOCATE( ztrdt(jpi,jpj,jpk) 72 ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) ) 71 73 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 72 74 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) … … 80 82 ! JMM avoid negative salinities near river outlet ! Ugly fix 81 83 ! JMM : restore negative salinities to small salinities: 82 WHERE( pts( :,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp84 WHERE( pts(A2D(0),:,jp_sal,Kaa) < 0._wp ) pts(A2D(0),:,jp_sal,Kaa) = 0.1_wp 83 85 !!gm 84 86 85 87 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 86 DO jk = 1, jpk m188 DO jk = 1, jpk 87 89 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) & 88 90 & - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & … … 94 96 & - ztrds(:,:,jk) 95 97 END DO 98 ! NOTE: [tiling-comms-merge] The diagnostic results change along the north fold if this is removed 96 99 !!gm this should be moved in trdtra.F90 and done on all trends 97 100 CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1.0_wp , ztrds, 'T', 1.0_wp ) … … 140 143 INTEGER :: ji, jj, jk, jn ! dummy loop indices 141 144 REAL(wp) :: zrhs, zzwi, zzws ! local scalars 142 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwi, zwt, zwd, zws145 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwt, zwd, zws 143 146 !!--------------------------------------------------------------------- 144 147 ! … … 154 157 ! 155 158 ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers 156 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt(:,:,2:jpk) 157 ELSE ; zwt(:,:,2:jpk) = avs(:,:,2:jpk) 159 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 160 DO_3D( 1, 1, 1, 1, 2, jpk ) 161 zwt(ji,jj,jk) = avt(ji,jj,jk) 162 END_3D 163 ELSE 164 DO_3D( 1, 1, 1, 1, 2, jpk ) 165 zwt(ji,jj,jk) = avs(ji,jj,jk) 166 END_3D 158 167 ENDIF 159 168 zwt(:,:,1) = 0._wp … … 222 231 END_2D 223 232 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 224 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) & 233 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) & 225 234 & + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side 226 235 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) -
NEMO/trunk/src/OCE/TRA/zpshde.F90
r13497 r13982 17 17 USE oce ! ocean: dynamics and tracers variables 18 18 USE dom_oce ! domain: ocean variables 19 USE domutl, ONLY : is_tile 19 20 USE phycst ! physical constants 20 21 USE eosbn2 ! ocean equation of state … … 40 41 CONTAINS 41 42 42 SUBROUTINE zps_hde( kt, Kmm, kjpt, pta, pgtu, pgtv, & 43 & prd, pgru, pgrv ) 43 SUBROUTINE zps_hde( kt, Kmm, kjpt, pta, pgtu, pgtv, & 44 & prd, pgru, pgrv ) 45 !! 46 INTEGER , INTENT(in ) :: kt ! ocean time-step index 47 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 48 INTEGER , INTENT(in ) :: kjpt ! number of tracers 49 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pta ! 4D tracers fields 50 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 51 REAL(wp), DIMENSION(:,:,:) , INTENT(inout), OPTIONAL :: prd ! 3D density anomaly fields 52 REAL(wp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 53 ! 54 INTEGER :: itrd, itgr 55 !! 56 IF( PRESENT(prd) ) THEN ; itrd = is_tile(prd) ; ELSE ; itrd = 0 ; ENDIF 57 IF( PRESENT(pgru) ) THEN ; itgr = is_tile(pgru) ; ELSE ; itgr = 0 ; ENDIF 58 59 CALL zps_hde_t( kt, Kmm, kjpt, pta, is_tile(pta), pgtu, pgtv, is_tile(pgtu), & 60 & prd, itrd, pgru, pgrv, itgr ) 61 END SUBROUTINE zps_hde 62 63 64 SUBROUTINE zps_hde_t( kt, Kmm, kjpt, pta, ktta, pgtu, pgtv, ktgt, & 65 & prd, ktrd, pgru, pgrv, ktgr ) 44 66 !!---------------------------------------------------------------------- 45 67 !! *** ROUTINE zps_hde *** … … 85 107 !! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 86 108 !!---------------------------------------------------------------------- 87 INTEGER , INTENT(in ) :: kt ! ocean time-step index 88 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 89 INTEGER , INTENT(in ) :: kjpt ! number of tracers 90 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 91 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 92 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 93 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 109 INTEGER , INTENT(in ) :: kt ! ocean time-step index 110 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 111 INTEGER , INTENT(in ) :: kjpt ! number of tracers 112 INTEGER , INTENT(in ) :: ktta, ktgt, ktrd, ktgr 113 REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(inout) :: pta ! 4D tracers fields 114 REAL(wp), DIMENSION(A2D_T(ktgt) ,KJPT), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 115 REAL(wp), DIMENSION(A2D_T(ktrd),JPK ), INTENT(inout), OPTIONAL :: prd ! 3D density anomaly fields 116 REAL(wp), DIMENSION(A2D_T(ktgr) ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 94 117 ! 95 118 INTEGER :: ji, jj, jn ! Dummy loop indices 96 119 INTEGER :: iku, ikv, ikum1, ikvm1 ! partial step level (ocean bottom level) at u- and v-points 97 120 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! local scalars 98 REAL(wp), DIMENSION( jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos99 REAL(wp), DIMENSION( jpi,jpj,kjpt) :: zti, ztj !121 REAL(wp), DIMENSION(A2D(nn_hls)) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 122 REAL(wp), DIMENSION(A2D(nn_hls),kjpt) :: zti, ztj ! 100 123 !!---------------------------------------------------------------------- 101 124 ! 102 125 IF( ln_timing ) CALL timing_start( 'zps_hde') 126 ! NOTE: [tiling-comms-merge] Some lbc_lnks in tra_adv and tra_ldf can be taken out in the zps case, because this lbc_lnk is called when zps_hde is called in the stp routine. In the zco case they are still needed. 127 IF (nn_hls.EQ.2) THEN 128 CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp) 129 IF(PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp) 130 END IF 103 131 ! 104 132 pgtu(:,:,:) = 0._wp ; zti (:,:,:) = 0._wp ; zhi (:,:) = 0._wp … … 107 135 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 108 136 ! 109 DO_2D( 1, 0, 1, 0 )137 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! Gradient of density at the last level 110 138 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 111 139 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 … … 146 174 END DO 147 175 ! 148 CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond.176 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond. 149 177 ! 150 178 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) 151 179 pgru(:,:) = 0._wp 152 180 pgrv(:,:) = 0._wp ! depth of the partial step level 153 DO_2D( 1, 0, 1, 0)181 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 154 182 iku = mbku(ji,jj) 155 183 ikv = mbkv(ji,jj) … … 167 195 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 168 196 ! 169 DO_2D( 1, 0, 1, 0) ! Gradient of density at the last level197 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! Gradient of density at the last level 170 198 iku = mbku(ji,jj) 171 199 ikv = mbkv(ji,jj) … … 179 207 ENDIF 180 208 END_2D 181 CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions209 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions 182 210 ! 183 211 END IF … … 185 213 IF( ln_timing ) CALL timing_stop( 'zps_hde') 186 214 ! 187 END SUBROUTINE zps_hde 215 END SUBROUTINE zps_hde_t 188 216 189 217 190 218 SUBROUTINE zps_hde_isf( kt, Kmm, kjpt, pta, pgtu, pgtv, pgtui, pgtvi, & 191 & prd, pgru, pgrv, pgrui, pgrvi ) 219 & prd, pgru, pgrv, pgrui, pgrvi ) 220 !! 221 INTEGER , INTENT(in ) :: kt ! ocean time-step index 222 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 223 INTEGER , INTENT(in ) :: kjpt ! number of tracers 224 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pta ! 4D tracers fields 225 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 226 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 227 REAL(wp), DIMENSION(:,:,:) , INTENT(inout), OPTIONAL :: prd ! 3D density anomaly fields 228 REAL(wp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 229 REAL(wp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 230 ! 231 INTEGER :: itrd, itgr, itgri 232 !! 233 IF( PRESENT(prd) ) THEN ; itrd = is_tile(prd) ; ELSE ; itrd = 0 ; ENDIF 234 IF( PRESENT(pgru) ) THEN ; itgr = is_tile(pgru) ; ELSE ; itgr = 0 ; ENDIF 235 IF( PRESENT(pgrui) ) THEN ; itgri = is_tile(pgrui) ; ELSE ; itgri = 0 ; ENDIF 236 237 CALL zps_hde_isf_t( kt, Kmm, kjpt, pta, is_tile(pta), pgtu, pgtv, is_tile(pgtu), pgtui, pgtvi, is_tile(pgtui), & 238 & prd, itrd, pgru, pgrv, itgr, pgrui, pgrvi, itgri ) 239 END SUBROUTINE zps_hde_isf 240 241 242 SUBROUTINE zps_hde_isf_t( kt, Kmm, kjpt, pta, ktta, pgtu, pgtv, ktgt, pgtui, pgtvi, ktgti, & 243 & prd, ktrd, pgru, pgrv, ktgr, pgrui, pgrvi, ktgri ) 192 244 !!---------------------------------------------------------------------- 193 245 !! *** ROUTINE zps_hde_isf *** … … 236 288 !! - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 237 289 !!---------------------------------------------------------------------- 238 INTEGER , INTENT(in ) :: kt ! ocean time-step index 239 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 240 INTEGER , INTENT(in ) :: kjpt ! number of tracers 241 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 242 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 243 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 244 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 245 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 246 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 290 INTEGER , INTENT(in ) :: kt ! ocean time-step index 291 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 292 INTEGER , INTENT(in ) :: kjpt ! number of tracers 293 INTEGER , INTENT(in ) :: ktta, ktgt, ktgti, ktrd, ktgr, ktgri 294 REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(inout) :: pta ! 4D tracers fields 295 REAL(wp), DIMENSION(A2D_T(ktgt) ,KJPT), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 296 REAL(wp), DIMENSION(A2D_T(ktgti) ,KJPT), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 297 REAL(wp), DIMENSION(A2D_T(ktrd),JPK ), INTENT(inout), OPTIONAL :: prd ! 3D density anomaly fields 298 REAL(wp), DIMENSION(A2D_T(ktgr) ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 299 REAL(wp), DIMENSION(A2D_T(ktgri) ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 247 300 ! 248 301 INTEGER :: ji, jj, jn ! Dummy loop indices 249 302 INTEGER :: iku, ikv, ikum1, ikvm1,ikup1, ikvp1 ! partial step level (ocean bottom level) at u- and v-points 250 303 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars 251 REAL(wp), DIMENSION( jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos252 REAL(wp), DIMENSION( jpi,jpj,kjpt) :: zti, ztj !304 REAL(wp), DIMENSION(A2D(nn_hls)) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 305 REAL(wp), DIMENSION(A2D(nn_hls),kjpt) :: zti, ztj ! 253 306 !!---------------------------------------------------------------------- 254 307 ! 255 308 IF( ln_timing ) CALL timing_start( 'zps_hde_isf') 256 309 ! 310 IF (nn_hls.EQ.2) THEN 311 CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp) 312 IF (PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp) 313 END IF 314 257 315 pgtu (:,:,:) = 0._wp ; pgtv (:,:,:) =0._wp 258 316 pgtui(:,:,:) = 0._wp ; pgtvi(:,:,:) =0._wp … … 262 320 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 263 321 ! 264 DO_2D( 1, 0, 1, 0)322 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 265 323 266 324 iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points … … 302 360 END DO 303 361 ! 304 CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond.362 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond. 305 363 306 364 ! horizontal derivative of density anomalies (rd) … … 308 366 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 309 367 ! 310 DO_2D( 1, 0, 1, 0)368 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 311 369 312 370 iku = mbku(ji,jj) … … 329 387 CALL eos( ztj, zhj, zrj ) 330 388 331 DO_2D( 1, 0, 1, 0 ) ! Gradient of density at the last level389 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 332 390 iku = mbku(ji,jj) 333 391 ikv = mbkv(ji,jj) … … 344 402 END_2D 345 403 346 CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions404 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions 347 405 ! 348 406 END IF … … 351 409 ! 352 410 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! ! 353 DO_2D( 1, 0, 1, 0)411 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 354 412 iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 355 413 ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 … … 395 453 ! 396 454 END DO 397 CALL lbc_lnk_multi( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond.455 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond. 398 456 399 457 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) 400 458 ! 401 459 pgrui(:,:) =0.0_wp; pgrvi(:,:) =0.0_wp; 402 DO_2D( 1, 0, 1, 0)460 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 403 461 404 462 iku = miku(ji,jj) … … 420 478 CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj 421 479 ! 422 DO_2D( 1, 0, 1, 0 ) ! Gradient of density at the last level480 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 423 481 iku = miku(ji,jj) 424 482 ikv = mikv(ji,jj) … … 434 492 435 493 END_2D 436 CALL lbc_lnk_multi( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp ) ! Lateral boundary conditions494 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp ) ! Lateral boundary conditions 437 495 ! 438 496 END IF … … 440 498 IF( ln_timing ) CALL timing_stop( 'zps_hde_isf') 441 499 ! 442 END SUBROUTINE zps_hde_isf 500 END SUBROUTINE zps_hde_isf_t 443 501 444 502 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.