Changeset 12448 for branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
- Timestamp:
- 2020-02-24T17:39:27+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package_collate_spectral_optics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r8058 r12448 24 24 USE phycst ! physical constants 25 25 USE zdfmxl ! mixed layer 26 USE sbcwave 27 ! 26 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 29 USE lib_mpp ! MPP manager … … 46 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 47 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points 50 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_tke1 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_tke3 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_psi1 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rtrans 48 55 49 56 ! !! ** Namelist namzdf_gls ** … … 59 66 REAL(wp) :: rn_emin ! minimum value of TKE (m2/s2) 60 67 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 61 REAL(wp) :: rn_crban 68 REAL(wp) :: rn_crban_default ! Craig and Banner constant for surface breaking waves mixing 62 69 REAL(wp) :: rn_hsro ! Minimum surface roughness 63 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1)70 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met = 1, 3) 64 71 65 72 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters … … 90 97 REAL(wp) :: rm7 = 0.0_wp 91 98 REAL(wp) :: rm8 = 0.318_wp 92 REAL(wp) :: rtrans 99 REAL(wp) :: rtrans_default = 0.1_wp 93 100 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 94 REAL(wp) :: rsbc_tke1 , rsbc_tke2, rfact_tke ! - - - -95 REAL(wp) :: rsbc_psi 1, rsbc_psi2, rfact_psi ! - - - -101 REAL(wp) :: rsbc_tke1_default, rsbc_tke2, rfact_tke ! - - - - 102 REAL(wp) :: rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 96 103 REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - - 97 104 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - … … 116 123 !!---------------------------------------------------------------------- 117 124 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 118 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 125 & rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) , & 126 & rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), & 127 & ustars2(jpi,jpj), ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 119 128 ! 120 129 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 161 170 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 162 171 172 ! variable initialization 173 IF( ln_phioc ) THEN 174 IF( nn_wmix==jp_janssen ) THEN 175 rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53.Dirichlet + Wave breaking 176 rsbc_tke3(:,:) = rdt * rn_crban(:,:) ! Neumann + Wave breaking 177 rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 178 ELSE IF( nn_wmix==jp_craigbanner ) THEN 179 rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf 180 rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf 181 rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 182 ENDIF 183 ENDIF 184 163 185 IF( kt /= nit000 ) THEN ! restore before value to compute tke 164 186 avt (:,:,:) = avt_k (:,:,:) … … 197 219 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 198 220 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 199 ! 221 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 222 zhsro = MAX(rn_frac_hs * hsw, rn_hsro) 200 223 END SELECT 201 224 … … 311 334 ! --------------------------------- 312 335 ! 336 IF( ln_phioc ) THEN 337 SELECT CASE ( nn_bc_surf ) 338 ! 339 CASE ( 0 ) ! Dirichlet case 340 IF( nn_wmix==jp_janssen ) THEN 341 ! First level 342 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 343 z_elem_a(:,:,1) = en(:,:,1) 344 z_elem_c(:,:,1) = 0._wp 345 z_elem_b(:,:,1) = 1._wp 346 347 ! One level below 348 en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 349 z_elem_a(:,:,2) = 0._wp 350 z_elem_c(:,:,2) = 0._wp 351 z_elem_b(:,:,2) = 1._wp 352 ELSE IF( nn_wmix==jp_craigbanner ) THEN 353 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 354 en(:,:,1) = MAX(en(:,:,1), rn_emin) 355 z_elem_a(:,:,1) = en(:,:,1) 356 z_elem_c(:,:,1) = 0._wp 357 z_elem_b(:,:,1) = 1._wp 358 ! 359 ! One level below 360 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) & 361 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 362 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 363 z_elem_a(:,:,2) = 0._wp 364 z_elem_c(:,:,2) = 0._wp 365 z_elem_b(:,:,2) = 1._wp 366 ! 367 ENDIF 368 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 369 IF( nn_wmix==jp_janssen ) THEN 370 ! Dirichlet conditions at k=1 371 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 372 z_elem_a(:,:,1) = en(:,:,1) 373 z_elem_c(:,:,1) = 0._wp 374 z_elem_b(:,:,1) = 1._wp 375 ! at k=2, set de/dz=Fw 376 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 377 z_elem_a(:,:,2) = 0._wp 378 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 379 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 380 ELSE IF( nn_wmix==jp_craigbanner ) THEN 381 ! Dirichlet conditions at k=1 382 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 383 en(:,:,1) = MAX(en(:,:,1), rn_emin) 384 z_elem_a(:,:,1) = en(:,:,1) 385 z_elem_c(:,:,1) = 0._wp 386 z_elem_b(:,:,1) = 1._wp 387 ! 388 ! at k=2, set de/dz=Fw 389 !cbr 390 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 391 z_elem_a(:,:,2) = 0._wp 392 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) )) 393 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) & 394 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 395 396 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 397 ! 398 ENDIF 399 END SELECT 400 ELSE 313 401 SELECT CASE ( nn_bc_surf ) 314 402 ! 315 403 CASE ( 0 ) ! Dirichlet case 316 404 ! First level 317 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 )**(2._wp/3._wp)405 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 318 406 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 407 z_elem_a(:,:,1) = en(:,:,1) … … 322 410 ! 323 411 ! One level below 324 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 325 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 412 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 326 413 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 327 414 z_elem_a(:,:,2) = 0._wp … … 331 418 ! 332 419 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 333 !334 420 ! Dirichlet conditions at k=1 335 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 )**(2._wp/3._wp)421 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 336 422 en(:,:,1) = MAX(en(:,:,1), rn_emin) 337 423 z_elem_a(:,:,1) = en(:,:,1) … … 343 429 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 344 430 z_elem_a(:,:,2) = 0._wp 345 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 346 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 347 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 431 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) )) 432 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 348 433 349 434 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) … … 351 436 ! 352 437 END SELECT 438 ENDIF ! ln_phioc 353 439 354 440 ! Bottom boundary condition on tke … … 536 622 ! 537 623 CASE ( 0 ) ! Dirichlet boundary conditions 538 ! 539 ! Surface value 540 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 541 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 542 z_elem_a(:,:,1) = psi(:,:,1) 543 z_elem_c(:,:,1) = 0._wp 544 z_elem_b(:,:,1) = 1._wp 545 ! 546 ! One level below 547 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 548 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 549 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 550 z_elem_a(:,:,2) = 0._wp 551 z_elem_c(:,:,2) = 0._wp 552 z_elem_b(:,:,2) = 1._wp 553 ! 554 ! 624 IF( ln_phioc ) THEN ! Wave induced mixing case 625 ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 626 ! balance between the production and the 627 ! dissipation terms including the wave effect 628 IF( nn_wmix==jp_janssen ) THEN 629 ! Surface value 630 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 631 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 632 z_elem_a(:,:,1) = psi(:,:,1) 633 z_elem_c(:,:,1) = 0._wp 634 z_elem_b(:,:,1) = 1._wp 635 ! 636 ! One level below 637 zex1 = (rmm*ra_sf+rnn) 638 zex2 = (rmm*ra_sf) 639 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 640 psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 641 z_elem_a(:,:,2) = 0._wp 642 z_elem_c(:,:,2) = 0._wp 643 z_elem_b(:,:,2) = 1._wp 644 ! 645 ! 646 ELSE IF( nn_wmix==jp_craigbanner ) THEN 647 ! Surface value 648 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 649 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 650 z_elem_a(:,:,1) = psi(:,:,1) 651 z_elem_c(:,:,1) = 0._wp 652 z_elem_b(:,:,1) = 1._wp 653 ! 654 ! One level below 655 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) ))) 656 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 657 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 658 z_elem_a(:,:,2) = 0._wp 659 z_elem_c(:,:,2) = 0._wp 660 z_elem_b(:,:,2) = 1._wp 661 ! 662 ! 663 ENDIF 664 ELSE ! Wave induced mixing case with default values 665 ! Surface value 666 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 667 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 668 z_elem_a(:,:,1) = psi(:,:,1) 669 z_elem_c(:,:,1) = 0._wp 670 z_elem_b(:,:,1) = 1._wp 671 ! 672 ! One level below 673 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) ))) 674 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 675 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 676 z_elem_a(:,:,2) = 0._wp 677 z_elem_c(:,:,2) = 0._wp 678 z_elem_b(:,:,2) = 1._wp 679 ! 680 ! 681 ENDIF 555 682 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 556 ! 557 ! Surface value: Dirichlet 558 zdep(:,:) = zhsro(:,:) * rl_sf 559 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 560 z_elem_a(:,:,1) = psi(:,:,1) 561 z_elem_c(:,:,1) = 0._wp 562 z_elem_b(:,:,1) = 1._wp 563 ! 564 ! Neumann condition at k=2 565 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 566 z_elem_a(:,:,2) = 0._wp 567 ! 568 ! Set psi vertical flux at the surface: 569 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 570 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 571 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 572 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 573 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 574 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 575 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 576 577 ! 578 ! 683 IF( ln_phioc ) THEN ! Wave induced mixing case with forced/coupled fields 684 IF( nn_wmix==jp_janssen ) THEN 685 ! 686 zdep(:,:) = rl_sf * zhsro(:,:) 687 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 688 z_elem_a(:,:,1) = psi(:,:,1) 689 z_elem_c(:,:,1) = 0._wp 690 z_elem_b(:,:,1) = 1._wp 691 ! 692 ! Neumann condition at k=2 693 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 694 z_elem_a(:,:,2) = 0._wp 695 ! 696 ! Set psi vertical flux at the surface: 697 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 698 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 699 & * en(:,:,1)**rmm * zdep 700 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 701 ! 702 ELSE IF( nn_wmix==jp_craigbanner ) THEN 703 ! Surface value: Dirichlet 704 zdep(:,:) = zhsro(:,:) * rl_sf 705 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 706 z_elem_a(:,:,1) = psi(:,:,1) 707 z_elem_c(:,:,1) = 0._wp 708 z_elem_b(:,:,1) = 1._wp 709 ! 710 ! Neumann condition at k=2 711 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 712 z_elem_a(:,:,2) = 0._wp 713 ! 714 ! Set psi vertical flux at the surface: 715 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 716 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 717 zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * & 718 (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 719 zdep(:,:) = rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 720 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 721 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 722 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 723 ! 724 ENDIF 725 ELSE ! Wave induced mixing case with default values 726 ! Surface value: Dirichlet 727 zdep(:,:) = zhsro(:,:) * rl_sf 728 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 729 z_elem_a(:,:,1) = psi(:,:,1) 730 z_elem_c(:,:,1) = 0._wp 731 z_elem_b(:,:,1) = 1._wp 732 ! 733 ! Neumann condition at k=2 734 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 735 z_elem_a(:,:,2) = 0._wp 736 ! 737 ! Set psi vertical flux at the surface: 738 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 739 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 740 zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 741 (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 742 zdep(:,:) = rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 743 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 744 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 745 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 746 ! 747 ! 748 ENDIF 579 749 END SELECT 580 750 … … 867 1037 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 868 1038 & rn_clim_galp, ln_sigpsi, rn_hsro, & 869 & rn_crban , rn_charn, rn_frac_hs,&1039 & rn_crban_default, rn_charn, rn_frac_hs,& 870 1040 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 871 1041 & nn_stab_func, nn_clos … … 895 1065 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 896 1066 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 897 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban1067 WRITE(numout,*) ' Craig and Banner coefficient (default) rn_crban = ', rn_crban_default 898 1068 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 899 1069 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met … … 909 1079 910 1080 ! !* Check of some namelist values 911 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 912 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 913 IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 914 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 915 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 1081 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 1082 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 1083 IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 1084 IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' ) 1085 IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN 1086 CALL ctl_warn('W A R N I N G: ln_rough=.TRUE. but nn_z0_met is not 3 - resetting nn_z0_met to 3') 1087 nn_z0_met = 3 1088 ENDIF 1089 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 1090 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 916 1091 917 1092 SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure … … 1085 1260 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1086 1261 1087 IF ( rn_crban==0._wp ) THEN1262 IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN 1088 1263 rl_sf = vkarmn 1089 1264 ELSE … … 1124 1299 rc03 = rc02 * rc0 1125 1300 rc04 = rc03 * rc0 1126 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking1127 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking1128 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )1129 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer1301 rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 1302 rsbc_tke2 = rdt * rn_crban_default / rl_sf 1303 zcr = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1304 rtrans_default = 0.2_wp / zcr 1130 1305 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1131 1306 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness
Note: See TracChangeset
for help on using the changeset viewer.
