- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r4624 r6225 8 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 9 !!---------------------------------------------------------------------- 10 #if defined key_zdfgls || defined key_esopa10 #if defined key_zdfgls 11 11 !!---------------------------------------------------------------------- 12 12 !! 'key_zdfgls' Generic Length Scale vertical physics … … 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 … … 41 42 LOGICAL , PUBLIC, PARAMETER :: lk_zdfgls = .TRUE. !: TKE vertical mixing flag 42 43 ! 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy44 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mxln !: now mixing length 45 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k ! not enhanced Kz47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avm_k ! not enhanced Kz48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avmu_k ! not enhanced Kz49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avmv_k ! not enhanced Kz50 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 51 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points 52 48 53 49 ! !! ** Namelist namzdf_gls ** 54 LOGICAL :: ln_crban ! =T use Craig and Banner scheme55 50 LOGICAL :: ln_length_lim ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 56 51 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) 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 61 55 INTEGER :: nn_stab_func ! stability functions G88, KC or Canuto (=0/1/2) 62 56 INTEGER :: nn_clos ! closure 0/1/2/3 MY82/k-eps/k-w/gen … … 66 60 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 67 61 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) 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 71 65 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 72 66 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 … … 96 90 REAL(wp) :: rm7 = 0.0_wp 97 91 REAL(wp) :: rm8 = 0.318_wp 98 92 REAL(wp) :: rtrans = 0.1_wp 99 93 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 ! - - - - 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 ! - - - - 104 97 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - 105 98 REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 ! - - - - … … 109 102 110 103 !! * Substitutions 111 # include "domzgr_substitute.h90"112 104 # include "vectopt_loop_substitute.h90" 113 105 !!---------------------------------------------------------------------- … … 122 114 !! *** FUNCTION zdf_gls_alloc *** 123 115 !!---------------------------------------------------------------------- 124 ALLOCATE( en(jpi,jpj,jpk), mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 125 & avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk), & 126 & avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), & 127 & ustars2(jpi,jpj), ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 116 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 117 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 128 118 ! 129 119 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 147 137 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 148 138 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 139 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 149 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 150 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) … … 153 144 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 154 145 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 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 147 REAL(wp), POINTER, DIMENSION(:,:,:) :: psi ! psi at time now 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a ! element of the first matrix diagonal 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 157 151 !!-------------------------------------------------------------------- 158 152 ! 159 153 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 160 154 ! 161 CALL wrk_alloc( jpi,jpj, zdep, zflxs, 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 155 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 156 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 157 164 158 ! Preliminary computing 165 159 … … 174 168 175 169 ! 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 170 DO jj = 2, jpjm1 171 DO ji = fs_2, fs_jpim1 ! vector opt. 172 ! 173 ! surface friction 182 174 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 175 ! 176 ! bottom friction (explicit before friction) 177 ! Note that we chose here not to bound the friction as in dynbfr) 178 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 179 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 180 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 181 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 182 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 183 END DO 184 END DO 185 186 ! Set surface roughness length 187 SELECT CASE ( nn_z0_met ) 188 ! 189 CASE ( 0 ) ! Constant roughness 190 zhsro(:,:) = rn_hsro 191 CASE ( 1 ) ! Standard Charnock formula 192 zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 193 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 194 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 195 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 196 ! 197 END SELECT 199 198 200 199 ! Compute shear and dissipation rate … … 204 203 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 205 204 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 206 & / ( fse3uw_n(ji,jj,jk) & 207 & * fse3uw_b(ji,jj,jk) ) 205 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 208 206 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 209 207 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 210 & / ( fse3vw_n(ji,jj,jk) & 211 & * fse3vw_b(ji,jj,jk) ) 208 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 212 209 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) 213 210 END DO … … 226 223 DO jj = 2, jpjm1 227 224 DO ji = fs_2, fs_jpim1 ! vector opt. 228 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)229 zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbkt(ji,jj)+1) )225 zup = mxln(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 226 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 230 227 zcoef = ( zup / MAX( zdown, rsmall ) ) 231 228 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) … … 284 281 ! lower diagonal 285 282 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 286 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )283 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 287 284 ! 288 285 ! upper diagonal 289 286 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 290 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )287 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 291 288 ! 292 289 ! diagonal … … 303 300 ! 304 301 ! 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 302 zwall_psi(:,:,1) = zwall_psi(:,:,2) 303 zwall_psi(:,:,jpk) = 1. 304 ! 314 305 ! Surface boundary condition on tke 315 306 ! --------------------------------- 316 307 ! 317 SELECT CASE ( nn_ tkebc_surf )308 SELECT CASE ( nn_bc_surf ) 318 309 ! 319 310 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 ! 311 ! First level 312 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 313 en(:,:,1) = MAX(en(:,:,1), rn_emin) 314 z_elem_a(:,:,1) = en(:,:,1) 315 z_elem_c(:,:,1) = 0._wp 316 z_elem_b(:,:,1) = 1._wp 317 ! 318 ! One level below 319 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 320 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 321 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 322 z_elem_a(:,:,2) = 0._wp 323 z_elem_c(:,:,2) = 0._wp 324 z_elem_b(:,:,2) = 1._wp 325 ! 326 ! 351 327 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 ! 328 ! 329 ! Dirichlet conditions at k=1 330 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 331 en(:,:,1) = MAX(en(:,:,1), 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 ! at k=2, set de/dz=Fw 337 !cbr 338 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 339 z_elem_a(:,:,2) = 0._wp 340 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 341 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 342 & * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf) 343 344 en(:,:,2) = en(:,:,2) + zflxs(:,:)/e3w_n(:,:,2) 345 ! 346 ! 379 347 END SELECT 380 348 … … 382 350 ! -------------------------------- 383 351 ! 384 SELECT CASE ( nn_ tkebc_bot )352 SELECT CASE ( nn_bc_bot ) 385 353 ! 386 354 CASE ( 0 ) ! Dirichlet 387 355 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 388 356 ! ! Balance between the production and the dissipation terms 389 !CDIR NOVERRCHK 390 DO jj = 2, jpjm1 391 !CDIR NOVERRCHK 357 DO jj = 2, jpjm1 392 358 DO ji = fs_2, fs_jpim1 ! vector opt. 393 359 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 410 376 CASE ( 1 ) ! Neumman boundary condition 411 377 ! 412 !CDIR NOVERRCHK 413 DO jj = 2, jpjm1 414 !CDIR NOVERRCHK 378 DO jj = 2, jpjm1 415 379 DO ji = fs_2, fs_jpim1 ! vector opt. 416 380 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 457 421 ! ! set the minimum value of tke 458 422 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 459 423 460 424 !!----------------------------------------!! 461 425 !! Solve prognostic equation for psi !! … … 541 505 ! lower diagonal 542 506 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 543 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )507 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 544 508 ! upper diagonal 545 509 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 546 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )510 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 547 511 ! diagonal 548 512 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & … … 560 524 ! --------------------------------- 561 525 ! 562 SELECT CASE ( nn_ psibc_surf )526 SELECT CASE ( nn_bc_surf ) 563 527 ! 564 528 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 ! 529 ! 530 ! Surface value 531 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 532 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 533 z_elem_a(:,:,1) = psi(:,:,1) 534 z_elem_c(:,:,1) = 0._wp 535 z_elem_b(:,:,1) = 1._wp 536 ! 537 ! One level below 538 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 539 zdep(:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 540 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 541 z_elem_a(:,:,2) = 0._wp 542 z_elem_c(:,:,2) = 0._wp 543 z_elem_b(:,:,2) = 1._wp 544 ! 545 ! 603 546 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 ! 547 ! 548 ! Surface value: Dirichlet 549 zdep(:,:) = zhsro(:,:) * rl_sf 550 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 551 z_elem_a(:,:,1) = psi(:,:,1) 552 z_elem_c(:,:,1) = 0._wp 553 z_elem_b(:,:,1) = 1._wp 554 ! 555 ! Neumann condition at k=2 556 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 557 z_elem_a(:,:,2) = 0._wp 558 ! 559 ! Set psi vertical flux at the surface: 560 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 561 zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 562 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 563 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 564 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 565 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 566 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 567 568 ! 569 ! 642 570 END SELECT 643 571 … … 645 573 ! -------------------------------- 646 574 ! 647 SELECT CASE ( nn_psibc_bot ) 575 SELECT CASE ( nn_bc_bot ) 576 ! 648 577 ! 649 578 CASE ( 0 ) ! Dirichlet 650 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro579 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 651 580 ! ! Balance between the production and the dissipation terms 652 !CDIR NOVERRCHK 653 DO jj = 2, jpjm1 654 !CDIR NOVERRCHK 581 DO jj = 2, jpjm1 655 582 DO ji = fs_2, fs_jpim1 ! vector opt. 656 583 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 657 584 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 658 zdep(ji,jj) = vkarmn * hbro585 zdep(ji,jj) = vkarmn * rn_bfrz0 659 586 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 660 587 z_elem_a(ji,jj,ibot) = 0._wp … … 663 590 ! 664 591 ! Just above last level, Dirichlet condition again (GOTM like) 665 zdep(ji,jj) = vkarmn * ( hbro + fse3t(ji,jj,ibotm1) )592 zdep(ji,jj) = vkarmn * ( rn_bfrz0 + e3t_n(ji,jj,ibotm1) ) 666 593 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 667 594 z_elem_a(ji,jj,ibotm1) = 0._wp … … 673 600 CASE ( 1 ) ! Neumman boundary condition 674 601 ! 675 !CDIR NOVERRCHK 676 DO jj = 2, jpjm1 677 !CDIR NOVERRCHK 602 DO jj = 2, jpjm1 678 603 DO ji = fs_2, fs_jpim1 ! vector opt. 679 604 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 681 606 ! 682 607 ! Bottom level Dirichlet condition: 683 zdep(ji,jj) = vkarmn * hbro608 zdep(ji,jj) = vkarmn * rn_bfrz0 684 609 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 685 610 ! … … 693 618 ! 694 619 ! Set psi vertical flux at the bottom: 695 zdep(ji,jj) = hbro + 0.5_wp*fse3t(ji,jj,ibotm1)620 zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1) 696 621 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) & 697 622 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 698 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)623 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 699 624 END DO 700 625 END DO … … 736 661 DO jj = 2, jpjm1 737 662 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)663 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 664 END DO 740 665 END DO … … 783 708 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 784 709 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))710 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 711 END DO 787 712 END DO … … 847 772 ! Boundary conditions on stability functions for momentum (Neumann): 848 773 ! Lines below are useless if GOTM style Dirichlet conditions are used 849 zcoef = rcm_sf / SQRT( 2._wp ) 774 775 avmv(:,:,1) = avmv(:,:,2) 776 850 777 DO jj = 2, jpjm1 851 778 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 779 avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj)) 859 780 END DO 860 781 END DO … … 900 821 avmv_k(:,:,:) = avmv(:,:,:) 901 822 ! 902 CALL wrk_dealloc( jpi,jpj, zdep, zflxs, zhsro )903 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )823 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 824 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 904 825 ! 905 826 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls') … … 932 853 !! 933 854 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, & 855 & rn_clim_galp, ln_sigpsi, rn_hsro, & 856 & rn_crban, rn_charn, rn_frac_hs, & 857 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 938 858 & nn_stab_func, nn_clos 939 859 !!---------------------------------------------------------- … … 955 875 WRITE(numout,*) '~~~~~~~~~~~~' 956 876 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 877 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 878 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 879 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 880 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 881 WRITE(numout,*) ' TKE Surface boundary condition nn_bc_surf = ', nn_bc_surf 882 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 883 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 967 884 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban 968 885 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 886 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met 887 WRITE(numout,*) ' Wave height frac. (used if nn_z0_met=2) rn_frac_hs = ', rn_frac_hs 969 888 WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func 970 889 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 890 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 891 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 974 892 ENDIF 975 893 … … 978 896 979 897 ! !* 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' ) 898 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 899 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 900 IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 984 901 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 902 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) … … 1001 918 SELECT CASE ( nn_stab_func ) 1002 919 CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions 1003 CASE( 2 ) ; rpsi3m = 2. 38_wp ! Canuto A stability functions920 CASE( 2 ) ; rpsi3m = 2.62_wp ! Canuto A stability functions 1004 921 CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified) 1005 922 END SELECT … … 1012 929 rnn = -1._wp 1013 930 rsc_tke = 1._wp 1014 rsc_psi = 1. 3_wp ! Schmidt number for psi931 rsc_psi = 1.2_wp ! Schmidt number for psi 1015 932 rpsi1 = 1.44_wp 1016 933 rpsi3p = 1._wp … … 1140 1057 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1141 1058 ! ! 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 ) 1059 IF( ln_sigpsi ) THEN 1060 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1061 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1062 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work 1063 ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn 1064 rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.)) 1147 1065 ELSE 1148 1066 rsc_psi0 = rsc_psi … … 1151 1069 ! !* Shear free turbulence parameters 1152 1070 ! 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 ) ) 1071 ra_sf = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) & 1072 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1073 1074 IF ( rn_crban==0._wp ) THEN 1075 rl_sf = vkarmn 1076 ELSE 1077 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke & 1078 & + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) & 1079 & *SQRT(rsc_tke*(rsc_tke & 1080 & + 24._wp*rsc_psi0*rpsi2)) ) & 1081 & /(12._wp*rnn**2.) & 1082 & ) 1083 ENDIF 1160 1084 1161 1085 ! … … 1187 1111 rc03 = rc02 * rc0 1188 1112 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 1113 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1114 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1115 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1116 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1117 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1118 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1119 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1120 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1121 1122 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1123 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1201 1124 1202 1125 ! !* Wall proximity function … … 1257 1180 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1258 1181 en (:,:,:) = rn_emin 1259 mxln(:,:,:) = 0.001 1182 mxln(:,:,:) = 0.05 1183 avt_k (:,:,:) = avt (:,:,:) 1184 avm_k (:,:,:) = avm (:,:,:) 1185 avmu_k(:,:,:) = avmu(:,:,:) 1186 avmv_k(:,:,:) = avmv(:,:,:) 1260 1187 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1261 1188 ENDIF … … 1263 1190 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1264 1191 en (:,:,:) = rn_emin 1265 mxln(:,:,:) = 0.0 011192 mxln(:,:,:) = 0.05 1266 1193 ENDIF 1267 1194 ! … … 1269 1196 ! ! ------------------- 1270 1197 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1271 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1198 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1272 1199 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1273 1200 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1274 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1201 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1275 1202 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1276 1203 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln )
Note: See TracChangeset
for help on using the changeset viewer.
