- Timestamp:
- 2017-03-17T15:21:06+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7481 r7809 24 24 USE phycst ! physical constants 25 25 USE zdfmxl ! mixed layer 26 USE sbcwave , ONLY: hsw ! significant wave height26 USE sbcwave, ONLY: hsw,rn_crban 27 27 ! 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 48 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 49 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 50 54 51 55 ! !! ** Namelist namzdf_gls ** … … 61 65 REAL(wp) :: rn_emin ! minimum value of TKE (m2/s2) 62 66 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 63 REAL(wp) :: rn_crban 67 REAL(wp) :: rn_crban_default ! Craig and Banner constant for surface breaking waves mixing 64 68 REAL(wp) :: rn_hsro ! Minimum surface roughness 65 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met >1)69 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met = 1) 66 70 67 71 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters … … 94 98 REAL(wp) :: rtrans = 0.1_wp 95 99 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 96 REAL(wp) :: rsbc_tke1 , rsbc_tke2, rfact_tke! - - - -97 REAL(wp) :: rsbc_psi 1, rsbc_psi2, rfact_psi ! - - - -100 REAL(wp) :: rsbc_tke1_default, rsbc_tke2, rfact_tke ! - - - - 101 REAL(wp) :: rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 98 102 REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - - 99 103 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - … … 117 121 !! *** FUNCTION zdf_gls_alloc *** 118 122 !!---------------------------------------------------------------------- 119 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 120 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 123 ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) , & 124 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , & 125 & rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) , & 126 & rsbc_psi1(jpi,jpj), STAT= zdf_gls_alloc ) 121 127 ! 122 128 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 163 169 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 164 170 171 172 ! variable initialization 173 IF( ln_phioc ) THEN 174 rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53.Dirichlet + Wave breaking 175 rsbc_tke3(:,:) = rdt * rn_crban(:,:) ! Neumann + Wave breaking 176 rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 177 ENDIF 178 165 179 IF( kt /= nit000 ) THEN ! restore before value to compute tke 166 180 avt (:,:,:) = avt_k (:,:,:) … … 200 214 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 201 215 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 202 zhsro(:,:) = hsw(:,:) 216 WHERE( hsw == 0._wp ) ! surface roughness length according to Charnock formula when sign. wave height 0 217 zhsro = MAX(rn_charn / grav * ustars2, rn_hsro) 218 ELSEWHERE 219 zhsro = MAX(hsw, rn_hsro) 220 END WHERE 203 221 END SELECT 204 222 … … 317 335 ! 318 336 CASE ( 0 ) ! Dirichlet case 319 ! First level 320 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 321 en(:,:,1) = MAX(en(:,:,1), rn_emin) 322 z_elem_a(:,:,1) = en(:,:,1) 323 z_elem_c(:,:,1) = 0._wp 324 z_elem_b(:,:,1) = 1._wp 325 ! 326 ! One level below 327 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 328 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 329 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 330 z_elem_a(:,:,2) = 0._wp 331 z_elem_c(:,:,2) = 0._wp 332 z_elem_b(:,:,2) = 1._wp 333 ! 334 ! 337 IF( ln_phioc ) THEN ! wave induced mixing case with forced/coupled fields 338 ! First level 339 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 340 z_elem_a(:,:,1) = en(:,:,1) 341 z_elem_c(:,:,1) = 0._wp 342 z_elem_b(:,:,1) = 1._wp 343 344 ! One level below 345 en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 346 z_elem_a(:,:,2) = 0._wp 347 z_elem_c(:,:,2) = 0._wp 348 z_elem_b(:,:,2) = 1._wp 349 ELSE ! wave induced mixing case with default values 350 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 351 en(:,:,1) = MAX(en(:,:,1), rn_emin) 352 z_elem_a(:,:,1) = en(:,:,1) 353 z_elem_c(:,:,1) = 0._wp 354 z_elem_b(:,:,1) = 1._wp 355 ! 356 ! One level below 357 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) & 358 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 359 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 360 z_elem_a(:,:,2) = 0._wp 361 z_elem_c(:,:,2) = 0._wp 362 z_elem_b(:,:,2) = 1._wp 363 ! 364 ! 365 ENDIF 335 366 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 336 ! 337 ! Dirichlet conditions at k=1 338 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 339 en(:,:,1) = MAX(en(:,:,1), rn_emin) 340 z_elem_a(:,:,1) = en(:,:,1) 341 z_elem_c(:,:,1) = 0._wp 342 z_elem_b(:,:,1) = 1._wp 343 ! 344 ! at k=2, set de/dz=Fw 345 !cbr 346 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 347 z_elem_a(:,:,2) = 0._wp 348 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 349 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 350 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 351 352 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 353 ! 354 ! 367 IF( ln_phioc ) THEN ! Shear free case: d(e)/dz=Fw with forced/coupled fields 368 ! Dirichlet conditions at k=1 369 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 370 z_elem_a(:,:,1) = en(:,:,1) 371 z_elem_c(:,:,1) = 0._wp 372 z_elem_b(:,:,1) = 1._wp 373 ! at k=2, set de/dz=Fw 374 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 375 z_elem_a(:,:,2) = 0._wp 376 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 377 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 378 ELSE ! Shear free case: d(e)/dz=Fw with default values 379 ! Dirichlet conditions at k=1 380 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 381 en(:,:,1) = MAX(en(:,:,1), rn_emin) 382 z_elem_a(:,:,1) = en(:,:,1) 383 z_elem_c(:,:,1) = 0._wp 384 z_elem_b(:,:,1) = 1._wp 385 ! 386 ! at k=2, set de/dz=Fw 387 !cbr 388 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 389 z_elem_a(:,:,2) = 0._wp 390 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 391 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 392 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 393 394 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 395 ! 396 ! 397 ENDIF 355 398 END SELECT 356 399 … … 539 582 ! 540 583 CASE ( 0 ) ! Dirichlet boundary conditions 541 ! 542 ! Surface value 543 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 544 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 545 z_elem_a(:,:,1) = psi(:,:,1) 546 z_elem_c(:,:,1) = 0._wp 547 z_elem_b(:,:,1) = 1._wp 548 ! 549 ! One level below 550 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 551 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 552 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 553 z_elem_a(:,:,2) = 0._wp 554 z_elem_c(:,:,2) = 0._wp 555 z_elem_b(:,:,2) = 1._wp 556 ! 557 ! 584 IF( ln_phioc ) THEN ! Wave induced mixing case 585 ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 586 ! balance between the production and the 587 ! dissipation terms including the wave effect 588 ! Surface value 589 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 590 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 591 z_elem_a(:,:,1) = psi(:,:,1) 592 z_elem_c(:,:,1) = 0._wp 593 z_elem_b(:,:,1) = 1._wp 594 ! 595 ! One level below 596 zex1 = (rmm*ra_sf+rnn) 597 zex2 = (rmm*ra_sf) 598 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 599 psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 600 z_elem_a(:,:,2) = 0._wp 601 z_elem_c(:,:,2) = 0._wp 602 z_elem_b(:,:,2) = 1._wp 603 ! 604 ! 605 ELSE ! Wave induced mixing case with default values 606 ! Surface value 607 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 608 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 609 z_elem_a(:,:,1) = psi(:,:,1) 610 z_elem_c(:,:,1) = 0._wp 611 z_elem_b(:,:,1) = 1._wp 612 ! 613 ! One level below 614 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 615 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 616 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 617 z_elem_a(:,:,2) = 0._wp 618 z_elem_c(:,:,2) = 0._wp 619 z_elem_b(:,:,2) = 1._wp 620 ! 621 ! 622 ENDIF 558 623 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 559 ! 560 ! Surface value: Dirichlet 561 zdep(:,:) = zhsro(:,:) * rl_sf 562 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 563 z_elem_a(:,:,1) = psi(:,:,1) 564 z_elem_c(:,:,1) = 0._wp 565 z_elem_b(:,:,1) = 1._wp 566 ! 567 ! Neumann condition at k=2 568 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 569 z_elem_a(:,:,2) = 0._wp 570 ! 571 ! Set psi vertical flux at the surface: 572 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 573 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 574 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 575 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 576 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 577 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 578 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 579 580 ! 581 ! 624 IF( ln_phioc ) THEN ! Wave induced mixing case with forced/coupled fields 625 ! 626 zdep(:,:) = rl_sf * zhsro(:,:) 627 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 628 z_elem_a(:,:,1) = psi(:,:,1) 629 z_elem_c(:,:,1) = 0._wp 630 z_elem_b(:,:,1) = 1._wp 631 ! 632 ! Neumann condition at k=2 633 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 634 z_elem_a(:,:,2) = 0._wp 635 ! 636 ! Set psi vertical flux at the surface: 637 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 638 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 639 & * en(:,:,1)**rmm * zdep 640 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 641 ! 642 ELSE ! Wave induced mixing case with default values 643 ! Surface value: Dirichlet 644 zdep(:,:) = zhsro(:,:) * rl_sf 645 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 646 z_elem_a(:,:,1) = psi(:,:,1) 647 z_elem_c(:,:,1) = 0._wp 648 z_elem_b(:,:,1) = 1._wp 649 ! 650 ! Neumann condition at k=2 651 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 652 z_elem_a(:,:,2) = 0._wp 653 ! 654 ! Set psi vertical flux at the surface: 655 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 656 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 657 zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 658 (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 659 zdep(:,:) = rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 660 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 661 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 662 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 663 ! 664 ! 665 ENDIF 582 666 END SELECT 583 667 … … 870 954 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 871 955 & rn_clim_galp, ln_sigpsi, rn_hsro, & 872 & rn_crban , rn_charn, rn_frac_hs,&956 & rn_crban_default, rn_charn, rn_frac_hs,& 873 957 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 874 958 & nn_stab_func, nn_clos … … 898 982 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 899 983 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 900 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban984 WRITE(numout,*) ' Craig and Banner coefficient (default) rn_crban = ', rn_crban_default 901 985 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 902 986 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met … … 915 999 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' ) 916 1000 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' ) 917 IF( nn_z0_met == 3 .AND. .NOT.ln_ wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )1001 IF( nn_z0_met == 3 .AND. .NOT.ln_sdw ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_sdw=T' ) 918 1002 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' ) 919 1003 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) … … 1089 1173 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1090 1174 1091 IF ( rn_crban==0._wp ) THEN1175 IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN 1092 1176 rl_sf = vkarmn 1093 1177 ELSE … … 1128 1212 rc03 = rc02 * rc0 1129 1213 rc04 = rc03 * rc0 1130 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking1131 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking1132 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )1133 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer1214 rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 1215 rsbc_tke2 = rdt * rn_crban_default / rl_sf 1216 zcr = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1217 rtrans = 0.2_wp / zcr 1134 1218 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1135 1219 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.