Changeset 2299
- Timestamp:
- 2010-10-20T12:14:24+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r2293 r2299 59 59 REAL(wp) :: rn_charn = 2.e+5_wp ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 60 60 REAL(wp) :: rn_crban = 100._wp ! Craig and Banner constant for surface breaking waves mixing 61 REAL(wp) :: rn_cm_sf = 0.73_wp ! Shear free turbulence parameters 62 REAL(wp) :: r n_a_sf = -2.0_wp ! Must be negative -2 < rn_a_sf < -163 REAL(wp) :: r n_l_sf = 0.2_wp ! 0 <rn_l_sf<vkarmn64 REAL(wp) :: r n_ghmin = -0.2865 REAL(wp) :: r n_gh0 = 0.032966 REAL(wp) :: r n_ghcri = 0.0367 REAL(wp) :: r n_a1 = 0.92_wp68 REAL(wp) :: r n_a2 = 0.74_wp69 REAL(wp) :: r n_b1 = 16.60_wp70 REAL(wp) :: r n_b2 = 10.10_wp71 REAL(wp) :: r n_e1 = 1.80_wp72 REAL(wp) :: r n_e2= 1.33_wp73 REAL(wp) :: r n_l1= 0.107_wp74 REAL(wp) :: r n_l2= 0.0032_wp75 REAL(wp) :: r n_l3= 0.0864_wp76 REAL(wp) :: r n_l4= 0.12_wp77 REAL(wp) :: r n_l5= 11.9_wp78 REAL(wp) :: r n_l6= 0.4_wp79 REAL(wp) :: r n_l7= 0.0_wp80 REAL(wp) :: r n_l8= 0.48_wp81 REAL(wp) :: r n_m1= 0.127_wp82 REAL(wp) :: r n_m2= 0.00336_wp83 REAL(wp) :: r n_m3= 0.0906_wp84 REAL(wp) :: r n_m4= 0.101_wp85 REAL(wp) :: r n_m5= 11.2_wp86 REAL(wp) :: r n_m6= 0.4_wp87 REAL(wp) :: r n_m7= 0.0_wp88 REAL(wp) :: r n_m8= 0.318_wp89 REAL(wp) :: r n_c0290 REAL(wp) :: r n_c02r91 REAL(wp) :: r n_c0392 REAL(wp) :: r n_c0493 REAL(wp) :: r n_c03_sqrt2_galp94 REAL(wp) :: r n_sbc_mb95 REAL(wp) :: r n_sbc_std96 REAL(wp) :: r n_sbc_tke197 REAL(wp) :: r n_sbc_tke298 REAL(wp) :: r n_sbc_tke399 REAL(wp) :: r n_sbc_psi1100 REAL(wp) :: r n_sbc_psi2101 REAL(wp) :: r n_sbc_psi3102 REAL(wp) :: r n_sbc_zs103 REAL(wp) :: zfact_tke, zfact_psi104 REAL(wp) :: r n_c0, rn_c2, rn_c3, rn_f6, rn_cff, rn_c_diff105 REAL(wp) :: r n_s0, rn_s1, rn_s2, rn_s4, rn_s5, rn_s6106 REAL(wp) :: r n_d0, rn_d1, rn_d2, rn_d3, rn_d4, rn_d5107 REAL(wp) :: rn _sc_tke, rn_sc_psi, rn_psi1, rn_psi2, rn_psi3, rn_sc_psi0108 REAL(wp) :: r n_psi3m, rn_psi3p, zpp, zmm, znn, epslim61 62 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 63 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 64 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 65 REAL(wp) :: rghmin = -0.28 66 REAL(wp) :: rgh0 = 0.0329 67 REAL(wp) :: rghcri = 0.03 68 REAL(wp) :: ra1 = 0.92_wp 69 REAL(wp) :: ra2 = 0.74_wp 70 REAL(wp) :: rb1 = 16.60_wp 71 REAL(wp) :: rb2 = 10.10_wp 72 REAL(wp) :: re2 = 1.33_wp 73 REAL(wp) :: rl1 = 0.107_wp 74 REAL(wp) :: rl2 = 0.0032_wp 75 REAL(wp) :: rl3 = 0.0864_wp 76 REAL(wp) :: rl4 = 0.12_wp 77 REAL(wp) :: rl5 = 11.9_wp 78 REAL(wp) :: rl6 = 0.4_wp 79 REAL(wp) :: rl7 = 0.0_wp 80 REAL(wp) :: rl8 = 0.48_wp 81 REAL(wp) :: rm1 = 0.127_wp 82 REAL(wp) :: rm2 = 0.00336_wp 83 REAL(wp) :: rm3 = 0.0906_wp 84 REAL(wp) :: rm4 = 0.101_wp 85 REAL(wp) :: rm5 = 11.2_wp 86 REAL(wp) :: rm6 = 0.4_wp 87 REAL(wp) :: rm7 = 0.0_wp 88 REAL(wp) :: rm8 = 0.318_wp 89 REAL(wp) :: rc02 90 REAL(wp) :: rc02r 91 REAL(wp) :: rc03 92 REAL(wp) :: rc04 93 REAL(wp) :: rc03_sqrt2_galp 94 REAL(wp) :: rsbc_mb 95 REAL(wp) :: rsbc_std 96 REAL(wp) :: rsbc_tke1 97 REAL(wp) :: rsbc_tke2 98 REAL(wp) :: rsbc_tke3 99 REAL(wp) :: rsbc_psi1 100 REAL(wp) :: rsbc_psi2 101 REAL(wp) :: rsbc_psi3 102 REAL(wp) :: rsbc_zs 103 REAL(wp) :: rfact_tke, rfact_psi 104 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff 105 REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 106 REAL(wp) :: rd0, rd1, rd2, rd3, rd4, rd5 107 REAL(wp) :: rnc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0 108 REAL(wp) :: rpsi3m, rpsi3p, rpp, rmm, rnn 109 109 110 110 !! * Substitutions … … 169 169 zty2 = 2. * (vtau(ji ,jj-1)*vmask(ji,jj-1,1) + vtau(ji,jj)*vmask(ji,jj,1)) / & 170 170 & MAX(1., vmask(ji,jj-1,1) + vmask(ji,jj,1)) 171 ustars2(ji,jj) = r n_sbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)171 ustars2(ji,jj) = rsbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 172 172 ! 173 173 ! bottom friction … … 183 183 ! Compute surface roughness length according to Charnock formula: 184 184 IF (ln_crban) THEN 185 zhsro(:,:) = MAX(r n_sbc_zs * ustars2(:,:), hsro)185 zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 186 186 ELSE 187 187 zhsro(:,:) = hsro … … 200 200 & / ( fse3vw_n(ji,jj,jk) & 201 201 & * fse3vw_b(ji,jj,jk) ) 202 eps(ji,jj,jk) = r n_c03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)202 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) 203 203 ENDDO 204 204 ENDDO … … 218 218 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbathy(ji,jj)) 219 219 zdown = vkarmn * fsdepw(ji,jj,jk) * (-fsdepw(ji,jj,jk)+fsdepw(ji,jj,mbathy(ji,jj))) 220 zwall (ji,jj,jk) = ( 1. + r n_e2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk)220 zwall (ji,jj,jk) = ( 1. + re2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk) 221 221 ENDDO 222 222 ENDDO … … 256 256 zdiss = dir*(diss/en(ji,jj,jk)) +(1-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 257 257 ! 258 ! Compute a wall function from 1. to r n_sc_psi*zwall/rn_sc_psi0258 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 259 259 ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) 260 260 ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 261 ! Otherwise, this should be r n_sc_psi/rn_sc_psi0261 ! Otherwise, this should be rsc_psi/rsc_psi0 262 262 IF (ln_sigpsi) THEN 263 263 zsigpsi = MIN(1., zesh2/eps(ji,jj,jk)) ! 0. <= zsigpsi <= 1. 264 zwall_psi(ji,jj,jk) = r n_sc_psi / (zsigpsi * rn_sc_psi + (1.-zsigpsi) * rn_sc_psi0 / MAX(zwall(ji,jj,jk),1.))264 zwall_psi(ji,jj,jk) = rsc_psi / (zsigpsi * rsc_psi + (1.-zsigpsi) * rsc_psi0 / MAX(zwall(ji,jj,jk),1.)) 265 265 ELSE 266 266 zwall_psi(ji,jj,jk) = 1. … … 268 268 ! 269 269 ! building the matrix 270 zcof = zfact_tke * tmask(ji,jj,jk)270 zcof = rfact_tke * tmask(ji,jj,jk) 271 271 ! 272 272 ! lower diagonal … … 294 294 DO jj = 2, jpjm1 295 295 DO ji = fs_2, fs_jpim1 ! vector opt. 296 zwall_psi(ji,jj,1) = r n_sc_psi/rn_sc_psi0296 zwall_psi(ji,jj,1) = rsc_psi/rsc_psi0 297 297 END DO 298 298 END DO … … 310 310 ! ! balance between the production and the dissipation terms including the wave effect 311 311 ! 312 en(:,:,1) = MAX( r n_sbc_tke1 * ustars2(:,:), rn_emin )312 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 313 313 z_elem_a(:,:,1) = en(:,:,1) 314 314 z_elem_c(:,:,1) = 0. … … 316 316 ! 317 317 ! one level below 318 en(:,:,2) = MAX( r n_sbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**rn_a_sf, rn_emin )318 en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 319 319 z_elem_a(:,:,2) = 0. 320 320 z_elem_c(:,:,2) = 0. … … 325 325 ! ! balance between the production and the dissipation terms 326 326 ! 327 en(:,:,1) = MAX( r n_c02r * ustars2(:,:), rn_emin )327 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 328 328 z_elem_a(:,:,1) = en(:,:,1) 329 329 z_elem_c(:,:,1) = 0. … … 331 331 ! 332 332 ! one level below 333 en(:,:,2) = MAX( r n_c02r * ustars2(:,:), rn_emin )333 en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 334 334 z_elem_a(:,:,2) = 0. 335 335 z_elem_c(:,:,2) = 0. … … 343 343 ! 344 344 ! Dirichlet conditions at k=1 (Cosmetic) 345 en(:,:,1) = MAX( r n_sbc_tke1 * ustars2(:,:), rn_emin )345 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 346 346 z_elem_a(:,:,1) = en(:,:,1) 347 347 z_elem_c(:,:,1) = 0. … … 350 350 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 351 351 z_elem_a(:,:,2) = 0. 352 zflxs(:,:) = r n_sbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*rn_a_sf)352 zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*ra_sf) 353 353 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 354 354 ! … … 356 356 ! 357 357 ! Dirichlet conditions at k=1 (Cosmetic) 358 en(:,:,1) = MAX( r n_c02r * ustars2(:,:), rn_emin )358 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 359 359 z_elem_a(:,:,1) = en(:,:,1) 360 360 z_elem_c(:,:,1) = 0. … … 388 388 z_elem_c(ji,jj,ibot ) = 0. 389 389 z_elem_b(ji,jj,ibot ) = 1. 390 en(ji,jj,ibot ) = MAX( r n_c02r * ustarb2(ji,jj), rn_emin )390 en(ji,jj,ibot ) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 391 391 ! 392 392 ! Just above last level, Dirichlet condition again … … 394 394 z_elem_c(ji,jj,ibotm1) = 0. 395 395 z_elem_b(ji,jj,ibotm1) = 1. 396 en(ji,jj,ibotm1) = MAX( r n_c02r * ustarb2(ji,jj), rn_emin )396 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 397 397 END DO 398 398 END DO … … 411 411 z_elem_c(ji,jj,ibot) = 0. 412 412 z_elem_b(ji,jj,ibot) = 1. 413 en(ji,jj,ibot) = MAX( r n_c02r * ustarb2(ji,jj), rn_emin )413 en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 414 414 ! 415 415 ! Just above last level: Neumann condition … … 482 482 DO jj = 2, jpjm1 483 483 DO ji = fs_2, fs_jpim1 ! vector opt. 484 psi(ji,jj,jk) = SQRT(en(ji,jj,jk)) / (r n_c0 * mxln(ji,jj,jk))484 psi(ji,jj,jk) = SQRT(en(ji,jj,jk)) / (rc0 * mxln(ji,jj,jk)) 485 485 ENDDO 486 486 ENDDO … … 492 492 DO jj = 2, jpjm1 493 493 DO ji = fs_2, fs_jpim1 ! vector opt. 494 psi(ji,jj,jk) = r n_c02 * en(ji,jj,jk) * mxln(ji,jj,jk)**znn494 psi(ji,jj,jk) = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 495 495 ENDDO 496 496 ENDDO … … 516 516 dir = 0.5 + sign(0.5,rn2(ji,jj,jk)) 517 517 ! 518 r n_psi3 = dir*rn_psi3m+(1-dir)*rn_psi3p518 rpsi3 = dir*rpsi3m+(1-dir)*rpsi3p 519 519 ! 520 520 ! shear prod. - stratif. destruction 521 prod = r n_psi1 * zratio * shear(ji,jj,jk)521 prod = rpsi1 * zratio * shear(ji,jj,jk) 522 522 ! 523 523 ! stratif. destruction 524 buoy = r n_psi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk))524 buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk)) 525 525 ! 526 526 ! shear prod. - stratif. destruction 527 diss = r n_psi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)527 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 528 528 ! 529 529 dir = 0.5 + sign(0.5,prod+buoy) ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) … … 533 533 ! 534 534 ! building the matrix 535 zcof = zfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)535 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 536 536 ! lower diagonal 537 537 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & … … 563 563 ! ! balance between the production and the dissipation terms including the wave effect 564 564 ! 565 zdep(:,:) = r n_l_sf * zhsro(:,:)566 psi (:,:,1) = r n_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)565 zdep(:,:) = rl_sf * zhsro(:,:) 566 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 567 567 z_elem_a(:,:,1) = psi(:,:,1) 568 568 z_elem_c(:,:,1) = 0. … … 570 570 ! 571 571 ! one level below 572 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**( zmm*rn_a_sf+znn) ) &573 & / zhsro(:,:)**( zmm*rn_a_sf)574 psi (:,:,2) = r n_sbc_psi1 * ustars2(:,:)**zmm * zdep(:,:) * tmask(:,:,1)572 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**(rmm*ra_sf+rnn) ) & 573 & / zhsro(:,:)**(rmm*ra_sf) 574 psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 575 575 z_elem_a(:,:,2) = 0. 576 576 z_elem_c(:,:,2) = 0. … … 582 582 ! 583 583 zdep(:,:) = vkarmn * zhsro(:,:) 584 psi (:,:,1) = r n_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)584 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 585 585 z_elem_a(:,:,1) = psi(:,:,1) 586 586 z_elem_c(:,:,1) = 0. … … 589 589 ! one level below 590 590 zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 591 psi (:,:,2) = r n_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)591 psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 592 592 z_elem_a(:,:,2) = 0. 593 593 z_elem_c(:,:,2) = 0. … … 600 600 IF (ln_crban) THEN ! Wave induced mixing case 601 601 ! 602 zdep(:,:) = r n_l_sf * zhsro(:,:)603 psi (:,:,1) = r n_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)602 zdep(:,:) = rl_sf * zhsro(:,:) 603 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 604 604 z_elem_a(:,:,1) = psi(:,:,1) 605 605 z_elem_c(:,:,1) = 0. … … 611 611 ! 612 612 ! Set psi vertical flux at the surface: 613 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**( zmm*rn_a_sf+znn-1.) / zhsro(:,:)**(zmm*rn_a_sf)614 zflxs(:,:) = r n_sbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &615 & * en(:,:,1)** zmm * zdep613 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1.) / zhsro(:,:)**(rmm*ra_sf) 614 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 615 & * en(:,:,1)**rmm * zdep 616 616 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 617 617 ! … … 619 619 ! 620 620 zdep(:,:) = vkarmn * zhsro(:,:) 621 psi (:,:,1) = r n_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)621 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 622 622 z_elem_a(:,:,1) = psi(:,:,1) 623 623 z_elem_c(:,:,1) = 0. … … 630 630 ! Set psi vertical flux at the surface: 631 631 zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) 632 zflxs(:,:) = r n_sbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**zmm * zdep**(znn-1.)632 zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1.) 633 633 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 634 634 ! … … 653 653 ibotm1 = ibot-1 654 654 zdep(ji,jj) = vkarmn * rn_hbro 655 psi (ji,jj,ibot) = r n_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn655 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 656 656 z_elem_a(ji,jj,ibot) = 0. 657 657 z_elem_c(ji,jj,ibot) = 0. … … 660 660 ! Just above last level, Dirichlet condition again (GOTM like) 661 661 zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1)) 662 psi (ji,jj,ibotm1) = r n_c0**zpp * en(ji,jj,ibot )**zmm * zdep(ji,jj)**znn662 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 663 663 z_elem_a(ji,jj,ibotm1) = 0. 664 664 z_elem_c(ji,jj,ibotm1) = 0. … … 679 679 ! Bottom level Dirichlet condition: 680 680 zdep(ji,jj) = vkarmn * rn_hbro 681 psi (ji,jj,ibot) = r n_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn681 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 682 682 ! 683 683 z_elem_a(ji,jj,ibot) = 0. … … 691 691 ! Set psi vertical flux at the bottom: 692 692 zdep(ji,jj) = rn_hbro + 0.5*fse3t(ji,jj,ibotm1) 693 zflxb = r n_sbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * &694 & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))** zmm * zdep(ji,jj)**(znn-1.)693 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * & 694 & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1.) 695 695 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1) 696 696 END DO … … 734 734 DO jj = 2, jpjm1 735 735 DO ji = fs_2, fs_jpim1 ! vector opt. 736 eps(ji,jj,jk) = r n_c03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk)736 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk) 737 737 ENDDO 738 738 ENDDO … … 754 754 DO jj = 2, jpjm1 755 755 DO ji = fs_2, fs_jpim1 ! vector opt. 756 eps(ji,jj,jk) = r n_c04 * en(ji,jj,jk) * psi(ji,jj,jk)756 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 757 757 ENDDO 758 758 ENDDO … … 764 764 DO jj = 2, jpjm1 765 765 DO ji = fs_2, fs_jpim1 ! vector opt. 766 eps(ji,jj,jk) = r n_c0**(3.+zpp/znn) * en(ji,jj,jk)**(1.5+zmm/znn) * psi(ji,jj,jk)**(-1./znn)766 eps(ji,jj,jk) = rc0**(3.+rpp/rnn) * en(ji,jj,jk)**(1.5+rmm/rnn) * psi(ji,jj,jk)**(-1./rnn) 767 767 ENDDO 768 768 ENDDO … … 778 778 ! limitation 779 779 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 780 mxln(ji,jj,jk) = r n_c03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / eps(ji,jj,jk)780 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / eps(ji,jj,jk) 781 781 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 782 782 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) … … 801 801 ! Gh = -N²l²/q² 802 802 gh = - rn2(ji,jj,jk) * zcof 803 ! gh = MIN(gh, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri)) 804 gh = MIN( gh, rn_gh0 ) 805 gh = MAX( gh, rn_ghmin ) 803 gh = MIN( gh, rgh0 ) 804 gh = MAX( gh, rghmin ) 806 805 ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 807 sh = r n_a2*( 1.-6.*rn_a1/rn_b1 ) / ( 1.-3.*rn_a2*gh*(6.*rn_a1+rn_b2*( 1.-rn_c3 ) ) )808 sm = ( r n_b1**(-1./3.) + ( 18*rn_a1*rn_a1 + 9.*rn_a1*rn_a2*(1.-rn_c2) )*sh*gh ) / (1.-9.*rn_a1*rn_a2*gh)806 sh = ra2*( 1.-6.*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1.-rc3 ) ) ) 807 sm = ( rb1**(-1./3.) + ( 18*ra1*ra1 + 9.*ra1*ra2*(1.-rc2) )*sh*gh ) / (1.-9.*ra1*ra2*gh) 809 808 ! 810 809 ! Store stability function in avmu and avmv 811 avmu(ji,jj,jk) = r n_c_diff * sh * tmask(ji,jj,jk)812 avmv(ji,jj,jk) = r n_c_diff * sm * tmask(ji,jj,jk)810 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 811 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 813 812 END DO 814 813 END DO … … 824 823 ! Gh = -N²l²/q² 825 824 gh = - rn2(ji,jj,jk) * zcof 826 ! gh = MIN(gh, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri)) 827 gh = MIN( gh, rn_gh0 ) 828 gh = MAX( gh, rn_ghmin ) 829 gh = gh*rn_f6 825 gh = MIN( gh, rgh0 ) 826 gh = MAX( gh, rghmin ) 827 gh = gh*rf6 830 828 ! Gm = M²l²/q² Shear number 831 829 shr = shear(ji,jj,jk)/MAX(avm(ji,jj,jk), rsmall) 832 830 gm = shr * zcof 833 831 gm = MAX(gm, 1.e-10) 834 gm = gm*r n_f6835 gm = MIN ( (r n_d0 - rn_d1*gh + rn_d3*gh**2.)/(rn_d2-rn_d4*gh) , gm )832 gm = gm*rf6 833 gm = MIN ( (rd0 - rd1*gh + rd3*gh**2.)/(rd2-rd4*gh) , gm ) 836 834 ! Stability functions from Canuto 837 r n_cff = rn_d0 - rn_d1*gh +rn_d2*gm + rn_d3*gh**2. - rn_d4*gh*gm + rn_d5*gm**2.838 sm = (r n_s0 - rn_s1*gh + rn_s2*gm) / rn_cff839 sh = (r n_s4 - rn_s5*gh + rn_s6*gm) / rn_cff835 rcff = rd0 - rd1*gh +rd2*gm + rd3*gh**2. - rd4*gh*gm + rd5*gm**2. 836 sm = (rs0 - rs1*gh + rs2*gm) / rcff 837 sh = (rs4 - rs5*gh + rs6*gm) / rcff 840 838 ! 841 839 ! Store stability function in avmu and avmv 842 avmu(ji,jj,jk) = r n_c_diff * sh * tmask(ji,jj,jk)843 avmv(ji,jj,jk) = r n_c_diff * sm * tmask(ji,jj,jk)840 avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 841 avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 844 842 END DO 845 843 END DO … … 852 850 DO jj = 2, jpjm1 853 851 DO ji = fs_2, fs_jpim1 ! vector opt. 854 avmv(ji,jj,1) = r n_cm_sf / SQRT(2.)852 avmv(ji,jj,1) = rcm_sf / SQRT(2.) 855 853 END DO 856 854 END DO … … 859 857 ibot=mbathy(ji,jj) 860 858 ibotm1=ibot-1 861 avmv(ji,jj,ibot) = r n_c0 / SQRT(2.)859 avmv(ji,jj,ibot) = rc0 / SQRT(2.) 862 860 END DO 863 861 END DO … … 988 986 ! 989 987 IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada' 990 zpp = 0.991 zmm = 1.992 znn = 1.993 r n_sc_tke = 1.96994 r n_sc_psi = 1.96995 r n_psi1 = 0.9996 r n_psi3p = 1.997 r n_psi2 = 0.5988 rpp = 0. 989 rmm = 1. 990 rnn = 1. 991 rsc_tke = 1.96 992 rsc_psi = 1.96 993 rpsi1 = 0.9 994 rpsi3p = 1. 995 rpsi2 = 0.5 998 996 ! 999 997 SELECT CASE ( nn_stab_func ) 1000 998 ! 1001 999 CASE( 0, 1 ) ! G88 or KC stability functions 1002 r n_psi3m = 2.531000 rpsi3m = 2.53 1003 1001 CASE( 2 ) ! Canuto A stability functions 1004 r n_psi3m = 2.381002 rpsi3m = 2.38 1005 1003 CASE( 3 ) ! Canuto B stability functions 1006 r n_psi3m = 2.38 ! caution : constant not identified1004 rpsi3m = 2.38 ! caution : constant not identified 1007 1005 END SELECT 1008 1006 ! … … 1010 1008 ! 1011 1009 IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps' 1012 zpp = 3.1013 zmm = 1.51014 znn = -1.1015 r n_sc_tke = 1.1016 r n_sc_psi = 1.3 ! Schmidt number for psi1017 r n_psi1 = 1.441018 r n_psi3p = 1.1019 r n_psi2 = 1.921010 rpp = 3. 1011 rmm = 1.5 1012 rnn = -1. 1013 rsc_tke = 1. 1014 rsc_psi = 1.3 ! Schmidt number for psi 1015 rpsi1 = 1.44 1016 rpsi3p = 1. 1017 rpsi2 = 1.92 1020 1018 ! 1021 1019 SELECT CASE ( nn_stab_func ) 1022 1020 ! 1023 1021 CASE( 0, 1 ) ! G88 or KC stability functions 1024 r n_psi3m = -0.521022 rpsi3m = -0.52 1025 1023 CASE( 2 ) ! Canuto A stability functions 1026 r n_psi3m = -0.6291024 rpsi3m = -0.629 1027 1025 CASE( 3 ) ! Canuto B stability functions 1028 r n_psi3m = -0.5661026 rpsi3m = -0.566 1029 1027 END SELECT 1030 1028 ! … … 1032 1030 ! 1033 1031 IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega' 1034 zpp = -1.1035 zmm = 0.51036 znn = -1.1037 r n_sc_tke = 2.1038 r n_sc_psi = 2.1039 r n_psi1 = 0.5551040 r n_psi3p = 1.1041 r n_psi2 = 0.8331032 rpp = -1. 1033 rmm = 0.5 1034 rnn = -1. 1035 rsc_tke = 2. 1036 rsc_psi = 2. 1037 rpsi1 = 0.555 1038 rpsi3p = 1. 1039 rpsi2 = 0.833 1042 1040 ! 1043 1041 SELECT CASE ( nn_stab_func ) 1044 1042 ! 1045 1043 CASE( 0, 1 ) ! G88 or KC stability functions 1046 r n_psi3m = -0.581044 rpsi3m = -0.58 1047 1045 CASE( 2 ) ! Canuto A stability functions 1048 r n_psi3m = -0.641046 rpsi3m = -0.64 1049 1047 CASE( 3 ) ! Canuto B stability functions 1050 r n_psi3m = -0.64 ! caution : constant not identified1048 rpsi3m = -0.64 ! caution : constant not identified 1051 1049 END SELECT 1052 1050 ! … … 1054 1052 ! 1055 1053 IF(lwp) WRITE(numout,*) 'The choosen closure is generic' 1056 zpp = 2.1057 zmm = 1.1058 znn = -0.671059 r n_sc_tke = 0.81060 r n_sc_psi = 1.071061 r n_psi1 = 1.1062 r n_psi3p = 1.1063 r n_psi2 = 1.221054 rpp = 2. 1055 rmm = 1. 1056 rnn = -0.67 1057 rsc_tke = 0.8 1058 rsc_psi = 1.07 1059 rpsi1 = 1. 1060 rpsi3p = 1. 1061 rpsi2 = 1.22 1064 1062 ! 1065 1063 SELECT CASE ( nn_stab_func ) 1066 1064 ! 1067 1065 CASE( 0, 1 ) ! G88 or KC stability functions 1068 r n_psi3m = 0.11066 rpsi3m = 0.1 1069 1067 CASE( 2 ) ! Canuto A stability functions 1070 r n_psi3m = 0.051068 rpsi3m = 0.05 1071 1069 CASE( 3 ) ! Canuto B stability functions 1072 r n_psi3m = 0.05 ! caution : constant not identified1070 rpsi3m = 0.05 ! caution : constant not identified 1073 1071 END SELECT 1074 1072 ! … … 1083 1081 ! 1084 1082 IF(lwp) WRITE(numout,*) 'Stability functions from Galperin' 1085 r n_c2 = 0.1086 r n_c3 = 0.1087 r n_c_diff = 1.1088 r n_c0 = 0.55441089 r n_cm_sf = 0.98841090 r n_ghmin = -0.281091 r n_gh0 = 0.02331092 r n_ghcri = 0.021083 rc2 = 0. 1084 rc3 = 0. 1085 rc_diff = 1. 1086 rc0 = 0.5544 1087 rcm_sf = 0.9884 1088 rghmin = -0.28 1089 rgh0 = 0.0233 1090 rghcri = 0.02 1093 1091 ! 1094 1092 CASE ( 1 ) ! Kantha-Clayson stability functions 1095 1093 ! 1096 1094 IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson' 1097 r n_c2 = 0.71098 r n_c3 = 0.21099 r n_c_diff = 1.1100 r n_c0 = 0.55441101 r n_cm_sf = 0.98841102 r n_ghmin = -0.281103 r n_gh0 = 0.02331104 r n_ghcri = 0.021095 rc2 = 0.7 1096 rc3 = 0.2 1097 rc_diff = 1. 1098 rc0 = 0.5544 1099 rcm_sf = 0.9884 1100 rghmin = -0.28 1101 rgh0 = 0.0233 1102 rghcri = 0.02 1105 1103 ! 1106 1104 CASE ( 2 ) ! Canuto A stability functions 1107 1105 ! 1108 1106 IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A' 1109 r n_s0 = 1.5*rn_l1*rn_l5**2.1110 r n_s1 = -rn_l4*(rn_l6+rn_l7) + 2.*rn_l4*rn_l5*(rn_l1-(1./3.)*rn_l2-rn_l3)+1.5*rn_l1*rn_l5*rn_l81111 r n_s2 = -(3./8.)*rn_l1*(rn_l6**2.-rn_l7**2.)1112 r n_s4 = 2.*rn_l51113 r n_s5 = 2.*rn_l41114 r n_s6 = (2./3.)*rn_l5*(3.*rn_l3**2.-rn_l2**2.)-0.5*rn_l5*rn_l1*(3.*rn_l3-rn_l2)+0.75*rn_l1*(rn_l6-rn_l7)1115 r n_d0 = 3*rn_l5**2.1116 r n_d1 = rn_l5*(7.*rn_l4+3.*rn_l8)1117 r n_d2 = rn_l5**2.*(3.*rn_l3**2.-rn_l2**2.)-0.75*(rn_l6**2.-rn_l7**2.)1118 r n_d3 = rn_l4*(4.*rn_l4+3.*rn_l8)1119 r n_d4 = rn_l4*(rn_l2*rn_l6-3.*rn_l3*rn_l7-rn_l5*(rn_l2**2.-rn_l3**2.))+rn_l5*rn_l8*(3.*rn_l3**2.-rn_l2**2.)1120 r n_d5 = 0.25*(rn_l2**2.-3.*rn_l3**2.)*(rn_l6**2.-rn_l7**2.)1121 r n_c0 = 0.52681122 r n_f6 = 8. / (rn_c0**6.)1123 r n_c_diff = SQRT(2.)/(rn_c0**3.)1124 r n_cm_sf = 0.73101125 r n_ghmin = -0.281126 r n_gh0 = 0.03291127 r n_ghcri = 0.031107 rs0 = 1.5*rl1*rl5**2. 1108 rs1 = -rl4*(rl6+rl7) + 2.*rl4*rl5*(rl1-(1./3.)*rl2-rl3)+1.5*rl1*rl5*rl8 1109 rs2 = -(3./8.)*rl1*(rl6**2.-rl7**2.) 1110 rs4 = 2.*rl5 1111 rs5 = 2.*rl4 1112 rs6 = (2./3.)*rl5*(3.*rl3**2.-rl2**2.)-0.5*rl5*rl1*(3.*rl3-rl2)+0.75*rl1*(rl6-rl7) 1113 rd0 = 3*rl5**2. 1114 rd1 = rl5*(7.*rl4+3.*rl8) 1115 rd2 = rl5**2.*(3.*rl3**2.-rl2**2.)-0.75*(rl6**2.-rl7**2.) 1116 rd3 = rl4*(4.*rl4+3.*rl8) 1117 rd4 = rl4*(rl2*rl6-3.*rl3*rl7-rl5*(rl2**2.-rl3**2.))+rl5*rl8*(3.*rl3**2.-rl2**2.) 1118 rd5 = 0.25*(rl2**2.-3.*rl3**2.)*(rl6**2.-rl7**2.) 1119 rc0 = 0.5268 1120 rf6 = 8. / (rc0**6.) 1121 rc_diff = SQRT(2.)/(rc0**3.) 1122 rcm_sf = 0.7310 1123 rghmin = -0.28 1124 rgh0 = 0.0329 1125 rghcri = 0.03 1128 1126 ! 1129 1127 CASE ( 3 ) ! Canuto B stability functions 1130 1128 ! 1131 1129 IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B' 1132 r n_s0 = 1.5*rn_m1*rn_m5**2.1133 r n_s1 = -rn_m4*(rn_m6+rn_m7) + 2.*rn_m4*rn_m5*(rn_m1-(1./3.)*rn_m2-rn_m3)+1.5*rn_m1*rn_m5*rn_m81134 r n_s2 = -(3./8.)*rn_m1*(rn_m6**2.-rn_m7**2.)1135 r n_s4 = 2.*rn_m51136 r n_s5 = 2.*rn_m41137 r n_s6 = (2./3.)*rn_m5*(3.*rn_m3**2.-rn_m2**2.)-0.5*rn_m5*rn_m1*(3.*rn_m3-rn_m2)+0.75*rn_m1*(rn_m6-rn_m7)1138 r n_d0 = 3*rn_m5**2.1139 r n_d1 = rn_m5*(7.*rn_m4+3.*rn_m8)1140 r n_d2 = rn_m5**2.*(3.*rn_m3**2.-rn_m2**2.)-0.75*(rn_m6**2.-rn_m7**2.)1141 r n_d3 = rn_m4*(4.*rn_m4+3.*rn_m8)1142 r n_d4 = rn_m4*(rn_m2*rn_m6-3.*rn_m3*rn_m7-rn_m5*(rn_m2**2.-rn_m3**2.))+rn_m5*rn_m8*(3.*rn_m3**2.-rn_m2**2.)1143 r n_d5 = 0.25*(rn_m2**2.-3.*rn_m3**2.)*(rn_m6**2.-rn_m7**2.)1144 r n_c0 = 0.5268 !! rn_c0 = 0.5540 (Warner ...) to verify !1145 r n_f6 = 8. / (rn_c0**6.)1146 r n_c_diff = SQRT(2.)/(rn_c0**3.)1147 r n_cm_sf = 0.74701148 r n_ghmin = -0.281149 r n_gh0 = 0.04441150 r n_ghcri = 0.04141130 rs0 = 1.5*rm1*rm5**2. 1131 rs1 = -rm4*(rm6+rm7) + 2.*rm4*rm5*(rm1-(1./3.)*rm2-rm3)+1.5*rm1*rm5*rm8 1132 rs2 = -(3./8.)*rm1*(rm6**2.-rm7**2.) 1133 rs4 = 2.*rm5 1134 rs5 = 2.*rm4 1135 rs6 = (2./3.)*rm5*(3.*rm3**2.-rm2**2.)-0.5*rm5*rm1*(3.*rm3-rm2)+0.75*rm1*(rm6-rm7) 1136 rd0 = 3*rm5**2. 1137 rd1 = rm5*(7.*rm4+3.*rm8) 1138 rd2 = rm5**2.*(3.*rm3**2.-rm2**2.)-0.75*(rm6**2.-rm7**2.) 1139 rd3 = rm4*(4.*rm4+3.*rm8) 1140 rd4 = rm4*(rm2*rm6-3.*rm3*rm7-rm5*(rm2**2.-rm3**2.))+rm5*rm8*(3.*rm3**2.-rm2**2.) 1141 rd5 = 0.25*(rm2**2.-3.*rm3**2.)*(rm6**2.-rm7**2.) 1142 rc0 = 0.5268 !! rc0 = 0.5540 (Warner ...) to verify ! 1143 rf6 = 8. / (rc0**6.) 1144 rc_diff = SQRT(2.)/(rc0**3.) 1145 rcm_sf = 0.7470 1146 rghmin = -0.28 1147 rgh0 = 0.0444 1148 rghcri = 0.0414 1151 1149 ! 1152 1150 END SELECT … … 1157 1155 ! or equation (17) of Burchard, JPO, 31, 3133-3145, 2001 1158 1156 IF ((ln_sigpsi).AND.(ln_crban)) THEN 1159 zcr = SQRT(1.5*r n_sc_tke) * rn_cm_sf /vkarmn1160 r n_sc_psi0 = vkarmn**2/(rn_psi2*rn_cm_sf**2) * &1161 & ( znn**2 - 4./3.*zcr*znn*zmm - 1./3.*zcr*znn &1162 & + 2./9.* zmm*zcr**2 + 4./9.*zcr**2*zmm**2)1157 zcr = SQRT(1.5*rsc_tke) * rcm_sf /vkarmn 1158 rsc_psi0 = vkarmn**2/(rpsi2*rcm_sf**2) * & 1159 & ( rnn**2 - 4./3.*zcr*rnn*rmm - 1./3.*zcr*rnn & 1160 & + 2./9.*rmm*zcr**2 + 4./9.*zcr**2*rmm**2) 1163 1161 ELSE 1164 r n_sc_psi0 = rn_sc_psi1162 rsc_psi0 = rsc_psi 1165 1163 ENDIF 1166 1164 1167 1165 ! Shear free turbulence parameters: 1168 1166 ! 1169 r n_a_sf = -4.*znn*SQRT(rn_sc_tke) / ( (1.+4.*zmm)*SQRT(rn_sc_tke) &1170 & - SQRT(r n_sc_tke + 24.*rn_sc_psi0*rn_psi2 ) )1171 r n_l_sf = rn_c0 * SQRT(rn_c0/rn_cm_sf) * SQRT( ( (1. + 4.*zmm + 8.*zmm**2)*rn_sc_tke &1172 & + 12. * r n_sc_psi0*rn_psi2 - (1. + 4.*zmm) &1173 & *SQRT(r n_sc_tke*(rn_sc_tke &1174 & + 24.*r n_sc_psi0*rn_psi2)) ) &1175 & /(12.* znn**2.) &1167 ra_sf = -4.*rnn*SQRT(rsc_tke) / ( (1.+4.*rmm)*SQRT(rsc_tke) & 1168 & - SQRT(rsc_tke + 24.*rsc_psi0*rpsi2 ) ) 1169 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1. + 4.*rmm + 8.*rmm**2)*rsc_tke & 1170 & + 12. * rsc_psi0*rpsi2 - (1. + 4.*rmm) & 1171 & *SQRT(rsc_tke*(rsc_tke & 1172 & + 24.*rsc_psi0*rpsi2)) ) & 1173 & /(12.*rnn**2.) & 1176 1174 & ) 1177 1175 … … 1182 1180 WRITE(numout,*) 'Limit values' 1183 1181 WRITE(numout,*) '~~~~~~~~~~~~' 1184 WRITE(numout,*) 'Parameter m = ', zmm1185 WRITE(numout,*) 'Parameter n = ', znn1186 WRITE(numout,*) 'Parameter p = ', zpp1187 WRITE(numout,*) 'r n_psi1 = ',rn_psi11188 WRITE(numout,*) 'r n_psi2 = ',rn_psi21189 WRITE(numout,*) 'r n_psi3m = ',rn_psi3m1190 WRITE(numout,*) 'r n_psi3p = ',rn_psi3p1191 WRITE(numout,*) 'r n_sc_tke = ',rn_sc_tke1192 WRITE(numout,*) 'r n_sc_psi = ',rn_sc_psi1193 WRITE(numout,*) 'r n_sc_psi0 = ',rn_sc_psi01194 WRITE(numout,*) 'r n_c0 = ',rn_c01182 WRITE(numout,*) 'Parameter m = ',rmm 1183 WRITE(numout,*) 'Parameter n = ',rnn 1184 WRITE(numout,*) 'Parameter p = ',rpp 1185 WRITE(numout,*) 'rpsi1 = ',rpsi1 1186 WRITE(numout,*) 'rpsi2 = ',rpsi2 1187 WRITE(numout,*) 'rpsi3m = ',rpsi3m 1188 WRITE(numout,*) 'rpsi3p = ',rpsi3p 1189 WRITE(numout,*) 'rsc_tke = ',rsc_tke 1190 WRITE(numout,*) 'rsc_psi = ',rsc_psi 1191 WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0 1192 WRITE(numout,*) 'rc0 = ',rc0 1195 1193 WRITE(numout,*) 1196 1194 WRITE(numout,*) 'Shear free turbulence parameters:' 1197 WRITE(numout,*) 'r n_cm_sf = ',rn_cm_sf1198 WRITE(numout,*) 'r n_a_sf = ',rn_a_sf1199 WRITE(numout,*) 'r n_l_sf = ',rn_l_sf1195 WRITE(numout,*) 'rcm_sf = ',rcm_sf 1196 WRITE(numout,*) 'ra_sf = ',ra_sf 1197 WRITE(numout,*) 'rl_sf = ',rl_sf 1200 1198 WRITE(numout,*) 1201 1199 ENDIF 1202 1200 1203 1201 ! Constants initialization 1204 r n_c02r = 1. / rn_c0**2.1205 r n_c02 = rn_c0**2._wp1206 r n_c03 = rn_c0**3._wp1207 r n_c04 = rn_c0**4._wp1208 r n_c03_sqrt2_galp = rn_c03 / SQRT(2._wp) / rn_clim_galp1209 r n_sbc_mb = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg1210 r n_sbc_std = 3.75 ! Surf. bound. cond. standard (prod=diss)1211 r n_sbc_tke1 = (-rn_sc_tke*rn_crban/(rn_cm_sf*rn_a_sf*rn_l_sf))**(2./3.) ! k_eps = 53. Dirichlet + Wave breaking1212 r n_sbc_tke2 = 0.5 / rau01213 r n_sbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking1214 r n_sbc_zs = rn_charn/grav ! Charnock formula1215 r n_sbc_psi1 = rn_c0**zpp * rn_sbc_tke1**zmm * rn_l_sf**znn ! Dirichlet + Wave breaking1216 r n_sbc_psi2 = -0.5 * rdt * rn_c0**zpp * znn * vkarmn**znn / rn_sc_psi ! Neumann + NO Wave breaking1217 r n_sbc_psi3 = -0.5 * rdt * rn_c0**zpp * rn_l_sf**znn / rn_sc_psi * (znn + zmm*rn_a_sf) ! Neumann + Wave breaking1218 zfact_tke = -0.5 / rn_sc_tke * rdt ! Cst used for the Diffusion term of tke1219 zfact_psi = -0.5 / rn_sc_psi * rdt ! Cst used for the Diffusion term of tke1202 rc02r = 1. / rc0**2. 1203 rc02 = rc0**2._wp 1204 rc03 = rc0**3._wp 1205 rc04 = rc0**4._wp 1206 rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp 1207 rsbc_mb = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg 1208 rsbc_std = 3.75 ! Surf. bound. cond. standard (prod=diss) 1209 rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2./3.) ! k_eps = 53. Dirichlet + Wave breaking 1210 rsbc_tke2 = 0.5 / rau0 1211 rsbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking 1212 rsbc_zs = rn_charn/grav ! Charnock formula 1213 rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 1214 rsbc_psi2 = -0.5 * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1215 rsbc_psi3 = -0.5 * rdt * rc0**rpp * rl_sf**rnn / rsc_psi * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 1216 rfact_tke = -0.5 / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1217 rfact_psi = -0.5 / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1220 1218 1221 1219 ! Wall proximity function
Note: See TracChangeset
for help on using the changeset viewer.
