- Timestamp:
- 2015-02-20T20:11:47+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4990 r5098 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 28 29 !!---------------------------------------------------------------------- 29 30 #if defined key_zdftke || defined key_esopa … … 236 237 zfact3 = 0.5_wp * rn_ediss 237 238 ! 239 ! 238 240 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 239 241 ! ! Surface boundary condition on tke 240 242 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 243 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 244 DO ji = fs_2, fs_jpim1 ! vector opt. 245 en(ji,jj,mikt(ji,jj))=rn_emin 246 END DO 247 END DO 241 248 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 249 DO ji = fs_2, fs_jpim1 ! vector opt. 243 IF (mikt(ji,jj) .GT. 1) THEN 244 en(ji,jj,mikt(ji,jj))=rn_emin 245 ELSE 246 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 247 END IF 250 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 248 251 END DO 249 252 END DO … … 301 304 END DO 302 305 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 306 !CDIR NOVERRCHK 303 307 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 304 DO jj = 2, jpjm1 308 !CDIR NOVERRCHK 309 DO jj = 2, jpjm1 310 !CDIR NOVERRCHK 305 311 DO ji = fs_2, fs_jpim1 ! vector opt. 306 312 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift … … 309 315 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 310 316 ! ! TKE Langmuir circulation source term 311 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)317 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 312 318 END DO 313 319 END DO … … 338 344 END DO 339 345 ! 340 DO j j = 2, jpjm1341 DO j i = fs_2, fs_jpim1 ! vector opt.342 DO j k = mikt(ji,jj)+1, jpkm1 !* Matrix and right hand side in en346 DO jk = 2, jpkm1 !* Matrix and right hand side in en 347 DO jj = 2, jpjm1 348 DO ji = fs_2, fs_jpim1 ! vector opt. 343 349 zcof = zfact1 * tmask(ji,jj,jk) 344 350 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal … … 357 363 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 358 364 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 359 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 360 END DO 361 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 362 DO jk = mikt(ji,jj)+2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 365 & * wmask(ji,jj,jk) 366 END DO 367 END DO 368 END DO 369 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 370 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 371 DO jj = 2, jpjm1 372 DO ji = fs_2, fs_jpim1 ! vector opt. 363 373 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 364 374 END DO 365 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 366 zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj)) ! Surface boudary conditions on tke 367 ! 368 DO jk = mikt(ji,jj)+2, jpkm1 375 END DO 376 END DO 377 ! 378 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 379 DO jj = 2, jpjm1 380 DO ji = fs_2, fs_jpim1 ! vector opt. 381 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 382 END DO 383 END DO 384 DO jk = 3, jpkm1 385 DO jj = 2, jpjm1 386 DO ji = fs_2, fs_jpim1 ! vector opt. 369 387 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) 370 388 END DO 371 ! 372 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 389 END DO 390 END DO 391 ! 392 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 393 DO jj = 2, jpjm1 394 DO ji = fs_2, fs_jpim1 ! vector opt. 373 395 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 374 ! 375 DO jk = jpk-2, mikt(ji,jj)+1, -1 396 END DO 397 END DO 398 DO jk = jpk-2, 2, -1 399 DO jj = 2, jpjm1 400 DO ji = fs_2, fs_jpim1 ! vector opt. 376 401 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 377 402 END DO 378 ! 379 DO jk = mikt(ji,jj), jpkm1 ! set the minimum value of tke 380 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 403 END DO 404 END DO 405 DO jk = 2, jpkm1 ! set the minimum value of tke 406 DO jj = 2, jpjm1 407 DO ji = fs_2, fs_jpim1 ! vector opt. 408 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 381 409 END DO 382 410 END DO … … 391 419 DO ji = fs_2, fs_jpim1 ! vector opt. 392 420 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 393 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1)421 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 394 422 END DO 395 423 END DO … … 400 428 jk = nmln(ji,jj) 401 429 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 402 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1)430 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 403 431 END DO 404 432 END DO … … 416 444 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 417 445 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 418 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1)446 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 419 447 END DO 420 448 END DO … … 484 512 ! !* Buoyancy length scale: l=sqrt(2*e/n**2) 485 513 ! 514 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 515 zmxlm(:,:,:) = rmxl_min 516 zmxld(:,:,:) = rmxl_min 517 ! 486 518 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 487 519 DO jj = 2, jpjm1 488 520 DO ji = fs_2, fs_jpim1 489 IF (mikt(ji,jj) .GT. 1) THEN 490 zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 491 ELSE 492 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 493 zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 494 END IF 521 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 522 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 495 523 END DO 496 524 END DO 497 525 ELSE 498 DO jj = 2, jpjm1 499 DO ji = fs_2, fs_jpim1 ! surface set to the minimum value 500 zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 501 END DO 502 END DO 526 zmxlm(:,:,1) = rn_mxl0 503 527 ENDIF 504 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 505 ! 506 !CDIR NOVERRCHK 507 DO jj = 2, jpjm1 508 !CDIR NOVERRCHK 509 DO ji = fs_2, fs_jpim1 ! vector opt. 510 !CDIR NOVERRCHK 511 DO jk = mikt(ji,jj)+1, jpkm1 ! interior value : l=sqrt(2*e/n^2) 528 ! 529 !CDIR NOVERRCHK 530 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 531 !CDIR NOVERRCHK 532 DO jj = 2, jpjm1 533 !CDIR NOVERRCHK 534 DO ji = fs_2, fs_jpim1 ! vector opt. 512 535 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 513 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 514 END DO 515 zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj)) ! surface set to the minimum value 536 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 537 END DO 516 538 END DO 517 539 END DO … … 519 541 ! !* Physical limits for the mixing length 520 542 ! 521 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value543 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 522 544 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 523 545 ! 524 546 SELECT CASE ( nn_mxl ) 525 547 ! 548 ! where wmask = 0 set zmxlm == fse3w 526 549 CASE ( 0 ) ! bounded by the distance to surface and bottom 527 DO j j = 2, jpjm1528 DO j i = fs_2, fs_jpim1 ! vector opt.529 DO j k = mikt(ji,jj)+1, jpkm1550 DO jk = 2, jpkm1 551 DO jj = 2, jpjm1 552 DO ji = fs_2, fs_jpim1 ! vector opt. 530 553 zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 531 554 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 532 zmxlm(ji,jj,jk) = zemxl 533 zmxld(ji,jj,jk) = zemxl 555 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 556 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 557 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 534 558 END DO 535 559 END DO … … 537 561 ! 538 562 CASE ( 1 ) ! bounded by the vertical scale factor 539 DO j j = 2, jpjm1540 DO j i = fs_2, fs_jpim1 ! vector opt.541 DO j k = mikt(ji,jj)+1, jpkm1563 DO jk = 2, jpkm1 564 DO jj = 2, jpjm1 565 DO ji = fs_2, fs_jpim1 ! vector opt. 542 566 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 543 567 zmxlm(ji,jj,jk) = zemxl … … 548 572 ! 549 573 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 550 DO j j = 2, jpjm1551 DO j i = fs_2, fs_jpim1 ! vector opt.552 DO j k = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom :574 DO jk = 2, jpkm1 ! from the surface to the bottom : 575 DO jj = 2, jpjm1 576 DO ji = fs_2, fs_jpim1 ! vector opt. 553 577 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 554 578 END DO 555 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : 579 END DO 580 END DO 581 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 582 DO jj = 2, jpjm1 583 DO ji = fs_2, fs_jpim1 ! vector opt. 556 584 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 557 585 zmxlm(ji,jj,jk) = zemxl … … 562 590 ! 563 591 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 564 DO j j = 2, jpjm1565 DO j i = fs_2, fs_jpim1 ! vector opt.566 DO j k = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : lup592 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 593 DO jj = 2, jpjm1 594 DO ji = fs_2, fs_jpim1 ! vector opt. 567 595 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 568 596 END DO 569 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : ldown 597 END DO 598 END DO 599 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 600 DO jj = 2, jpjm1 601 DO ji = fs_2, fs_jpim1 ! vector opt. 570 602 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 571 603 END DO … … 604 636 zsqen = SQRT( en(ji,jj,jk) ) 605 637 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 606 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)607 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)638 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 639 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 608 640 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 609 641 END DO … … 612 644 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 613 645 ! 614 DO jj = 2, jpjm1 615 DO ji = fs_2, fs_jpim1 ! vector opt. 616 DO jk = miku(ji,jj)+1, jpkm1 !* vertical eddy viscosity at u- and v-points 617 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 618 END DO 619 DO jk = mikv(ji,jj)+1, jpkm1 620 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 646 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 647 DO jj = 2, jpjm1 648 DO ji = fs_2, fs_jpim1 ! vector opt. 649 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 650 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 621 651 END DO 622 652 END DO … … 625 655 ! 626 656 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 627 DO j j = 2, jpjm1628 DO j i = fs_2, fs_jpim1 ! vector opt.629 DO j k = mikt(ji,jj)+1, jpkm1657 DO jk = 2, jpkm1 658 DO jj = 2, jpjm1 659 DO ji = fs_2, fs_jpim1 ! vector opt. 630 660 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 631 661 ! ! shear … … 639 669 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 640 670 !!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) ) 641 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)671 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 642 672 # if defined key_c1d 643 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number644 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri673 e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 674 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 645 675 # endif 646 676 END DO … … 749 779 ! !* set vertical eddy coef. to the background value 750 780 DO jk = 1, jpk 751 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)752 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)753 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)754 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)781 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 782 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 783 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 784 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 755 785 END DO 756 786 dissl(:,:,:) = 1.e-12_wp … … 808 838 en(:,:,:) = rn_emin * tmask(:,:,:) 809 839 DO jk = 1, jpk ! set the Kz to the background value 810 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)811 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)812 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)813 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)840 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 841 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 842 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 843 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 814 844 END DO 815 845 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
