Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4624 r5965 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 IF ( ln_isfcav ) THEN 244 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 245 DO ji = fs_2, fs_jpim1 ! vector opt. 246 en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 247 END DO 248 END DO 249 END IF 241 250 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 251 DO ji = fs_2, fs_jpim1 ! vector opt. … … 291 300 END DO 292 301 ! ! finite LC depth 293 # if defined key_vectopt_loop294 DO jj = 1, 1295 DO ji = 1, jpij ! vector opt. (forced unrolling)296 # else297 302 DO jj = 1, jpj 298 303 DO ji = 1, jpi 299 # endif300 304 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 301 305 END DO … … 313 317 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 314 318 ! ! TKE Langmuir circulation source term 315 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)319 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 316 320 END DO 317 321 END DO … … 332 336 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 333 337 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 334 & / ( fse3uw_n(ji,jj,jk)&335 & * fse3uw_b(ji,jj,jk))338 & / ( fse3uw_n(ji,jj,jk) & 339 & * fse3uw_b(ji,jj,jk) ) 336 340 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 337 341 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & … … 360 364 ! ! right hand side in en 361 365 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 362 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) * tmask(ji,jj,jk) 366 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 367 & * wmask(ji,jj,jk) 363 368 END DO 364 369 END DO … … 372 377 END DO 373 378 END DO 374 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 375 DO ji = fs_2, fs_jpim1 ! vector opt. 379 ! 380 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 381 DO jj = 2, jpjm1 382 DO ji = fs_2, fs_jpim1 ! vector opt. 376 383 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 377 384 END DO … … 384 391 END DO 385 392 END DO 386 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 387 DO ji = fs_2, fs_jpim1 ! vector opt. 393 ! 394 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 395 DO jj = 2, jpjm1 396 DO ji = fs_2, fs_jpim1 ! vector opt. 388 397 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 389 398 END DO … … 399 408 DO jj = 2, jpjm1 400 409 DO ji = fs_2, fs_jpim1 ! vector opt. 401 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)410 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 402 411 END DO 403 412 END DO … … 412 421 DO ji = fs_2, fs_jpim1 ! vector opt. 413 422 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 414 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)423 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 415 424 END DO 416 425 END DO … … 421 430 jk = nmln(ji,jj) 422 431 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 423 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)432 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 424 433 END DO 425 434 END DO … … 433 442 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 434 443 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 435 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! module of the mean stress444 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 436 445 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 437 446 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 438 447 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 439 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)448 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 440 449 END DO 441 450 END DO … … 505 514 ! !* Buoyancy length scale: l=sqrt(2*e/n**2) 506 515 ! 516 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 517 zmxlm(:,:,:) = rmxl_min 518 zmxld(:,:,:) = rmxl_min 519 ! 507 520 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 508 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 509 zmxlm(:,:,1) = MAX( rn_mxl0, zraug * taum(:,:) ) 510 ELSE ! surface set to the minimum value 521 DO jj = 2, jpjm1 522 DO ji = fs_2, fs_jpim1 523 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 524 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 525 END DO 526 END DO 527 ELSE 511 528 zmxlm(:,:,1) = rn_mxl0 512 529 ENDIF 513 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value514 530 ! 515 531 !CDIR NOVERRCHK … … 520 536 DO ji = fs_2, fs_jpim1 ! vector opt. 521 537 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 522 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) 538 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 523 539 END DO 524 540 END DO … … 527 543 ! !* Physical limits for the mixing length 528 544 ! 529 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value545 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 530 546 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 531 547 ! 532 548 SELECT CASE ( nn_mxl ) 533 549 ! 550 ! where wmask = 0 set zmxlm == fse3w 534 551 CASE ( 0 ) ! bounded by the distance to surface and bottom 535 552 DO jk = 2, jpkm1 536 553 DO jj = 2, jpjm1 537 554 DO ji = fs_2, fs_jpim1 ! vector opt. 538 zemxl = MIN( fsdepw(ji,jj,jk) , zmxlm(ji,jj,jk), &555 zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 539 556 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 540 zmxlm(ji,jj,jk) = zemxl 541 zmxld(ji,jj,jk) = zemxl 557 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 558 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 559 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 542 560 END DO 543 561 END DO … … 620 638 zsqen = SQRT( en(ji,jj,jk) ) 621 639 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 622 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)623 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)640 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 641 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 624 642 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 625 643 END DO … … 628 646 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 629 647 ! 630 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- andv-points648 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 631 649 DO jj = 2, jpjm1 632 650 DO ji = fs_2, fs_jpim1 ! vector opt. 633 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk)634 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk)651 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 652 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 635 653 END DO 636 654 END DO … … 653 671 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 654 672 !!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) ) 655 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)673 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 656 674 # if defined key_c1d 657 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)! c1d configuration : save masked Prandlt number658 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)! c1d config. : save Ri675 e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 676 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 659 677 # endif 660 678 END DO … … 740 758 ! 741 759 ! !* Check of some namelist values 742 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 743 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 744 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 745 #if ! key_coupled 746 IF( nn_etau == 3 ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 747 #endif 760 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 761 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 762 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 763 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 748 764 749 765 IF( ln_mxl0 ) THEN … … 765 781 ! !* set vertical eddy coef. to the background value 766 782 DO jk = 1, jpk 767 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)768 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)769 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)770 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)783 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 784 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 785 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 786 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 771 787 END DO 772 788 dissl(:,:,:) = 1.e-12_wp … … 819 835 en (:,:,:) = rn_emin * tmask(:,:,:) 820 836 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 837 ! 838 avt_k (:,:,:) = avt (:,:,:) 839 avm_k (:,:,:) = avm (:,:,:) 840 avmu_k(:,:,:) = avmu(:,:,:) 841 avmv_k(:,:,:) = avmv(:,:,:) 842 ! 821 843 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 822 844 ENDIF … … 824 846 en(:,:,:) = rn_emin * tmask(:,:,:) 825 847 DO jk = 1, jpk ! set the Kz to the background value 826 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)827 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)828 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)829 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)848 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 849 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 850 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 851 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 830 852 END DO 831 853 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
