Changeset 5109
- Timestamp:
- 2015-03-02T15:05:44+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM
- Files:
-
- 3 edited
-
CONFIG/AMM12/EXP00/namelist_cfg (modified) (1 diff)
-
CONFIG/SHARED/namelist_ref (modified) (1 diff)
-
NEMO/OPA_SRC/ZDF/zdfgls.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
r4839 r5109 355 355 !----------------------------------------------------------------------- 356 356 rn_charn = 100000. ! Charnock constant for wb induced roughness length 357 nn_z0_met = 1 ! Method for surface roughness computation (0/1/2) 357 358 / 358 359 !----------------------------------------------------------------------- -
trunk/NEMOGCM/CONFIG/SHARED/namelist_ref
r5102 r5109 961 961 &namzdf_gls ! GLS vertical diffusion ("key_zdfgls") 962 962 !----------------------------------------------------------------------- 963 rn_emin = 1.e- 6! minimum value of e [m2/s2]963 rn_emin = 1.e-7 ! minimum value of e [m2/s2] 964 964 rn_epsmin = 1.e-12 ! minimum value of eps [m2/s3] 965 965 ln_length_lim = .true. ! limit on the dissipation rate under stable stratification (Galperin et al., 1988) 966 rn_clim_galp = 0.53 ! galperin limit 967 ln_crban = .true. ! Use Craig & Banner (1994) surface wave mixing parametrisation 966 rn_clim_galp = 0.267 ! galperin limit 968 967 ln_sigpsi = .true. ! Activate or not Burchard 2001 mods on psi schmidt number in the wb case 969 968 rn_crban = 100. ! Craig and Banner 1994 constant for wb tke flux 970 969 rn_charn = 70000. ! Charnock constant for wb induced roughness length 971 nn_tkebc_surf = 1 ! surface tke condition (0/1/2=Dir/Neum/Dir Mellor-Blumberg) 972 nn_tkebc_bot = 1 ! bottom tke condition (0/1=Dir/Neum) 973 nn_psibc_surf = 1 ! surface psi condition (0/1/2=Dir/Neum/Dir Mellor-Blumberg) 974 nn_psibc_bot = 1 ! bottom psi condition (0/1=Dir/Neum) 975 nn_stab_func = 2 ! stability function (0=Galp, 1= KC94, 2=CanutoA, 3=CanutoB) 976 nn_clos = 1 ! predefined closure type (0=MY82, 1=k-eps, 2=k-w, 3=Gen) 970 rn_hsro = 0.02 ! Minimum surface roughness 971 rn_frac_hs = 1.3 ! Fraction of wave height as roughness (if nn_z0_met=2) 972 nn_z0_met = 2 ! Method for surface roughness computation (0/1/2) 973 nn_bc_surf = 1 ! surface condition (0/1=Dir/Neum) 974 nn_bc_bot = 1 ! bottom condition (0/1=Dir/Neum) 975 nn_stab_func = 2 ! stability function (0=Galp, 1= KC94, 2=CanutoA, 3=CanutoB) 976 nn_clos = 1 ! predefined closure type (0=MY82, 1=k-eps, 2=k-w, 3=Gen) 977 977 / 978 978 !----------------------------------------------------------------------- -
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r4839 r5109 20 20 USE domvvl ! ocean space and time domain : variable volume layer 21 21 USE zdf_oce ! ocean vertical physics 22 USE zdfbfr ! bottom friction (only for rn_bfrz0) 22 23 USE sbc_oce ! surface boundary condition: ocean 23 24 USE phycst ! physical constants … … 52 53 53 54 ! !! ** Namelist namzdf_gls ** 54 LOGICAL :: ln_crban ! =T use Craig and Banner scheme55 55 LOGICAL :: ln_length_lim ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 56 56 LOGICAL :: ln_sigpsi ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 57 INTEGER :: nn_tkebc_surf ! TKE surface boundary condition (=0/1) 58 INTEGER :: nn_tkebc_bot ! TKE bottom boundary condition (=0/1) 59 INTEGER :: nn_psibc_surf ! PSI surface boundary condition (=0/1) 60 INTEGER :: nn_psibc_bot ! PSI bottom boundary condition (=0/1) 57 INTEGER :: nn_bc_surf ! surface boundary condition (=0/1) 58 INTEGER :: nn_bc_bot ! bottom boundary condition (=0/1) 59 INTEGER :: nn_z0_met ! Method for surface roughness computation 61 60 INTEGER :: nn_stab_func ! stability functions G88, KC or Canuto (=0/1/2) 62 61 INTEGER :: nn_clos ! closure 0/1/2/3 MY82/k-eps/k-w/gen … … 66 65 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 67 66 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 68 69 REAL(wp) :: hsro = 0.003_wp ! Minimum surface roughness70 REAL(wp) :: hbro = 0.003_wp ! Bottom roughness (m) 67 REAL(wp) :: rn_hsro ! Minimum surface roughness 68 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 69 71 70 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 72 71 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 … … 96 95 REAL(wp) :: rm7 = 0.0_wp 97 96 REAL(wp) :: rm8 = 0.318_wp 98 97 REAL(wp) :: rtrans = 0.1_wp 99 98 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 100 REAL(wp) :: rc03_sqrt2_galp ! - - - - 101 REAL(wp) :: rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke ! - - - - 102 REAL(wp) :: rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 103 REAL(wp) :: rsbc_mb , rsbc_std , rsbc_zs ! - - - - 99 REAL(wp) :: rsbc_tke1, rsbc_tke2, rfact_tke ! - - - - 100 REAL(wp) :: rsbc_psi1, rsbc_psi2, rfact_psi ! - - - - 101 REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - - 104 102 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - 105 103 REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 ! - - - - … … 147 145 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 148 146 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 147 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 149 148 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 150 149 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) … … 153 152 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 154 153 REAL(wp), POINTER, DIMENSION(:,:,:) :: eps ! dissipation rate 155 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T) 156 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a, z_elem_b, z_elem_c, psi 154 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 155 REAL(wp), POINTER, DIMENSION(:,:,:) :: psi ! psi at time now 156 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a ! element of the first matrix diagonal 157 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 158 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 157 159 !!-------------------------------------------------------------------- 158 160 ! 159 161 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 160 162 ! 161 CALL wrk_alloc( jpi,jpj, zdep, z flxs, zhsro )162 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )163 163 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 164 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 165 164 166 ! Preliminary computing 165 167 … … 174 176 175 177 ! Compute surface and bottom friction at T-points 176 !CDIR NOVERRCHK 177 DO jj = 2, jpjm1 178 !CDIR NOVERRCHK 179 DO ji = fs_2, fs_jpim1 ! vector opt. 180 ! 181 ! surface friction 178 !CDIR NOVERRCHK 179 DO jj = 2, jpjm1 180 !CDIR NOVERRCHK 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 ! 183 ! surface friction 182 184 ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 183 ! 184 ! bottom friction (explicit before friction) 185 ! Note that we chose here not to bound the friction as in dynbfr) 186 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 187 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 188 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 189 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 190 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 191 END DO 192 END DO 193 194 ! In case of breaking surface waves mixing, 195 ! Compute surface roughness length according to Charnock formula: 196 IF( ln_crban ) THEN ; zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 197 ELSE ; zhsro(:,:) = hsro 198 ENDIF 185 ! 186 ! bottom friction (explicit before friction) 187 ! Note that we chose here not to bound the friction as in dynbfr) 188 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 189 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 190 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 191 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 192 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 193 END DO 194 END DO 195 196 ! Set surface roughness length 197 SELECT CASE ( nn_z0_met ) 198 ! 199 CASE ( 0 ) ! Constant roughness 200 zhsro(:,:) = rn_hsro 201 CASE ( 1 ) ! Standard Charnock formula 202 zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 203 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 204 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 205 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 206 ! 207 END SELECT 199 208 200 209 ! Compute shear and dissipation rate … … 303 312 ! 304 313 ! Set surface condition on zwall_psi (1 at the bottom) 305 IF( ln_sigpsi ) THEN 306 zcoef = rsc_psi / rsc_psi0 307 DO jj = 2, jpjm1 308 DO ji = fs_2, fs_jpim1 ! vector opt. 309 zwall_psi(ji,jj,1) = zcoef 310 END DO 311 END DO 312 ENDIF 313 314 zwall_psi(:,:,1) = zwall_psi(:,:,2) 315 zwall_psi(:,:,jpk) = 1. 316 ! 314 317 ! Surface boundary condition on tke 315 318 ! --------------------------------- 316 319 ! 317 SELECT CASE ( nn_ tkebc_surf )320 SELECT CASE ( nn_bc_surf ) 318 321 ! 319 322 CASE ( 0 ) ! Dirichlet case 320 ! 321 IF (ln_crban) THEN ! Wave induced mixing case 322 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 323 ! ! balance between the production and the dissipation terms including the wave effect 324 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 325 z_elem_a(:,:,1) = en(:,:,1) 326 z_elem_c(:,:,1) = 0._wp 327 z_elem_b(:,:,1) = 1._wp 328 ! 329 ! one level below 330 en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 331 z_elem_a(:,:,2) = 0._wp 332 z_elem_c(:,:,2) = 0._wp 333 z_elem_b(:,:,2) = 1._wp 334 ! 335 ELSE ! No wave induced mixing case 336 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 337 ! ! balance between the production and the dissipation terms 338 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 339 z_elem_a(:,:,1) = en(:,:,1) 340 z_elem_c(:,:,1) = 0._wp 341 z_elem_b(:,:,1) = 1._wp 342 ! 343 ! one level below 344 en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 345 z_elem_a(:,:,2) = 0._wp 346 z_elem_c(:,:,2) = 0._wp 347 z_elem_b(:,:,2) = 1._wp 348 ! 349 ENDIF 350 ! 323 ! First level 324 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 325 en(:,:,1) = MAX(en(:,:,1), rn_emin) 326 z_elem_a(:,:,1) = en(:,:,1) 327 z_elem_c(:,:,1) = 0._wp 328 z_elem_b(:,:,1) = 1._wp 329 ! 330 ! One level below 331 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 332 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 333 z_elem_a(:,:,2) = 0._wp 334 z_elem_c(:,:,2) = 0._wp 335 z_elem_b(:,:,2) = 1._wp 336 ! 337 ! 351 338 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 352 ! 353 IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw 354 ! 355 ! Dirichlet conditions at k=1 (Cosmetic) 356 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 357 z_elem_a(:,:,1) = en(:,:,1) 358 z_elem_c(:,:,1) = 0._wp 359 z_elem_b(:,:,1) = 1._wp 360 ! at k=2, set de/dz=Fw 361 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 362 z_elem_a(:,:,2) = 0._wp 363 zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 364 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 365 ! 366 ELSE ! No wave induced mixing case: d(e)/dz=0. 367 ! 368 ! Dirichlet conditions at k=1 (Cosmetic) 369 en(:,:,1) = MAX( rc02r * 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=0.: 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 ! 377 ENDIF 378 ! 339 ! 340 ! Dirichlet conditions at k=1 341 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 342 en(:,:,1) = MAX(en(:,:,1), 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 ! at k=2, set de/dz=Fw 348 !cbr 349 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 350 z_elem_a(:,:,2) = 0._wp 351 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 352 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 353 354 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 355 ! 356 ! 379 357 END SELECT 380 358 … … 382 360 ! -------------------------------- 383 361 ! 384 SELECT CASE ( nn_ tkebc_bot )362 SELECT CASE ( nn_bc_bot ) 385 363 ! 386 364 CASE ( 0 ) ! Dirichlet … … 457 435 ! ! set the minimum value of tke 458 436 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 459 437 460 438 !!----------------------------------------!! 461 439 !! Solve prognostic equation for psi !! … … 560 538 ! --------------------------------- 561 539 ! 562 SELECT CASE ( nn_ psibc_surf )540 SELECT CASE ( nn_bc_surf ) 563 541 ! 564 542 CASE ( 0 ) ! Dirichlet boundary conditions 565 ! 566 IF( ln_crban ) THEN ! Wave induced mixing case 567 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 568 ! ! balance between the production and the dissipation terms including the wave effect 569 zdep(:,:) = rl_sf * zhsro(:,:) 570 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 571 z_elem_a(:,:,1) = psi(:,:,1) 572 z_elem_c(:,:,1) = 0._wp 573 z_elem_b(:,:,1) = 1._wp 574 ! 575 ! one level below 576 zex1 = (rmm*ra_sf+rnn) 577 zex2 = (rmm*ra_sf) 578 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 579 psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 580 z_elem_a(:,:,2) = 0._wp 581 z_elem_c(:,:,2) = 0._wp 582 z_elem_b(:,:,2) = 1._wp 583 ! 584 ELSE ! No wave induced mixing case 585 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 586 ! ! balance between the production and the dissipation terms 587 ! 588 zdep(:,:) = vkarmn * zhsro(:,:) 589 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 590 z_elem_a(:,:,1) = psi(:,:,1) 591 z_elem_c(:,:,1) = 0._wp 592 z_elem_b(:,:,1) = 1._wp 593 ! 594 ! one level below 595 zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 596 psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 597 z_elem_a(:,:,2) = 0._wp 598 z_elem_c(:,:,2) = 0._wp 599 z_elem_b(:,:,2) = 1. 600 ! 601 ENDIF 602 ! 543 ! 544 ! Surface value 545 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 546 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 547 z_elem_a(:,:,1) = psi(:,:,1) 548 z_elem_c(:,:,1) = 0._wp 549 z_elem_b(:,:,1) = 1._wp 550 ! 551 ! One level below 552 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 553 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 554 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 555 z_elem_a(:,:,2) = 0._wp 556 z_elem_c(:,:,2) = 0._wp 557 z_elem_b(:,:,2) = 1._wp 558 ! 559 ! 603 560 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 604 ! 605 IF( ln_crban ) THEN ! Wave induced mixing case 606 ! 607 zdep(:,:) = rl_sf * zhsro(:,:) 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 ! Neumann condition at k=2 614 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 615 z_elem_a(:,:,2) = 0._wp 616 ! 617 ! Set psi vertical flux at the surface: 618 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 619 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 620 & * en(:,:,1)**rmm * zdep 621 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 622 ! 623 ELSE ! No wave induced mixing 624 ! 625 zdep(:,:) = vkarmn * zhsro(:,:) 626 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 627 z_elem_a(:,:,1) = psi(:,:,1) 628 z_elem_c(:,:,1) = 0._wp 629 z_elem_b(:,:,1) = 1._wp 630 ! 631 ! Neumann condition at k=2 632 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 633 z_elem_a(ji,jj,2) = 0._wp 634 ! 635 ! Set psi vertical flux at the surface: 636 zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) 637 zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp) 638 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 639 ! 640 ENDIF 641 ! 561 ! 562 ! Surface value: Dirichlet 563 zdep(:,:) = zhsro(:,:) * rl_sf 564 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 565 z_elem_a(:,:,1) = psi(:,:,1) 566 z_elem_c(:,:,1) = 0._wp 567 z_elem_b(:,:,1) = 1._wp 568 ! 569 ! Neumann condition at k=2 570 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 571 z_elem_a(:,:,2) = 0._wp 572 ! 573 ! Set psi vertical flux at the surface: 574 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 575 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 576 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 577 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 578 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 579 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 580 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 581 582 ! 583 ! 642 584 END SELECT 643 585 … … 645 587 ! -------------------------------- 646 588 ! 647 SELECT CASE ( nn_psibc_bot ) 589 SELECT CASE ( nn_bc_bot ) 590 ! 648 591 ! 649 592 CASE ( 0 ) ! Dirichlet 650 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro593 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 651 594 ! ! Balance between the production and the dissipation terms 652 595 !CDIR NOVERRCHK … … 656 599 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 657 600 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 658 zdep(ji,jj) = vkarmn * hbro601 zdep(ji,jj) = vkarmn * rn_bfrz0 659 602 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 660 603 z_elem_a(ji,jj,ibot) = 0._wp … … 663 606 ! 664 607 ! Just above last level, Dirichlet condition again (GOTM like) 665 zdep(ji,jj) = vkarmn * ( hbro+ fse3t(ji,jj,ibotm1) )608 zdep(ji,jj) = vkarmn * ( rn_bfrz0 + fse3t(ji,jj,ibotm1) ) 666 609 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 667 610 z_elem_a(ji,jj,ibotm1) = 0._wp … … 681 624 ! 682 625 ! Bottom level Dirichlet condition: 683 zdep(ji,jj) = vkarmn * hbro626 zdep(ji,jj) = vkarmn * rn_bfrz0 684 627 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 685 628 ! … … 693 636 ! 694 637 ! Set psi vertical flux at the bottom: 695 zdep(ji,jj) = hbro+ 0.5_wp*fse3t(ji,jj,ibotm1)638 zdep(ji,jj) = rn_bfrz0 + 0.5_wp*fse3t(ji,jj,ibotm1) 696 639 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) & 697 640 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) … … 736 679 DO jj = 2, jpjm1 737 680 DO ji = fs_2, fs_jpim1 ! vector opt. 738 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk)681 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 739 682 END DO 740 683 END DO … … 783 726 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 784 727 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 785 mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk))728 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 786 729 END DO 787 730 END DO … … 847 790 ! Boundary conditions on stability functions for momentum (Neumann): 848 791 ! Lines below are useless if GOTM style Dirichlet conditions are used 849 zcoef = rcm_sf / SQRT( 2._wp ) 792 793 avmv(:,:,1) = avmv(:,:,2) 794 850 795 DO jj = 2, jpjm1 851 796 DO ji = fs_2, fs_jpim1 ! vector opt. 852 avmv(ji,jj,1) = zcoef 853 END DO 854 END DO 855 zcoef = rc0 / SQRT( 2._wp ) 856 DO jj = 2, jpjm1 857 DO ji = fs_2, fs_jpim1 ! vector opt. 858 avmv(ji,jj,mbkt(ji,jj)+1) = zcoef 797 avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj)) 859 798 END DO 860 799 END DO … … 900 839 avmv_k(:,:,:) = avmv(:,:,:) 901 840 ! 902 CALL wrk_dealloc( jpi,jpj, zdep, z flxs, zhsro )841 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 903 842 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 904 843 ! … … 932 871 !! 933 872 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 934 & rn_clim_galp, ln_crban, ln_sigpsi, & 935 & rn_crban, rn_charn, & 936 & nn_tkebc_surf, nn_tkebc_bot, & 937 & nn_psibc_surf, nn_psibc_bot, & 873 & rn_clim_galp, ln_sigpsi, rn_hsro, & 874 & rn_crban, rn_charn, rn_frac_hs, & 875 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 938 876 & nn_stab_func, nn_clos 939 877 !!---------------------------------------------------------- … … 955 893 WRITE(numout,*) '~~~~~~~~~~~~' 956 894 WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' 957 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 958 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 959 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 960 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 961 WRITE(numout,*) ' TKE Surface boundary condition nn_tkebc_surf = ', nn_tkebc_surf 962 WRITE(numout,*) ' TKE Bottom boundary condition nn_tkebc_bot = ', nn_tkebc_bot 963 WRITE(numout,*) ' PSI Surface boundary condition nn_psibc_surf = ', nn_psibc_surf 964 WRITE(numout,*) ' PSI Bottom boundary condition nn_psibc_bot = ', nn_psibc_bot 965 WRITE(numout,*) ' Craig and Banner scheme ln_crban = ', ln_crban 966 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 895 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 896 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 897 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 898 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 899 WRITE(numout,*) ' TKE Surface boundary condition nn_bc_surf = ', nn_bc_surf 900 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 901 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 967 902 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban 968 903 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 904 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met 905 WRITE(numout,*) ' Wave height frac. (used if nn_z0_met=2) rn_frac_hs = ', rn_frac_hs 969 906 WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func 970 907 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 971 WRITE(numout,*) ' Hard coded parameters' 972 WRITE(numout,*) ' Surface roughness (m) hsro = ', hsro 973 WRITE(numout,*) ' Bottom roughness (m) hbro = ', hbro 908 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 909 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 974 910 ENDIF 975 911 … … 978 914 979 915 ! !* Check of some namelist values 980 IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' ) 981 IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' ) 982 IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) 983 IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) 916 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 917 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 918 IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 984 919 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 985 920 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) … … 1001 936 SELECT CASE ( nn_stab_func ) 1002 937 CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions 1003 CASE( 2 ) ; rpsi3m = 2. 38_wp ! Canuto A stability functions938 CASE( 2 ) ; rpsi3m = 2.62_wp ! Canuto A stability functions 1004 939 CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified) 1005 940 END SELECT … … 1012 947 rnn = -1._wp 1013 948 rsc_tke = 1._wp 1014 rsc_psi = 1. 3_wp ! Schmidt number for psi949 rsc_psi = 1.2_wp ! Schmidt number for psi 1015 950 rpsi1 = 1.44_wp 1016 951 rpsi3p = 1._wp … … 1140 1075 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1141 1076 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 1142 IF( ln_sigpsi .AND. ln_crban ) THEN 1143 zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn 1144 rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf ) & 1145 & * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn & 1146 & + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm ) 1077 IF( ln_sigpsi ) THEN 1078 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1079 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1080 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work 1081 ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn 1082 rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.)) 1147 1083 ELSE 1148 1084 rsc_psi0 = rsc_psi … … 1151 1087 ! !* Shear free turbulence parameters 1152 1088 ! 1153 ra_sf = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke ) & 1154 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1155 rl_sf = rc0 * SQRT( rc0 / rcm_sf ) & 1156 & * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke & 1157 & + 12._wp * rsc_psi0 * rpsi2 & 1158 & - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) ) ) & 1159 & / ( 12._wp*rnn*rnn ) ) 1089 ra_sf = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) & 1090 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1091 1092 IF ( rn_crban==0._wp ) THEN 1093 rl_sf = vkarmn 1094 ELSE 1095 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke & 1096 & + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) & 1097 & *SQRT(rsc_tke*(rsc_tke & 1098 & + 24._wp*rsc_psi0*rpsi2)) ) & 1099 & /(12._wp*rnn**2.) & 1100 & ) 1101 ENDIF 1160 1102 1161 1103 ! … … 1187 1129 rc03 = rc02 * rc0 1188 1130 rc04 = rc03 * rc0 1189 rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp 1190 rsbc_mb = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp) ! Surf. bound. cond. from Mellor and Blumberg 1191 rsbc_std = 3.75_wp ! Surf. bound. cond. standard (prod=diss) 1192 rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53. Dirichlet + Wave breaking 1193 rsbc_tke2 = 0.5_wp / rau0 1194 rsbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking 1195 rsbc_zs = rn_charn / grav ! Charnock formula 1196 rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 1197 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1198 rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 1199 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1200 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1131 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1132 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1133 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1134 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1135 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1136 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1137 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1138 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1139 1140 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1141 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1201 1142 1202 1143 ! !* Wall proximity function … … 1257 1198 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1258 1199 en (:,:,:) = rn_emin 1259 mxln(:,:,:) = 0.0 011200 mxln(:,:,:) = 0.05 1260 1201 avt_k (:,:,:) = avt (:,:,:) 1261 1202 avm_k (:,:,:) = avm (:,:,:) … … 1267 1208 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1268 1209 en (:,:,:) = rn_emin 1269 mxln(:,:,:) = 0.0 011210 mxln(:,:,:) = 0.05 1270 1211 ENDIF 1271 1212 ! … … 1273 1214 ! ! ------------------- 1274 1215 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1275 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1216 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1276 1217 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1277 1218 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1278 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1219 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1279 1220 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1280 1221 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln )
Note: See TracChangeset
for help on using the changeset viewer.
