Changeset 7698 for trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2017-02-18T10:02:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r6497 r7698 171 171 !!---------------------------------------------------------------------- 172 172 INTEGER, INTENT(in) :: kt ! ocean time step 173 INTEGER :: jk, jj, ji 173 174 !!---------------------------------------------------------------------- 174 175 ! … … 179 180 ! 180 181 IF( kt /= nit000 ) THEN ! restore before value to compute tke 181 avt (:,:,:) = avt_k (:,:,:) 182 avm (:,:,:) = avm_k (:,:,:) 183 avmu(:,:,:) = avmu_k(:,:,:) 184 avmv(:,:,:) = avmv_k(:,:,:) 182 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 183 DO jk = 1, jpk 184 DO jj = 1, jpj 185 DO ji = 1, jpi 186 avt (ji,jj,jk) = avt_k (ji,jj,jk) 187 avm (ji,jj,jk) = avm_k (ji,jj,jk) 188 avmu(ji,jj,jk) = avmu_k(ji,jj,jk) 189 avmv(ji,jj,jk) = avmv_k(ji,jj,jk) 190 END DO 191 END DO 192 END DO 185 193 ENDIF 186 194 ! … … 189 197 CALL tke_avn ! now avt, avm, avmu, avmv 190 198 ! 191 avt_k (:,:,:) = avt (:,:,:) 192 avm_k (:,:,:) = avm (:,:,:) 193 avmu_k(:,:,:) = avmu(:,:,:) 194 avmv_k(:,:,:) = avmv(:,:,:) 199 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 200 DO jk = 1, jpk 201 DO jj = 1, jpj 202 DO ji = 1, jpi 203 avt_k (ji,jj,jk) = avt (ji,jj,jk) 204 avm_k (ji,jj,jk) = avm (ji,jj,jk) 205 avmu_k(ji,jj,jk) = avmu(ji,jj,jk) 206 avmv_k(ji,jj,jk) = avmv(ji,jj,jk) 207 END DO 208 END DO 209 END DO 195 210 ! 196 211 #if defined key_agrif … … 253 268 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 254 269 IF ( ln_isfcav ) THEN 270 !$OMP PARALLEL DO schedule(static) private(jj, ji) 255 271 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 256 272 DO ji = fs_2, fs_jpim1 ! vector opt. … … 259 275 END DO 260 276 END IF 277 !$OMP PARALLEL DO schedule(static) private(jj, ji) 261 278 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 262 279 DO ji = fs_2, fs_jpim1 ! vector opt. … … 293 310 ! 294 311 ! !* total energy produce by LC : cumulative sum over jk 295 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1) 312 !$OMP PARALLEL 313 !$OMP DO schedule(static) private(jj, ji) 314 DO jj =1, jpj 315 DO ji=1, jpi 316 zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * gdepw_n(ji,jj,1) * e3w_n(ji,jj,1) 317 END DO 318 END DO 296 319 DO jk = 2, jpk 297 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk) 320 !$OMP DO schedule(static) private(jj, ji) 321 DO jj =1, jpj 322 DO ji=1, jpi 323 zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk) 324 END DO 325 END DO 298 326 END DO 299 327 ! !* finite Langmuir Circulation depth 300 328 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 301 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 329 !$OMP DO schedule(static) private(jj,ji) 330 DO jj = 1, jpj 331 DO ji = 1, jpi 332 imlc(ji,jj) = mbkt(ji,jj) + 1 ! Initialization to the number of w ocean point (=2 over land) 333 END DO 334 END DO 302 335 DO jk = jpkm1, 2, -1 336 !$OMP DO schedule(static) private(jj, ji, zus) 303 337 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 304 338 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) … … 309 343 END DO 310 344 ! ! finite LC depth 345 !$OMP DO schedule(static) private(jj, ji) 311 346 DO jj = 1, jpj 312 347 DO ji = 1, jpi … … 315 350 END DO 316 351 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 352 !$OMP DO schedule(static) private(jk, jj, ji, zus, zind, zwlc) 317 353 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 318 354 DO jj = 2, jpjm1 … … 328 364 END DO 329 365 END DO 366 !$OMP END PARALLEL 330 367 ! 331 368 ENDIF … … 338 375 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 339 376 ! 377 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji) 340 378 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 341 379 DO jj = 1, jpjm1 … … 356 394 ! Note that zesh2 is also computed in the next loop. 357 395 ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 396 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zesh2, zri) 358 397 DO jk = 2, jpkm1 359 398 DO jj = 2, jpjm1 … … 372 411 ENDIF 373 412 ! 413 !$OMP PARALLEL 414 !$OMP DO schedule(static) private(jk, jj, ji, zcof, zzd_up, zzd_lw, zesh2) 374 415 DO jk = 2, jpkm1 !* Matrix and right hand side in en 375 416 DO jj = 2, jpjm1 … … 405 446 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 406 447 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 448 !$OMP DO schedule(static) private(jj, ji) 407 449 DO jj = 2, jpjm1 408 450 DO ji = fs_2, fs_jpim1 ! vector opt. … … 411 453 END DO 412 454 END DO 455 !$OMP DO schedule(static) private(jj, ji) 413 456 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 414 457 DO ji = fs_2, fs_jpim1 ! vector opt. … … 417 460 END DO 418 461 DO jk = 3, jpkm1 462 !$OMP DO schedule(static) private(jj, ji) 419 463 DO jj = 2, jpjm1 420 464 DO ji = fs_2, fs_jpim1 ! vector opt. … … 423 467 END DO 424 468 END DO 469 !$OMP DO schedule(static) private(jj, ji) 425 470 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 426 471 DO ji = fs_2, fs_jpim1 ! vector opt. … … 429 474 END DO 430 475 DO jk = jpk-2, 2, -1 476 !$OMP DO schedule(static) private(jj, ji) 431 477 DO jj = 2, jpjm1 432 478 DO ji = fs_2, fs_jpim1 ! vector opt. … … 435 481 END DO 436 482 END DO 483 !$OMP DO schedule(static) private(jk,jj, ji) 437 484 DO jk = 2, jpkm1 ! set the minimum value of tke 438 485 DO jj = 2, jpjm1 … … 442 489 END DO 443 490 END DO 491 !$OMP END PARALLEL 444 492 445 493 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 450 498 451 499 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 500 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji) 452 501 DO jk = 2, jpkm1 453 502 DO jj = 2, jpjm1 … … 459 508 END DO 460 509 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 510 !$OMP PARALLEL DO schedule(static) private(jj, ji, jk) 461 511 DO jj = 2, jpjm1 462 512 DO ji = fs_2, fs_jpim1 ! vector opt. … … 467 517 END DO 468 518 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 519 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ztx2, zty2, ztau, zdif) 469 520 DO jk = 2, jpkm1 470 521 DO jj = 2, jpjm1 … … 545 596 ! 546 597 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 547 zmxlm(:,:,:) = rmxl_min 548 zmxld(:,:,:) = rmxl_min 598 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 599 DO jk = 1, jpk 600 DO jj = 1, jpj 601 DO ji = 1, jpi 602 zmxlm(ji,jj,jk) = rmxl_min 603 zmxld(ji,jj,jk) = rmxl_min 604 END DO 605 END DO 606 END DO 549 607 ! 550 608 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 609 !$OMP PARALLEL DO schedule(static) private(jj, ji, zraug) 551 610 DO jj = 2, jpjm1 552 611 DO ji = fs_2, fs_jpim1 … … 556 615 END DO 557 616 ELSE 558 zmxlm(:,:,1) = rn_mxl0 617 !$OMP PARALLEL DO schedule(static) private(jj,ji) 618 DO jj = 1, jpj 619 DO ji = 1, jpi 620 zmxlm(ji,jj,1) = rn_mxl0 621 END DO 622 END DO 559 623 ENDIF 560 624 ! 625 !$OMP PARALLEL 626 !$OMP DO schedule(static) private(jk, jj, ji, zrn2) 561 627 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 562 628 DO jj = 2, jpjm1 … … 570 636 ! !* Physical limits for the mixing length 571 637 ! 572 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 573 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 638 !$OMP DO schedule(static) private(jj,ji) 639 DO jj = 1, jpj 640 DO ji = 1, jpi 641 zmxld(ji,jj, 1 ) = zmxlm(ji,jj,1) ! surface set to the minimum value 642 zmxld(ji,jj,jpk) = rmxl_min ! last level set to the minimum value 643 END DO 644 END DO 645 !$OMP END PARALLEL 574 646 ! 575 647 SELECT CASE ( nn_mxl ) … … 578 650 ! where wmask = 0 set zmxlm == e3w_n 579 651 CASE ( 0 ) ! bounded by the distance to surface and bottom 652 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl) 580 653 DO jk = 2, jpkm1 581 654 DO jj = 2, jpjm1 … … 591 664 ! 592 665 CASE ( 1 ) ! bounded by the vertical scale factor 666 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl) 593 667 DO jk = 2, jpkm1 594 668 DO jj = 2, jpjm1 … … 602 676 ! 603 677 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 678 !$OMP PARALLEL 604 679 DO jk = 2, jpkm1 ! from the surface to the bottom : 680 !$OMP DO schedule(static) private(jj, ji) 605 681 DO jj = 2, jpjm1 606 682 DO ji = fs_2, fs_jpim1 ! vector opt. … … 610 686 END DO 611 687 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 688 !$OMP DO schedule(static) private(jj, ji, zemxl) 612 689 DO jj = 2, jpjm1 613 690 DO ji = fs_2, fs_jpim1 ! vector opt. … … 618 695 END DO 619 696 END DO 697 !$OMP END PARALLEL 620 698 ! 621 699 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 700 !$OMP PARALLEL 622 701 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 702 !$OMP DO schedule(static) private(jj, ji) 623 703 DO jj = 2, jpjm1 624 704 DO ji = fs_2, fs_jpim1 ! vector opt. … … 628 708 END DO 629 709 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 710 !$OMP DO schedule(static) private(jj, ji) 630 711 DO jj = 2, jpjm1 631 712 DO ji = fs_2, fs_jpim1 ! vector opt. … … 634 715 END DO 635 716 END DO 717 !$OMP DO schedule(static) private(jk, jj, ji, zemlm, zemlp) 636 718 DO jk = 2, jpkm1 637 719 DO jj = 2, jpjm1 … … 644 726 END DO 645 727 END DO 728 !$OMP END PARALLEL 646 729 ! 647 730 END SELECT 648 731 ! 649 732 # if defined key_c1d 650 e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales 651 e_mix(:,:,:) = zmxlm(:,:,:) 733 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 734 DO jk = 1, jpk 735 DO jj = 1, jpj 736 DO ji = 1, jpi 737 e_dis(ji,jj,jk) = zmxld(ji,jj,jk) ! c1d configuration : save mixing and dissipation turbulent length scales 738 e_mix(ji,jj,jk) = zmxlm(ji,jj,jk) 739 END DO 740 END DO 741 END DO 652 742 # endif 653 743 … … 655 745 ! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt) 656 746 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 747 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zsqen, zav) 657 748 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 658 749 DO jj = 2, jpjm1 … … 668 759 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 669 760 ! 761 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji) 670 762 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 671 763 DO jj = 2, jpjm1 … … 679 771 ! 680 772 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 773 !$OMP PARALLEL DO schedule(static) private(jk, jj, ji) 681 774 DO jk = 2, jpkm1 682 775 DO jj = 2, jpjm1 … … 798 891 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 799 892 CASE( 0 ) ! constant depth penetration (here 10 meters) 800 htau(:,:) = 10._wp 893 !$OMP PARALLEL DO schedule(static) private(jj,ji) 894 DO jj = 1, jpj 895 DO ji = 1, jpi 896 htau(ji,jj) = 10._wp 897 END DO 898 END DO 801 899 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 802 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 900 !$OMP PARALLEL DO schedule(static) private(jj,ji) 901 DO jj = 1, jpj 902 DO ji = 1, jpi 903 htau(ji,jj) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) ) 904 END DO 905 END DO 803 906 END SELECT 804 907 ENDIF 805 908 ! !* set vertical eddy coef. to the background value 909 !$OMP PARALLEL 910 !$OMP DO schedule(static) private(jk,jj,ji) 806 911 DO jk = 1, jpk 807 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 808 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 809 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 810 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 811 END DO 812 dissl(:,:,:) = 1.e-12_wp 912 DO jj = 1, jpj 913 DO ji = 1, jpi 914 avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk) 915 avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk) 916 avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk) 917 avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk) 918 END DO 919 END DO 920 END DO 921 !$OMP END DO NOWAIT 922 !$OMP DO schedule(static) private(jk,jj,ji) 923 DO jk = 1, jpk 924 DO jj = 1, jpj 925 DO ji = 1, jpi 926 dissl(ji,jj,jk) = 1.e-12_wp 927 END DO 928 END DO 929 END DO 930 !$OMP END PARALLEL 813 931 ! 814 932 CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files … … 830 948 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 831 949 ! 832 INTEGER :: jit, jk ! dummy loop indices950 INTEGER :: jit, jk, jj, ji ! dummy loop indices 833 951 INTEGER :: id1, id2, id3, id4, id5, id6 ! local integers 834 952 !!---------------------------------------------------------------------- … … 857 975 ELSE ! No TKE array found: initialisation 858 976 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 859 en (:,:,:) = rn_emin * tmask(:,:,:) 977 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 978 DO jk = 1, jpk 979 DO jj = 1, jpj 980 DO ji = 1, jpi 981 en (ji,jj,jk) = rn_emin * tmask(ji,jj,jk) 982 END DO 983 END DO 984 END DO 860 985 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 861 986 ! 862 avt_k (:,:,:) = avt (:,:,:) 863 avm_k (:,:,:) = avm (:,:,:) 864 avmu_k(:,:,:) = avmu(:,:,:) 865 avmv_k(:,:,:) = avmv(:,:,:) 987 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 988 DO jk = 1, jpk 989 DO jj = 1, jpj 990 DO ji = 1, jpi 991 avt_k (ji,jj,jk) = avt (ji,jj,jk) 992 avm_k (ji,jj,jk) = avm (ji,jj,jk) 993 avmu_k(ji,jj,jk) = avmu(ji,jj,jk) 994 avmv_k(ji,jj,jk) = avmv(ji,jj,jk) 995 END DO 996 END DO 997 END DO 866 998 ! 867 999 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 868 1000 ENDIF 869 1001 ELSE !* Start from rest 870 en(:,:,:) = rn_emin * tmask(:,:,:) 1002 !$OMP PARALLEL 1003 !$OMP DO schedule(static) private(jk,jj,ji) 1004 DO jk = 1, jpk 1005 DO jj = 1, jpj 1006 DO ji = 1, jpi 1007 en(ji,jj,jk) = rn_emin * tmask(ji,jj,jk) 1008 END DO 1009 END DO 1010 END DO 1011 !$OMP END DO NOWAIT 1012 !$OMP DO schedule(static) private(jk) 871 1013 DO jk = 1, jpk ! set the Kz to the background value 872 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 873 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 874 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 875 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 1014 DO jj = 1, jpj 1015 DO ji = 1, jpi 1016 avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk) 1017 avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk) 1018 avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk) 1019 avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk) 1020 END DO 1021 END DO 876 1022 END DO 1023 !$OMP END PARALLEL 877 1024 ENDIF 878 1025 !
Note: See TracChangeset
for help on using the changeset viewer.
