Changeset 14834 for NEMO/trunk/src/OCE/ZDF/zdfgls.F90
- Timestamp:
- 2021-05-11T11:24:44+02:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ZDF/zdfgls.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdfgls.F90
r14156 r14834 137 137 USE zdf_oce , ONLY : en, avtb, avmb ! ocean vertical physics 138 138 !! 139 INTEGER , INTENT(in ) :: kt ! ocean time step140 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices141 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: p_sh2 ! shear production term142 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)139 INTEGER , INTENT(in ) :: kt ! ocean time step 140 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 141 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: p_sh2 ! shear production term 142 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 143 143 ! 144 144 INTEGER :: ji, jj, jk ! dummy loop arguments … … 151 151 REAL(wp) :: gh, gm, shr, dif, zsqen, zavt, zavm ! - - 152 152 REAL(wp) :: zmsku, zmskv ! - - 153 REAL(wp), DIMENSION( jpi,jpj) :: zdep154 REAL(wp), DIMENSION( jpi,jpj) :: zkar155 REAL(wp), DIMENSION( jpi,jpj) :: zflxs! Turbulence fluxed induced by internal waves156 REAL(wp), DIMENSION( jpi,jpj) :: zhsro! Surface roughness (surface waves)157 REAL(wp), DIMENSION( jpi,jpj) :: zice_fra! Tapering of wave breaking under sea ice158 REAL(wp), DIMENSION( jpi,jpj,jpk) :: eb! tke at time before159 REAL(wp), DIMENSION( jpi,jpj,jpk) :: hmxl_b! mixing length at time before160 REAL(wp), DIMENSION( jpi,jpj,jpk) :: eps! dissipation rate161 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwall_psi! Wall function use in the wb case (ln_sigpsi)162 REAL(wp), DIMENSION( jpi,jpj,jpk) :: psi! psi at time now163 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix164 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zstt, zstm! stability function on tracer and momentum153 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdep 154 REAL(wp), DIMENSION(A2D(nn_hls)) :: zkar 155 REAL(wp), DIMENSION(A2D(nn_hls)) :: zflxs ! Turbulence fluxed induced by internal waves 156 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhsro ! Surface roughness (surface waves) 157 REAL(wp), DIMENSION(A2D(nn_hls)) :: zice_fra ! Tapering of wave breaking under sea ice 158 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: eb ! tke at time before 159 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: hmxl_b ! mixing length at time before 160 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: eps ! dissipation rate 161 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 162 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: psi ! psi at time now 163 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zd_lw, zd_up, zdiag ! lower, upper and diagonal of the matrix 164 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zstt, zstm ! stability function on tracer and momentum 165 165 !!-------------------------------------------------------------------- 166 166 ! 167 167 ! Preliminary computing 168 169 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 170 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 171 ustar2_bot (:,:) = 0._wp 168 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 169 ustar2_surf(ji,jj) = 0._wp ; ustar2_top(ji,jj) = 0._wp ; ustar2_bot(ji,jj) = 0._wp 170 END_2D 171 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 172 psi(ji,jj,jk) = 0._wp ; zwall_psi(ji,jj,jk) = 0._wp 173 END_3D 172 174 173 175 SELECT CASE ( nn_z0_ice ) 174 176 CASE( 0 ) ; zice_fra(:,:) = 0._wp 175 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i( :,:) * 10._wp )176 CASE( 2 ) ; zice_fra(:,:) = fr_i( :,:)177 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i( :,:) , 1._wp )177 CASE( 1 ) ; zice_fra(:,:) = TANH( fr_i(A2D(nn_hls)) * 10._wp ) 178 CASE( 2 ) ; zice_fra(:,:) = fr_i(A2D(nn_hls)) 179 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp ) 178 180 END SELECT 179 181 180 182 ! Compute surface, top and bottom friction at T-points 181 DO_2D ( 0, 0, 0, 0) !== surface ocean friction183 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== surface ocean friction 182 184 ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) ! surface friction 183 185 END_2D … … 186 188 ! 187 189 IF( .NOT.ln_drg_OFF ) THEN !== top/bottom friction (explicit before friction) 188 DO_2D ( 0, 0, 0, 0 )! bottom friction (explicit before friction)190 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! bottom friction (explicit before friction) 189 191 zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 190 192 zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) … … 193 195 END_2D 194 196 IF( ln_isfcav ) THEN 195 DO_2D ( 0, 0, 0, 0) ! top friction197 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! top friction 196 198 zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 197 199 zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) … … 206 208 zhsro(:,:) = rn_hsro 207 209 CASE ( 1 ) ! Standard Charnock formula 208 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf( :,:) , rn_hsro )210 zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(A2D(nn_hls)) , rn_hsro ) 209 211 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 210 212 !!gm faster coding : the 2 comment lines should be used 211 213 !!gm zcof = 2._wp * 0.6_wp / 28._wp 212 214 !!gm zdep(:,:) = 30._wp * TANH( zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) ) ) ! Wave age (eq. 10) 213 zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) ) ! Wave age (eq. 10) 214 zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 215 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 216 zcof = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(ji,jj),rsmall))) ) ! Wave age (eq. 10) 217 zhsro(ji,jj) = MAX(rsbc_zs2 * ustar2_surf(ji,jj) * zcof**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 218 END_2D 215 219 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 216 zhsro(:,:) = MAX(rn_frac_hs * hsw( :,:), rn_hsro) ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )220 zhsro(:,:) = MAX(rn_frac_hs * hsw(A2D(nn_hls)), rn_hsro) ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 ) 217 221 END SELECT 218 222 ! 219 223 ! adapt roughness where there is sea ice 220 zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1) + (1._wp - tmask(:,:,1))*rn_hsro 221 ! 222 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== Compute dissipation rate ==! 224 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 225 zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1) + & 226 & (1._wp - tmask(ji,jj,1))*rn_hsro 227 END_2D 228 ! 229 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !== Compute dissipation rate ==! 223 230 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 224 231 END_3D 225 232 226 233 ! Save tke at before time step 227 eb (:,:,:) = en (:,:,:) 228 hmxl_b(:,:,:) = hmxl_n(:,:,:) 234 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 235 eb (ji,jj,jk) = en (ji,jj,jk) 236 hmxl_b(ji,jj,jk) = hmxl_n(ji,jj,jk) 237 END_3D 229 238 230 239 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 231 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )240 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 232 241 zup = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 233 242 zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) … … 250 259 ! Warning : after this step, en : right hand side of the matrix 251 260 252 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )261 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 253 262 ! 254 263 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction … … 303 312 ! 304 313 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 305 ! First level 306 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 ) 307 zd_lw(:,:,1) = en(:,:,1) 308 zd_up(:,:,1) = 0._wp 309 zdiag(:,:,1) = 1._wp 310 ! 311 ! One level below 312 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 313 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 314 zd_lw(:,:,2) = 0._wp 315 zd_up(:,:,2) = 0._wp 316 zdiag(:,:,2) = 1._wp 317 ! 318 ! 314 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 315 ! First level 316 en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3 ) 317 zd_lw(ji,jj,1) = en(ji,jj,1) 318 zd_up(ji,jj,1) = 0._wp 319 zdiag(ji,jj,1) = 1._wp 320 ! 321 ! One level below 322 en (ji,jj,2) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1 & 323 & * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )**(1.5_wp*ra_sf) )**r2_3 ) 324 zd_lw(ji,jj,2) = 0._wp 325 zd_up(ji,jj,2) = 0._wp 326 zdiag(ji,jj,2) = 1._wp 327 END_2D 328 ! 329 ! 319 330 CASE ( 1 ) ! Neumann boundary condition (set d(e)/dz) 320 ! 321 ! Dirichlet conditions at k=1 322 en (:,:,1) = MAX( rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin ) 323 zd_lw(:,:,1) = en(:,:,1) 324 zd_up(:,:,1) = 0._wp 325 zdiag(:,:,1) = 1._wp 326 ! 327 ! at k=2, set de/dz=Fw 328 !cbr 329 DO_2D( 0, 0, 0, 0 ) ! zdiag zd_lw not defined/used on the halo 330 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 331 zd_lw(ji,jj,2) = 0._wp 332 END_2D 333 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 334 zflxs(:,:) = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 335 & * ( ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:) )**(1.5_wp*ra_sf) 331 ! 332 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 333 ! Dirichlet conditions at k=1 334 en (ji,jj,1) = MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3 ) 335 zd_lw(ji,jj,1) = en(ji,jj,1) 336 zd_up(ji,jj,1) = 0._wp 337 zdiag(ji,jj,1) = 1._wp 338 ! 339 ! at k=2, set de/dz=Fw 340 !cbr 341 ! zdiag zd_lw not defined/used on the halo 342 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 343 zd_lw(ji,jj,2) = 0._wp 344 ! 345 zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) )) 346 zflxs(ji,jj) = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) & 347 & * ( ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj) )**(1.5_wp*ra_sf) 336 348 !!gm why not : * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 337 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 338 ! 339 ! 349 en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 350 END_2D 351 ! 352 ! 340 353 END SELECT 341 354 … … 348 361 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 349 362 ! ! Balance between the production and the dissipation terms 350 DO_2D ( 0, 0, 0, 0)363 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 351 364 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 352 365 !! With thick deep ocean level thickness, this may be quite large, no ??? … … 365 378 END_2D 366 379 ! 380 ! NOTE: ctl_stop with ln_isfcav when using GLS 367 381 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 368 DO_2D ( 0, 0, 0, 0)382 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 369 383 itop = mikt(ji,jj) ! k top w-point 370 384 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 384 398 CASE ( 1 ) ! Neumman boundary condition 385 399 ! 386 DO_2D ( 0, 0, 0, 0)400 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 387 401 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 388 402 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 398 412 en (ji,jj,ibot) = z_en 399 413 END_2D 414 ! NOTE: ctl_stop with ln_isfcav when using GLS 400 415 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 401 DO_2D ( 0, 0, 0, 0)416 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 402 417 itop = mikt(ji,jj) ! k top w-point 403 418 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one … … 420 435 ! ---------------------------------------------------------- 421 436 ! 422 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1437 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 423 438 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 424 439 END_3D 425 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1440 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 426 441 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 427 442 END_3D 428 DO_3DS ( 0, 0, 0, 0, jpkm1, 2, -1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk443 DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 429 444 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 430 445 END_3D 431 446 ! ! set the minimum value of tke 432 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 447 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 448 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) 449 END_3D 433 450 434 451 !!----------------------------------------!! … … 441 458 ! 442 459 CASE( 0 ) ! k-kl (Mellor-Yamada) 443 DO_3D( 0, 0, 0, 0, 2, jpkm1 )460 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 444 461 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 445 462 END_3D 446 463 ! 447 464 CASE( 1 ) ! k-eps 448 DO_3D( 0, 0, 0, 0, 2, jpkm1 )465 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 449 466 psi(ji,jj,jk) = eps(ji,jj,jk) 450 467 END_3D 451 468 ! 452 469 CASE( 2 ) ! k-w 453 DO_3D( 0, 0, 0, 0, 2, jpkm1 )470 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 454 471 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 455 472 END_3D 456 473 ! 457 474 CASE( 3 ) ! generic 458 DO_3D( 0, 0, 0, 0, 2, jpkm1 )475 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 459 476 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 460 477 END_3D … … 469 486 ! Warning : after this step, en : right hand side of the matrix 470 487 471 DO_3D( 0, 0, 0, 0, 2, jpkm1 )488 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 472 489 ! 473 490 ! psi / k … … 516 533 CASE ( 0 ) ! Dirichlet boundary conditions 517 534 ! 518 ! Surface value 519 zdep (:,:) = zhsro(:,:) * rl_sf ! Cosmetic 520 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 521 zd_lw(:,:,1) = psi(:,:,1) 522 zd_up(:,:,1) = 0._wp 523 zdiag(:,:,1) = 1._wp 524 ! 525 ! One level below 526 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 527 zdep (:,:) = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 528 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 529 zd_lw(:,:,2) = 0._wp 530 zd_up(:,:,2) = 0._wp 531 zdiag(:,:,2) = 1._wp 535 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 536 ! Surface value 537 zdep (ji,jj) = zhsro(ji,jj) * rl_sf ! Cosmetic 538 psi (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 539 zd_lw(ji,jj,1) = psi(ji,jj,1) 540 zd_up(ji,jj,1) = 0._wp 541 zdiag(ji,jj,1) = 1._wp 542 ! 543 ! One level below 544 zkar (ji,jj) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) ))) 545 zdep (ji,jj) = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj) 546 psi (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 547 zd_lw(ji,jj,2) = 0._wp 548 zd_up(ji,jj,2) = 0._wp 549 zdiag(ji,jj,2) = 1._wp 550 END_2D 532 551 ! 533 552 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 534 553 ! 535 ! Surface value: Dirichlet536 zdep (:,:) = zhsro(:,:) * rl_sf537 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)538 zd_lw(:,:,1) = psi(:,:,1)539 zd_up(:,:,1) = 0._wp540 zdiag(:,:,1) = 1._wp541 !542 ! Neumann condition at k=2543 DO_2D( 0, 0, 0, 0 ) !zdiag zd_lw not defined/used on the halo554 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 555 ! Surface value: Dirichlet 556 zdep (ji,jj) = zhsro(ji,jj) * rl_sf 557 psi (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1) 558 zd_lw(ji,jj,1) = psi(ji,jj,1) 559 zd_up(ji,jj,1) = 0._wp 560 zdiag(ji,jj,1) = 1._wp 561 ! 562 ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo 544 563 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) ! Remove zd_lw from zdiag 545 564 zd_lw(ji,jj,2) = 0._wp 565 ! 566 ! Set psi vertical flux at the surface: 567 zkar (ji,jj) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope 568 zdep (ji,jj) = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf) 569 zflxs(ji,jj) = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) & 570 & *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp) 571 zdep (ji,jj) = rsbc_psi1 * (zwall_psi(ji,jj,1)*p_avm(ji,jj,1)+zwall_psi(ji,jj,2)*p_avm(ji,jj,2)) * & 572 & ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.) 573 zflxs(ji,jj) = zdep(ji,jj) * zflxs(ji,jj) 574 psi (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm) 546 575 END_2D 547 !548 ! Set psi vertical flux at the surface:549 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope550 zdep (:,:) = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf)551 zflxs(:,:) = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) &552 & *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)553 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &554 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.)555 zflxs(:,:) = zdep(:,:) * zflxs(:,:)556 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)557 576 ! 558 577 END SELECT … … 569 588 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 570 589 ! ! Balance between the production and the dissipation terms 571 DO_2D( 0, 0, 0, 0)590 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 572 591 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 573 592 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 588 607 CASE ( 1 ) ! Neumman boundary condition 589 608 ! 590 DO_2D( 0, 0, 0, 0)609 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 591 610 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 592 611 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 … … 616 635 ! ---------------- 617 636 ! 618 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1637 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 619 638 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 620 639 END_3D 621 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1640 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 622 641 zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1) 623 642 END_3D 624 DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk643 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 ) ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 625 644 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 626 645 END_3D … … 632 651 ! 633 652 CASE( 0 ) ! k-kl (Mellor-Yamada) 634 DO_3D( 0, 0, 0, 0, 1, jpkm1 )653 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 635 654 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) 636 655 END_3D 637 656 ! 638 657 CASE( 1 ) ! k-eps 639 DO_3D( 0, 0, 0, 0, 1, jpkm1 )658 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 640 659 eps(ji,jj,jk) = psi(ji,jj,jk) 641 660 END_3D 642 661 ! 643 662 CASE( 2 ) ! k-w 644 DO_3D( 0, 0, 0, 0, 1, jpkm1 )663 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 645 664 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 646 665 END_3D … … 650 669 zex1 = ( 1.5_wp + rmm/rnn ) 651 670 zex2 = -1._wp / rnn 652 DO_3D( 0, 0, 0, 0, 1, jpkm1 )671 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 653 672 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 654 673 END_3D … … 658 677 ! Limit dissipation rate under stable stratification 659 678 ! -------------------------------------------------- 660 DO_3D ( 0, 0, 0, 0, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time679 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! Note that this set boundary conditions on hmxl_n at the same time 661 680 ! limitation 662 681 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 663 682 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 664 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 665 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 666 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 667 END_3D 683 END_3D 684 IF( ln_length_lim ) THEN ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 685 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 686 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 687 hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) 688 END_3D 689 ENDIF 668 690 669 691 ! … … 674 696 ! 675 697 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 676 DO_3D( 0, 0, 0, 0, 2, jpkm1 )698 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 677 699 ! zcof = l²/q² 678 700 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) … … 691 713 ! 692 714 CASE ( 2, 3 ) ! Canuto stability functions 693 DO_3D( 0, 0, 0, 0, 2, jpkm1 )715 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 694 716 ! zcof = l²/q² 695 717 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) … … 723 745 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 724 746 zstm(:,:,jpk) = 0. 725 DO_2D( 0, 0, 0, 0) ! update bottom with good values747 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! update bottom with good values 726 748 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 727 749 END_2D 728 750 729 zstt(:,:, 1) = wmask( :,:, 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)730 zstt(:,:,jpk) = wmask( :,:,jpk) ! default value not needed but avoid a bug when looking for undefined values (-fpe0)751 zstt(:,:, 1) = wmask(A2D(nn_hls), 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 752 zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) 731 753 732 754 !!gm should be done for ISF (top boundary cond.) … … 738 760 ! later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 739 761 ! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 740 DO_3D ( 0, 0, 0, 0, 1, jpk )762 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 741 763 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 742 764 zavt = zsqen * zstt(ji,jj,jk) … … 745 767 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 746 768 END_3D 747 p_avt( :,:,1) = 0._wp769 p_avt(A2D(nn_hls),1) = 0._wp 748 770 ! 749 771 IF(sn_cfctl%l_prtctl) THEN
Note: See TracChangeset
for help on using the changeset viewer.
