- Timestamp:
- 2011-12-11T16:00:26+01:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r2715 r3211 82 82 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 83 83 84 !! DCSE_NEMO: en is public because it is used by asmtrj 84 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] 85 86 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 86 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing leng htof dissipation87 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing length of dissipation 87 88 #if defined key_c1d 88 89 ! !!** 1D cfg only ** ('key_c1d') 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers 91 #endif 90 !! DCSE_NEMO: these arrays do not need to be public 91 ! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent length scales 92 ! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers 93 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent length scales 94 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers 95 #endif 96 97 !! * Control permutation of array indices 98 # include "zdftke_ftrans.h90" 99 # include "oce_ftrans.h90" 100 # include "dom_oce_ftrans.h90" 101 # include "domvvl_ftrans.h90" 102 # include "sbc_oce_ftrans.h90" 103 # include "zdf_oce_ftrans.h90" 104 !FTRANS dissl e_dis e_mix e_pdl e_ric :I :I :z 92 105 93 106 !! * Substitutions … … 195 208 USE wrk_nemo, ONLY: zhlc => wrk_2d_1 ! 2D REAL workspace 196 209 USE wrk_nemo, ONLY: zpelc => wrk_3d_1 ! 3D REAL workspace 210 211 !! DCSE_NEMO: need additional directives for renamed module variables 212 !FTRANS zdiag zd_up zd_lw :I :I :z 213 !FTRANS zpelc :I :I :z 197 214 ! 198 215 INTEGER :: ji, jj, jk ! dummy loop arguments … … 226 243 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 227 244 DO ji = fs_2, fs_jpim1 ! vector opt. 245 #if defined key_z_first 246 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask_1(ji,jj) 247 #else 228 248 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 249 #endif 229 250 END DO 230 251 END DO … … 260 281 ! 261 282 ! !* total energy produce by LC : cumulative sum over jk 283 #if defined key_z_first 284 DO jj = 1, jpj 285 DO ji = 1, jpi 286 zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * fsdepw(ji,jj,1) * fse3w(ji,jj,1) 287 DO jk = 2, jpk 288 zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + MAX( rn2b(ji,jj,jk), 0._wp ) * fsdepw(ji,jj,jk) * fse3w(ji,jj,jk) 289 END DO 290 END DO 291 END DO 292 #else 262 293 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1) 263 294 DO jk = 2, jpk 264 295 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 265 296 END DO 297 #endif 266 298 ! !* finite Langmuir Circulation depth 267 299 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 268 300 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 301 #if defined key_z_first 302 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 303 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 304 zus = zcof * taum(ji,jj) 305 DO jk = jpkm1, 2, -1 306 #else 269 307 DO jk = jpkm1, 2, -1 270 308 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 271 309 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 272 310 zus = zcof * taum(ji,jj) 311 #endif 273 312 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 274 313 END DO … … 287 326 END DO 288 327 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 328 #if defined key_z_first 329 DO jj = 2, jpjm1 !* TKE Langmuir circulation source term added to en 330 DO ji = 2, jpim1 331 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 332 DO jk = 2, jpkm1 333 #else 289 334 !CDIR NOVERRCHK 290 335 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en … … 294 339 DO ji = fs_2, fs_jpim1 ! vector opt. 295 340 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 341 #endif 296 342 ! ! vertical velocity due to LC 297 343 zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) … … 312 358 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 313 359 ! 360 #if defined key_z_first 361 !* Shear production at uw- and vw-points (energy conserving form) 362 ! here avmu, avmv used as workspace 363 DO jj = 1, jpj 364 DO ji = 1, jpi 365 DO jk = 2, jpkm1 366 #else 314 367 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 315 368 DO jj = 1, jpj ! here avmu, avmv used as workspace 316 369 DO ji = 1, jpi 370 #endif 317 371 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 318 372 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & … … 326 380 END DO 327 381 END DO 328 ! 382 383 ! 384 #if defined key_z_first 385 DO jj = 2, jpjm1 386 DO ji = 2, jpim1 387 DO jk = 2, jpkm1 !* Matrix and right hand side in en 388 #else 329 389 DO jk = 2, jpkm1 !* Matrix and right hand side in en 330 390 DO jj = 2, jpjm1 331 391 DO ji = fs_2, fs_jpim1 ! vector opt. 392 #endif 332 393 zcof = zfact1 * tmask(ji,jj,jk) 333 394 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal … … 350 411 END DO 351 412 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 413 #if defined key_z_first 414 DO jj = 2, jpjm1 415 DO ji = 2, jpim1 416 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 417 DO jk = 3, jpkm1 418 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 419 END DO 420 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 421 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 422 DO jk = 3, jpkm1 423 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 424 END DO 425 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 426 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 427 DO jk = jpk-2, 2, -1 428 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 429 END DO 430 DO jk = 2, jpkm1 ! set the minimum value of tke 431 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 432 END DO 433 END DO 434 END DO 435 #else 352 436 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 353 437 DO jj = 2, jpjm1 … … 388 472 END DO 389 473 END DO 474 #endif 390 475 391 476 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 393 478 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 394 479 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 480 #if defined key_z_first 481 DO jj = 2, jpjm1 482 DO ji = 2, jpim1 483 DO jk = 2, jpkm1 484 #else 395 485 DO jk = 2, jpkm1 396 486 DO jj = 2, jpjm1 397 487 DO ji = fs_2, fs_jpim1 ! vector opt. 488 #endif 398 489 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 399 490 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) … … 410 501 END DO 411 502 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 503 504 !! DCSE_NEMO: its probably not worth changing the order of these loops for level first indexing, 505 !! unless we also make zdif a 2-d (jpi,jpj) array 412 506 !CDIR NOVERRCHK 413 507 DO jk = 2, jpkm1 … … 435 529 END SUBROUTINE tke_tke 436 530 531 !! * Reset control of array index permutation 532 # include "zdftke_ftrans.h90" 533 # include "oce_ftrans.h90" 534 # include "dom_oce_ftrans.h90" 535 # include "domvvl_ftrans.h90" 536 # include "sbc_oce_ftrans.h90" 537 # include "zdf_oce_ftrans.h90" 538 !FTRANS dissl e_dis e_mix e_pdl e_ric :I :I :z 437 539 438 540 SUBROUTINE tke_avn … … 472 574 !!---------------------------------------------------------------------- 473 575 USE oce, ONLY: zmpdl => ua , zmxlm => va , zmxld => ta ! (ua,va,ta) used as workspace 576 !! DCSE_NEMO: need additional directives for renamed module variables 577 !FTRANS zmpdl zmxlm zmxld :I :I :z 474 578 !! 475 579 INTEGER :: ji, jj, jk ! dummy loop indices … … 491 595 zmxlm(:,:,1) = rn_mxl0 492 596 ENDIF 493 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 597 598 #if defined key_z_first 599 DO jj = 2, jpjm1 600 DO ji = 2, jpim1 601 zmxlm(ji,jj,jpk) = rmxl_min ! last level set to the interior minium value 602 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 603 #else 604 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 494 605 ! 495 606 !CDIR NOVERRCHK … … 499 610 !CDIR NOVERRCHK 500 611 DO ji = fs_2, fs_jpim1 ! vector opt. 612 #endif 501 613 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 502 614 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) … … 513 625 ! 514 626 CASE ( 0 ) ! bounded by the distance to surface and bottom 627 #if defined key_z_first 628 DO jj = 2, jpjm1 629 DO ji = 2, jpim1 630 DO jk = 2, jpkm1 631 #else 515 632 DO jk = 2, jpkm1 516 633 DO jj = 2, jpjm1 517 634 DO ji = fs_2, fs_jpim1 ! vector opt. 635 #endif 518 636 zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk), & 519 637 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) … … 525 643 ! 526 644 CASE ( 1 ) ! bounded by the vertical scale factor 645 #if defined key_z_first 646 DO jj = 2, jpjm1 647 DO ji = 2, jpim1 648 DO jk = 2, jpkm1 649 #else 527 650 DO jk = 2, jpkm1 528 651 DO jj = 2, jpjm1 529 652 DO ji = fs_2, fs_jpim1 ! vector opt. 653 #endif 530 654 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 531 655 zmxlm(ji,jj,jk) = zemxl … … 536 660 ! 537 661 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 662 #if defined key_z_first 663 DO jj = 2, jpjm1 664 DO ji = 2, jpim1 665 DO jk = 2, jpkm1 ! from the surface to the bottom : 666 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 667 END DO 668 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 669 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 670 zmxlm(ji,jj,jk) = zemxl 671 zmxld(ji,jj,jk) = zemxl 672 END DO 673 END DO 674 END DO 675 #else 538 676 DO jk = 2, jpkm1 ! from the surface to the bottom : 539 677 DO jj = 2, jpjm1 … … 552 690 END DO 553 691 END DO 692 #endif 554 693 ! 555 694 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 695 #if defined key_z_first 696 DO jj = 2, jpjm1 697 DO ji = 2, jpim1 698 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 699 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 700 END DO 701 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 702 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 703 END DO 704 DO jk = 2, jpkm1 705 zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) ) 706 zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) ) 707 zmxlm(ji,jj,jk) = zemlm 708 zmxld(ji,jj,jk) = zemlp 709 END DO 710 END DO 711 END DO 712 #else 556 713 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 557 714 DO jj = 2, jpjm1 … … 581 738 END DO 582 739 END DO 740 #endif 583 741 ! 584 742 END SELECT … … 592 750 ! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt) 593 751 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 752 #if defined key_z_first 753 DO jj = 2, jpjm1 754 DO ji = 2, jpim1 755 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points 756 #else 594 757 !CDIR NOVERRCHK 595 758 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points … … 598 761 !CDIR NOVERRCHK 599 762 DO ji = fs_2, fs_jpim1 ! vector opt. 763 #endif 600 764 zsqen = SQRT( en(ji,jj,jk) ) 601 765 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen … … 608 772 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 609 773 ! 774 #if defined key_z_first 775 DO jj = 2, jpjm1 776 DO ji = 2, jpim1 777 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 778 #else 610 779 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 611 780 DO jj = 2, jpjm1 612 781 DO ji = fs_2, fs_jpim1 ! vector opt. 782 #endif 613 783 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 614 784 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) … … 619 789 ! 620 790 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 791 #if defined key_z_first 792 DO jj = 2, jpjm1 793 DO ji = 2, jpim1 794 DO jk = 2, jpkm1 795 #else 621 796 DO jk = 2, jpkm1 622 797 DO jj = 2, jpjm1 623 798 DO ji = fs_2, fs_jpim1 ! vector opt. 799 #endif 624 800 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 625 801 ! ! shear … … 652 828 END SUBROUTINE tke_avn 653 829 830 !! * Reset control of array index permutation 831 # include "zdftke_ftrans.h90" 832 # include "oce_ftrans.h90" 833 # include "dom_oce_ftrans.h90" 834 # include "domvvl_ftrans.h90" 835 # include "sbc_oce_ftrans.h90" 836 # include "zdf_oce_ftrans.h90" 837 !FTRANS dissl e_dis e_mix e_pdl e_ric :I :I :z 654 838 655 839 SUBROUTINE zdf_tke_init … … 733 917 ENDIF 734 918 ! !* set vertical eddy coef. to the background value 919 #if defined key_z_first 920 DO jj = 1, jpj 921 DO ji = 1, jpi 922 avt (ji,jj,:) = avtb(:) * tmask(ji,jj,:) 923 avm (ji,jj,:) = avmb(:) * tmask(ji,jj,:) 924 avmu(ji,jj,:) = avmb(:) * umask(ji,jj,:) 925 avmv(ji,jj,:) = avmb(:) * vmask(ji,jj,:) 926 END DO 927 END DO 928 #else 735 929 DO jk = 1, jpk 736 930 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) … … 739 933 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 740 934 END DO 935 #endif 741 936 dissl(:,:,:) = 1.e-12_wp 742 937 ! … … 759 954 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 760 955 ! 761 INTEGER :: jit, j k ! dummy loop indices956 INTEGER :: jit, ji, jj, jk ! dummy loop indices 762 957 INTEGER :: id1, id2, id3, id4, id5, id6 ! local integers 763 958 !!---------------------------------------------------------------------- … … 792 987 ELSE !* Start from rest 793 988 en(:,:,:) = rn_emin * tmask(:,:,:) 989 #if defined key_z_first 990 DO jj = 1, jpj ! set the Kz to the background value 991 DO ji = 1, jpi 992 avt (ji,jj,:) = avtb(:) * tmask(ji,jj,:) 993 avm (ji,jj,:) = avmb(:) * tmask(ji,jj,:) 994 avmu(ji,jj,:) = avmb(:) * umask(ji,jj,:) 995 avmv(ji,jj,:) = avmb(:) * vmask(ji,jj,:) 996 END DO 997 END DO 998 #else 794 999 DO jk = 1, jpk ! set the Kz to the background value 795 1000 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) … … 798 1003 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 799 1004 END DO 1005 #endif 1006 800 1007 ENDIF 801 1008 !
Note: See TracChangeset
for help on using the changeset viewer.
