Changeset 4990 for trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2014-12-15T17:42:49+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4624 r4990 241 241 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 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 244 248 END DO 245 249 END DO … … 291 295 END DO 292 296 ! ! finite LC depth 293 # if defined key_vectopt_loop294 DO jj = 1, 1295 DO ji = 1, jpij ! vector opt. (forced unrolling)296 # else297 297 DO jj = 1, jpj 298 298 DO ji = 1, jpi 299 # endif300 299 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 301 300 END DO 302 301 END DO 303 302 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 304 !CDIR NOVERRCHK305 303 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 306 !CDIR NOVERRCHK307 304 DO jj = 2, jpjm1 308 !CDIR NOVERRCHK309 305 DO ji = fs_2, fs_jpim1 ! vector opt. 310 306 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift … … 342 338 END DO 343 339 ! 344 DO j k = 2, jpkm1 !* Matrix and right hand side in en345 DO j j = 2, jpjm1346 DO j i = fs_2, fs_jpim1 ! vector opt.340 DO jj = 2, jpjm1 341 DO ji = fs_2, fs_jpim1 ! vector opt. 342 DO jk = mikt(ji,jj)+1, jpkm1 !* Matrix and right hand side in en 347 343 zcof = zfact1 * tmask(ji,jj,jk) 348 344 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal … … 360 356 ! ! right hand side in en 361 357 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) 363 END DO 364 END DO 365 END DO 366 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 367 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 368 DO jj = 2, jpjm1 369 DO ji = fs_2, fs_jpim1 ! vector opt. 358 & + 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 370 363 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 371 364 END DO 372 END DO 373 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. 376 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 377 END DO 378 END DO 379 DO jk = 3, jpkm1 380 DO jj = 2, jpjm1 381 DO ji = fs_2, fs_jpim1 ! vector opt. 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 382 369 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) 383 370 END DO 384 END DO 385 END DO 386 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 387 DO ji = fs_2, fs_jpim1 ! vector opt. 371 ! 372 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 388 373 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 389 END DO 390 END DO 391 DO jk = jpk-2, 2, -1 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 ! vector opt. 374 ! 375 DO jk = jpk-2, mikt(ji,jj)+1, -1 394 376 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 395 377 END DO 396 END DO 397 END DO 398 DO jk = 2, jpkm1 ! set the minimum value of tke 399 DO jj = 2, jpjm1 400 DO ji = fs_2, fs_jpim1 ! vector opt. 378 ! 379 DO jk = mikt(ji,jj), jpkm1 ! set the minimum value of tke 401 380 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 402 381 END DO … … 412 391 DO ji = fs_2, fs_jpim1 ! vector opt. 413 392 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)393 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1) 415 394 END DO 416 395 END DO … … 421 400 jk = nmln(ji,jj) 422 401 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)402 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1) 424 403 END DO 425 404 END DO … … 433 412 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 434 413 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 435 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! module of the mean stress414 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 436 415 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 437 416 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 438 417 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)418 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 440 419 END DO 441 420 END DO … … 506 485 ! 507 486 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 511 zmxlm(:,:,1) = rn_mxl0 487 DO jj = 2, jpjm1 488 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 495 END DO 496 END DO 497 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 512 503 ENDIF 513 504 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 514 505 ! 515 506 !CDIR NOVERRCHK 516 DO j k = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2)517 !CDIR NOVERRCHK 518 DO j j = 2, jpjm1519 !CDIR NOVERRCHK520 DO j i = fs_2, fs_jpim1 ! vector opt.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) 521 512 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 522 513 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 523 514 END DO 515 zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj)) ! surface set to the minimum value 524 516 END DO 525 517 END DO … … 533 525 ! 534 526 CASE ( 0 ) ! bounded by the distance to surface and bottom 535 DO j k = 2, jpkm1536 DO j j = 2, jpjm1537 DO j i = fs_2, fs_jpim1 ! vector opt.538 zemxl = MIN( fsdepw(ji,jj,jk) , zmxlm(ji,jj,jk), &527 DO jj = 2, jpjm1 528 DO ji = fs_2, fs_jpim1 ! vector opt. 529 DO jk = mikt(ji,jj)+1, jpkm1 530 zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 539 531 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 540 532 zmxlm(ji,jj,jk) = zemxl … … 545 537 ! 546 538 CASE ( 1 ) ! bounded by the vertical scale factor 547 DO j k = 2, jpkm1548 DO j j = 2, jpjm1549 DO j i = fs_2, fs_jpim1 ! vector opt.539 DO jj = 2, jpjm1 540 DO ji = fs_2, fs_jpim1 ! vector opt. 541 DO jk = mikt(ji,jj)+1, jpkm1 550 542 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 551 543 zmxlm(ji,jj,jk) = zemxl … … 556 548 ! 557 549 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 558 DO j k = 2, jpkm1 ! from the surface to the bottom :559 DO j j = 2, jpjm1560 DO j i = fs_2, fs_jpim1 ! vector opt.550 DO jj = 2, jpjm1 551 DO ji = fs_2, fs_jpim1 ! vector opt. 552 DO jk = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : 561 553 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 562 554 END DO 563 END DO 564 END DO 565 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 566 DO jj = 2, jpjm1 567 DO ji = fs_2, fs_jpim1 ! vector opt. 555 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : 568 556 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 569 557 zmxlm(ji,jj,jk) = zemxl … … 574 562 ! 575 563 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 576 DO j k = 2, jpkm1 ! from the surface to the bottom : lup577 DO j j = 2, jpjm1578 DO j i = fs_2, fs_jpim1 ! vector opt.564 DO jj = 2, jpjm1 565 DO ji = fs_2, fs_jpim1 ! vector opt. 566 DO jk = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : lup 579 567 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 580 568 END DO 581 END DO 582 END DO 583 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 569 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : ldown 586 570 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 587 571 END DO … … 628 612 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 629 613 ! 630 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 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) 621 END DO 622 END DO 623 END DO 624 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 625 ! 626 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 631 627 DO jj = 2, jpjm1 632 628 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) 635 END DO 636 END DO 637 END DO 638 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 639 ! 640 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 641 DO jk = 2, jpkm1 642 DO jj = 2, jpjm1 643 DO ji = fs_2, fs_jpim1 ! vector opt. 629 DO jk = mikt(ji,jj)+1, jpkm1 644 630 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 645 631 ! ! shear … … 655 641 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 656 642 # if defined key_c1d 657 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) 658 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)! c1d config. : save Ri643 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 644 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri 659 645 # endif 660 646 END DO … … 740 726 ! 741 727 ! !* 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 728 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 729 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 730 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 731 IF( nn_etau == 3 .AND. .NOT. lk_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 748 732 749 733 IF( ln_mxl0 ) THEN
Note: See TracChangeset
for help on using the changeset viewer.