- Timestamp:
- 2015-07-10T13:28:53+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r5081 r5581 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 … … 47 48 48 49 ! !! ** Namelist namzdf_gls ** 49 LOGICAL :: ln_crban ! =T use Craig and Banner scheme50 50 LOGICAL :: ln_length_lim ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 51 51 LOGICAL :: ln_sigpsi ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 52 INTEGER :: nn_tkebc_surf ! TKE surface boundary condition (=0/1) 53 INTEGER :: nn_tkebc_bot ! TKE bottom boundary condition (=0/1) 54 INTEGER :: nn_psibc_surf ! PSI surface boundary condition (=0/1) 55 INTEGER :: nn_psibc_bot ! PSI bottom boundary condition (=0/1) 52 INTEGER :: nn_bc_surf ! surface boundary condition (=0/1) 53 INTEGER :: nn_bc_bot ! bottom boundary condition (=0/1) 54 INTEGER :: nn_z0_met ! Method for surface roughness computation 56 55 INTEGER :: nn_stab_func ! stability functions G88, KC or Canuto (=0/1/2) 57 56 INTEGER :: nn_clos ! closure 0/1/2/3 MY82/k-eps/k-w/gen … … 61 60 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 62 61 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 63 64 REAL(wp) :: hsro = 0.003_wp ! Minimum surface roughness65 REAL(wp) :: hbro = 0.003_wp ! Bottom roughness (m) 62 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) 64 66 65 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 67 66 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 … … 91 90 REAL(wp) :: rm7 = 0.0_wp 92 91 REAL(wp) :: rm8 = 0.318_wp 93 92 REAL(wp) :: rtrans = 0.1_wp 94 93 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 95 REAL(wp) :: rc03_sqrt2_galp ! - - - - 96 REAL(wp) :: rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke ! - - - - 97 REAL(wp) :: rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 98 REAL(wp) :: rsbc_mb , rsbc_std , rsbc_zs ! - - - - 94 REAL(wp) :: rsbc_tke1, rsbc_tke2, rfact_tke ! - - - - 95 REAL(wp) :: rsbc_psi1, rsbc_psi2, rfact_psi ! - - - - 96 REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - - 99 97 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - 100 98 REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 ! - - - - … … 140 138 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 141 139 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 142 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 143 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) … … 146 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 147 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: eps ! dissipation rate 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T) 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a, z_elem_b, z_elem_c, psi 147 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: psi ! psi at time now 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a ! element of the first matrix diagonal 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 151 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 150 152 !!-------------------------------------------------------------------- 151 153 ! 152 154 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 153 155 ! 154 CALL wrk_alloc( jpi,jpj, zdep, z flxs, zhsro )155 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )156 156 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 157 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 158 157 159 ! Preliminary computing 158 160 … … 167 169 168 170 ! Compute surface and bottom friction at T-points 169 !CDIR NOVERRCHK 170 DO jj = 2, jpjm1 171 !CDIR NOVERRCHK 172 DO ji = fs_2, fs_jpim1 ! vector opt. 173 ! 174 ! surface friction 171 !CDIR NOVERRCHK 172 DO jj = 2, jpjm1 173 !CDIR NOVERRCHK 174 DO ji = fs_2, fs_jpim1 ! vector opt. 175 ! 176 ! surface friction 175 177 ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 176 ! 177 ! bottom friction (explicit before friction) 178 ! Note that we chose here not to bound the friction as in dynbfr) 179 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 180 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 181 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 182 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 183 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 184 END DO 185 END DO 186 187 ! In case of breaking surface waves mixing, 188 ! Compute surface roughness length according to Charnock formula: 189 IF( ln_crban ) THEN ; zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 190 ELSE ; zhsro(:,:) = hsro 191 ENDIF 178 ! 179 ! bottom friction (explicit before friction) 180 ! Note that we chose here not to bound the friction as in dynbfr) 181 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 182 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 183 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 184 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 185 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 186 END DO 187 END DO 188 189 ! Set surface roughness length 190 SELECT CASE ( nn_z0_met ) 191 ! 192 CASE ( 0 ) ! Constant roughness 193 zhsro(:,:) = rn_hsro 194 CASE ( 1 ) ! Standard Charnock formula 195 zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 196 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 197 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 198 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 199 ! 200 END SELECT 192 201 193 202 ! Compute shear and dissipation rate … … 296 305 ! 297 306 ! Set surface condition on zwall_psi (1 at the bottom) 298 IF( ln_sigpsi ) THEN 299 zcoef = rsc_psi / rsc_psi0 300 DO jj = 2, jpjm1 301 DO ji = fs_2, fs_jpim1 ! vector opt. 302 zwall_psi(ji,jj,1) = zcoef 303 END DO 304 END DO 305 ENDIF 306 307 zwall_psi(:,:,1) = zwall_psi(:,:,2) 308 zwall_psi(:,:,jpk) = 1. 309 ! 307 310 ! Surface boundary condition on tke 308 311 ! --------------------------------- 309 312 ! 310 SELECT CASE ( nn_ tkebc_surf )313 SELECT CASE ( nn_bc_surf ) 311 314 ! 312 315 CASE ( 0 ) ! Dirichlet case 313 ! 314 IF (ln_crban) THEN ! Wave induced mixing case 315 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 316 ! ! balance between the production and the dissipation terms including the wave effect 317 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 318 z_elem_a(:,:,1) = en(:,:,1) 319 z_elem_c(:,:,1) = 0._wp 320 z_elem_b(:,:,1) = 1._wp 321 ! 322 ! one level below 323 en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 324 z_elem_a(:,:,2) = 0._wp 325 z_elem_c(:,:,2) = 0._wp 326 z_elem_b(:,:,2) = 1._wp 327 ! 328 ELSE ! No wave induced mixing case 329 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 330 ! ! balance between the production and the dissipation terms 331 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 332 z_elem_a(:,:,1) = en(:,:,1) 333 z_elem_c(:,:,1) = 0._wp 334 z_elem_b(:,:,1) = 1._wp 335 ! 336 ! one level below 337 en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 338 z_elem_a(:,:,2) = 0._wp 339 z_elem_c(:,:,2) = 0._wp 340 z_elem_b(:,:,2) = 1._wp 341 ! 342 ENDIF 343 ! 316 ! First level 317 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 318 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 z_elem_a(:,:,1) = en(:,:,1) 320 z_elem_c(:,:,1) = 0._wp 321 z_elem_b(:,:,1) = 1._wp 322 ! 323 ! One level below 324 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 325 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 326 z_elem_a(:,:,2) = 0._wp 327 z_elem_c(:,:,2) = 0._wp 328 z_elem_b(:,:,2) = 1._wp 329 ! 330 ! 344 331 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 345 ! 346 IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw 347 ! 348 ! Dirichlet conditions at k=1 (Cosmetic) 349 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 350 z_elem_a(:,:,1) = en(:,:,1) 351 z_elem_c(:,:,1) = 0._wp 352 z_elem_b(:,:,1) = 1._wp 353 ! at k=2, set de/dz=Fw 354 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 355 z_elem_a(:,:,2) = 0._wp 356 zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 357 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 358 ! 359 ELSE ! No wave induced mixing case: d(e)/dz=0. 360 ! 361 ! Dirichlet conditions at k=1 (Cosmetic) 362 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 363 z_elem_a(:,:,1) = en(:,:,1) 364 z_elem_c(:,:,1) = 0._wp 365 z_elem_b(:,:,1) = 1._wp 366 ! at k=2 set de/dz=0.: 367 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 368 z_elem_a(:,:,2) = 0._wp 369 ! 370 ENDIF 371 ! 332 ! 333 ! Dirichlet conditions at k=1 334 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 335 en(:,:,1) = MAX(en(:,:,1), rn_emin) 336 z_elem_a(:,:,1) = en(:,:,1) 337 z_elem_c(:,:,1) = 0._wp 338 z_elem_b(:,:,1) = 1._wp 339 ! 340 ! at k=2, set de/dz=Fw 341 !cbr 342 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 343 z_elem_a(:,:,2) = 0._wp 344 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 345 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 346 347 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 348 ! 349 ! 372 350 END SELECT 373 351 … … 375 353 ! -------------------------------- 376 354 ! 377 SELECT CASE ( nn_ tkebc_bot )355 SELECT CASE ( nn_bc_bot ) 378 356 ! 379 357 CASE ( 0 ) ! Dirichlet … … 450 428 ! ! set the minimum value of tke 451 429 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 452 430 453 431 !!----------------------------------------!! 454 432 !! Solve prognostic equation for psi !! … … 553 531 ! --------------------------------- 554 532 ! 555 SELECT CASE ( nn_ psibc_surf )533 SELECT CASE ( nn_bc_surf ) 556 534 ! 557 535 CASE ( 0 ) ! Dirichlet boundary conditions 558 ! 559 IF( ln_crban ) THEN ! Wave induced mixing case 560 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 561 ! ! balance between the production and the dissipation terms including the wave effect 562 zdep(:,:) = rl_sf * zhsro(:,:) 563 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 564 z_elem_a(:,:,1) = psi(:,:,1) 565 z_elem_c(:,:,1) = 0._wp 566 z_elem_b(:,:,1) = 1._wp 567 ! 568 ! one level below 569 zex1 = (rmm*ra_sf+rnn) 570 zex2 = (rmm*ra_sf) 571 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 572 psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 573 z_elem_a(:,:,2) = 0._wp 574 z_elem_c(:,:,2) = 0._wp 575 z_elem_b(:,:,2) = 1._wp 576 ! 577 ELSE ! No wave induced mixing case 578 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 579 ! ! balance between the production and the dissipation terms 580 ! 581 zdep(:,:) = vkarmn * zhsro(:,:) 582 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 583 z_elem_a(:,:,1) = psi(:,:,1) 584 z_elem_c(:,:,1) = 0._wp 585 z_elem_b(:,:,1) = 1._wp 586 ! 587 ! one level below 588 zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 589 psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 590 z_elem_a(:,:,2) = 0._wp 591 z_elem_c(:,:,2) = 0._wp 592 z_elem_b(:,:,2) = 1. 593 ! 594 ENDIF 595 ! 536 ! 537 ! Surface value 538 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 539 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 540 z_elem_a(:,:,1) = psi(:,:,1) 541 z_elem_c(:,:,1) = 0._wp 542 z_elem_b(:,:,1) = 1._wp 543 ! 544 ! One level below 545 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 546 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 547 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 548 z_elem_a(:,:,2) = 0._wp 549 z_elem_c(:,:,2) = 0._wp 550 z_elem_b(:,:,2) = 1._wp 551 ! 552 ! 596 553 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 597 ! 598 IF( ln_crban ) THEN ! Wave induced mixing case 599 ! 600 zdep(:,:) = rl_sf * zhsro(:,:) 601 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 602 z_elem_a(:,:,1) = psi(:,:,1) 603 z_elem_c(:,:,1) = 0._wp 604 z_elem_b(:,:,1) = 1._wp 605 ! 606 ! Neumann condition at k=2 607 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 608 z_elem_a(:,:,2) = 0._wp 609 ! 610 ! Set psi vertical flux at the surface: 611 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 612 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 613 & * en(:,:,1)**rmm * zdep 614 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 615 ! 616 ELSE ! No wave induced mixing 617 ! 618 zdep(:,:) = vkarmn * zhsro(:,:) 619 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 620 z_elem_a(:,:,1) = psi(:,:,1) 621 z_elem_c(:,:,1) = 0._wp 622 z_elem_b(:,:,1) = 1._wp 623 ! 624 ! Neumann condition at k=2 625 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 626 z_elem_a(ji,jj,2) = 0._wp 627 ! 628 ! Set psi vertical flux at the surface: 629 zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) 630 zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp) 631 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 632 ! 633 ENDIF 634 ! 554 ! 555 ! Surface value: Dirichlet 556 zdep(:,:) = zhsro(:,:) * rl_sf 557 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 558 z_elem_a(:,:,1) = psi(:,:,1) 559 z_elem_c(:,:,1) = 0._wp 560 z_elem_b(:,:,1) = 1._wp 561 ! 562 ! Neumann condition at k=2 563 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 564 z_elem_a(:,:,2) = 0._wp 565 ! 566 ! Set psi vertical flux at the surface: 567 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 568 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 569 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 570 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 571 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 572 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 573 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 574 575 ! 576 ! 635 577 END SELECT 636 578 … … 638 580 ! -------------------------------- 639 581 ! 640 SELECT CASE ( nn_psibc_bot ) 582 SELECT CASE ( nn_bc_bot ) 583 ! 641 584 ! 642 585 CASE ( 0 ) ! Dirichlet 643 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro586 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 644 587 ! ! Balance between the production and the dissipation terms 645 588 !CDIR NOVERRCHK … … 649 592 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 650 593 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 651 zdep(ji,jj) = vkarmn * hbro594 zdep(ji,jj) = vkarmn * rn_bfrz0 652 595 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 653 596 z_elem_a(ji,jj,ibot) = 0._wp … … 656 599 ! 657 600 ! Just above last level, Dirichlet condition again (GOTM like) 658 zdep(ji,jj) = vkarmn * ( hbro+ fse3t(ji,jj,ibotm1) )601 zdep(ji,jj) = vkarmn * ( rn_bfrz0 + fse3t(ji,jj,ibotm1) ) 659 602 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 660 603 z_elem_a(ji,jj,ibotm1) = 0._wp … … 674 617 ! 675 618 ! Bottom level Dirichlet condition: 676 zdep(ji,jj) = vkarmn * hbro619 zdep(ji,jj) = vkarmn * rn_bfrz0 677 620 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 678 621 ! … … 686 629 ! 687 630 ! Set psi vertical flux at the bottom: 688 zdep(ji,jj) = hbro+ 0.5_wp*fse3t(ji,jj,ibotm1)631 zdep(ji,jj) = rn_bfrz0 + 0.5_wp*fse3t(ji,jj,ibotm1) 689 632 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) & 690 633 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) … … 729 672 DO jj = 2, jpjm1 730 673 DO ji = fs_2, fs_jpim1 ! vector opt. 731 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk)674 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) 732 675 END DO 733 676 END DO … … 776 719 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 777 720 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 778 mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk))721 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 779 722 END DO 780 723 END DO … … 840 783 ! Boundary conditions on stability functions for momentum (Neumann): 841 784 ! Lines below are useless if GOTM style Dirichlet conditions are used 842 zcoef = rcm_sf / SQRT( 2._wp ) 785 786 avmv(:,:,1) = avmv(:,:,2) 787 843 788 DO jj = 2, jpjm1 844 789 DO ji = fs_2, fs_jpim1 ! vector opt. 845 avmv(ji,jj,1) = zcoef 846 END DO 847 END DO 848 zcoef = rc0 / SQRT( 2._wp ) 849 DO jj = 2, jpjm1 850 DO ji = fs_2, fs_jpim1 ! vector opt. 851 avmv(ji,jj,mbkt(ji,jj)+1) = zcoef 790 avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj)) 852 791 END DO 853 792 END DO … … 893 832 avmv_k(:,:,:) = avmv(:,:,:) 894 833 ! 895 CALL wrk_dealloc( jpi,jpj, zdep, z flxs, zhsro )834 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 896 835 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 897 836 ! … … 925 864 !! 926 865 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 927 & rn_clim_galp, ln_crban, ln_sigpsi, & 928 & rn_crban, rn_charn, & 929 & nn_tkebc_surf, nn_tkebc_bot, & 930 & nn_psibc_surf, nn_psibc_bot, & 866 & rn_clim_galp, ln_sigpsi, rn_hsro, & 867 & rn_crban, rn_charn, rn_frac_hs, & 868 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 931 869 & nn_stab_func, nn_clos 932 870 !!---------------------------------------------------------- … … 948 886 WRITE(numout,*) '~~~~~~~~~~~~' 949 887 WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' 950 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 951 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 952 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 953 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 954 WRITE(numout,*) ' TKE Surface boundary condition nn_tkebc_surf = ', nn_tkebc_surf 955 WRITE(numout,*) ' TKE Bottom boundary condition nn_tkebc_bot = ', nn_tkebc_bot 956 WRITE(numout,*) ' PSI Surface boundary condition nn_psibc_surf = ', nn_psibc_surf 957 WRITE(numout,*) ' PSI Bottom boundary condition nn_psibc_bot = ', nn_psibc_bot 958 WRITE(numout,*) ' Craig and Banner scheme ln_crban = ', ln_crban 959 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 888 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 889 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 890 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 891 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 892 WRITE(numout,*) ' TKE Surface boundary condition nn_bc_surf = ', nn_bc_surf 893 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 894 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 960 895 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban 961 896 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 897 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met 898 WRITE(numout,*) ' Wave height frac. (used if nn_z0_met=2) rn_frac_hs = ', rn_frac_hs 962 899 WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func 963 900 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 964 WRITE(numout,*) ' Hard coded parameters' 965 WRITE(numout,*) ' Surface roughness (m) hsro = ', hsro 966 WRITE(numout,*) ' Bottom roughness (m) hbro = ', hbro 901 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 902 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 967 903 ENDIF 968 904 … … 971 907 972 908 ! !* Check of some namelist values 973 IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' ) 974 IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' ) 975 IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) 976 IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) 909 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 910 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 911 IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 977 912 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 978 913 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) … … 994 929 SELECT CASE ( nn_stab_func ) 995 930 CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions 996 CASE( 2 ) ; rpsi3m = 2. 38_wp ! Canuto A stability functions931 CASE( 2 ) ; rpsi3m = 2.62_wp ! Canuto A stability functions 997 932 CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified) 998 933 END SELECT … … 1005 940 rnn = -1._wp 1006 941 rsc_tke = 1._wp 1007 rsc_psi = 1. 3_wp ! Schmidt number for psi942 rsc_psi = 1.2_wp ! Schmidt number for psi 1008 943 rpsi1 = 1.44_wp 1009 944 rpsi3p = 1._wp … … 1133 1068 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1134 1069 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 1135 IF( ln_sigpsi .AND. ln_crban ) THEN 1136 zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn 1137 rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf ) & 1138 & * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn & 1139 & + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm ) 1070 IF( ln_sigpsi ) THEN 1071 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1072 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1073 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work 1074 ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn 1075 rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.)) 1140 1076 ELSE 1141 1077 rsc_psi0 = rsc_psi … … 1144 1080 ! !* Shear free turbulence parameters 1145 1081 ! 1146 ra_sf = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke ) & 1147 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1148 rl_sf = rc0 * SQRT( rc0 / rcm_sf ) & 1149 & * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke & 1150 & + 12._wp * rsc_psi0 * rpsi2 & 1151 & - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) ) ) & 1152 & / ( 12._wp*rnn*rnn ) ) 1082 ra_sf = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) & 1083 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1084 1085 IF ( rn_crban==0._wp ) THEN 1086 rl_sf = vkarmn 1087 ELSE 1088 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke & 1089 & + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) & 1090 & *SQRT(rsc_tke*(rsc_tke & 1091 & + 24._wp*rsc_psi0*rpsi2)) ) & 1092 & /(12._wp*rnn**2.) & 1093 & ) 1094 ENDIF 1153 1095 1154 1096 ! … … 1180 1122 rc03 = rc02 * rc0 1181 1123 rc04 = rc03 * rc0 1182 rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp 1183 rsbc_mb = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp) ! Surf. bound. cond. from Mellor and Blumberg 1184 rsbc_std = 3.75_wp ! Surf. bound. cond. standard (prod=diss) 1185 rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53. Dirichlet + Wave breaking 1186 rsbc_tke2 = 0.5_wp / rau0 1187 rsbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking 1188 rsbc_zs = rn_charn / grav ! Charnock formula 1189 rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 1190 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1191 rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 1192 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1193 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1124 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1125 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1126 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1127 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1128 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1129 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1130 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1131 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1132 1133 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1134 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1194 1135 1195 1136 ! !* Wall proximity function … … 1250 1191 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1251 1192 en (:,:,:) = rn_emin 1252 mxln(:,:,:) = 0.001 1193 mxln(:,:,:) = 0.05 1194 avt_k (:,:,:) = avt (:,:,:) 1195 avm_k (:,:,:) = avm (:,:,:) 1196 avmu_k(:,:,:) = avmu(:,:,:) 1197 avmv_k(:,:,:) = avmv(:,:,:) 1253 1198 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1254 1199 ENDIF … … 1256 1201 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1257 1202 en (:,:,:) = rn_emin 1258 mxln(:,:,:) = 0.0 011203 mxln(:,:,:) = 0.05 1259 1204 ENDIF 1260 1205 ! … … 1262 1207 ! ! ------------------- 1263 1208 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1264 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1209 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1265 1210 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1266 1211 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1267 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1212 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1268 1213 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1269 1214 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln )
Note: See TracChangeset
for help on using the changeset viewer.
