- Timestamp:
- 2020-09-15T12:49:18+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/OCE/ZDF/zdfgls.F90
r13466 r13469 177 177 178 178 ! Compute surface, top and bottom friction at T-points 179 DO jj = 2, jpjm1 !== surface ocean friction 180 DO ji = fs_2, fs_jpim1 ! vector opt. 181 ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 182 END DO 183 END DO 179 DO_2D_00_00 180 ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 181 END_2D 184 182 ! 185 183 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 186 184 ! 187 185 IF( .NOT.ln_drg_OFF ) THEN !== top/bottom friction (explicit before friction) 188 DO jj = 2, jpjm1 ! bottom friction 189 DO ji = fs_2, fs_jpim1 ! vector opt. 190 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 191 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 192 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 193 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 194 END DO 195 END DO 186 DO_2D_00_00 187 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 188 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 189 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Nnn)+uu(ji-1,jj,mbkt(ji,jj),Nnn) ) )**2 & 190 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Nnn)+vv(ji,jj-1,mbkt(ji,jj),Nnn) ) )**2 ) 191 END_2D 196 192 IF( ln_isfcav ) THEN !top friction 197 DO jj = 2, jpjm1 198 DO ji = fs_2, fs_jpim1 ! vector opt. 199 zmsku = ( 2._wp - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 200 zmskv = ( 2._wp - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 201 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 & 202 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 203 END DO 204 END DO 193 DO_2D_00_00 194 zmsku = ( 2._wp - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 195 zmskv = ( 2._wp - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 196 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Nnn)+uu(ji-1,jj,mikt(ji,jj),Nnn) ) )**2 & 197 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Nnn)+vv(ji,jj-1,mikt(ji,jj),Nnn) ) )**2 ) 198 END_2D 205 199 ENDIF 206 200 ENDIF … … 224 218 zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1) + (1._wp - tmask(:,:,1))*rn_hsro 225 219 ! 226 DO jk = 2, jpkm1 !== Compute dissipation rate ==! 227 DO jj = 1, jpjm1 228 DO ji = 1, jpim1 229 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 230 END DO 231 END DO 232 END DO 220 DO_3D_10_10( 2, jpkm1 ) 221 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk) 222 END_3D 233 223 234 224 ! Save tke at before time step … … 237 227 238 228 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 239 DO jk = 2, jpkm1 240 DO jj = 2, jpjm1 241 DO ji = fs_2, fs_jpim1 ! vector opt. 242 zup = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 243 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 244 zcoef = ( zup / MAX( zdown, rsmall ) ) 245 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 229 DO_3D_00_00( 2, jpkm1 ) 230 zup = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 231 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 232 zcoef = ( zup / MAX( zdown, rsmall ) ) 233 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 234 END_3D 249 235 ENDIF 250 236 … … 262 248 ! Warning : after this step, en : right hand side of the matrix 263 249 264 DO jk = 2, jpkm1 265 DO jj = 2, jpjm1 266 DO ji = 2, jpim1 267 ! 268 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 269 ! 270 diss = eps(ji,jj,jk) ! dissipation 271 ! 272 zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 273 ! 274 zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk) ! production term 275 zdiss = zdir*(diss/en(ji,jj,jk)) +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 250 DO_3D_00_00( 2, jpkm1 ) 251 ! 252 buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk) ! stratif. destruction 253 ! 254 diss = eps(ji,jj,jk) ! dissipation 255 ! 256 zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 257 ! 258 zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk) ! production term 259 zdiss = zdir*(diss/en(ji,jj,jk)) +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term 276 260 !!gm better coding, identical results 277 261 ! zesh2 = p_sh2(ji,jj,jk) + zdir*buoy ! production term 278 262 ! zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term 279 263 !!gm 280 ! 281 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 282 ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) 283 ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 284 ! Otherwise, this should be rsc_psi/rsc_psi0 285 IF( ln_sigpsi ) THEN 286 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 287 zwall_psi(ji,jj,jk) = rsc_psi / & 288 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 289 ELSE 290 zwall_psi(ji,jj,jk) = 1._wp 291 ENDIF 292 ! 293 ! building the matrix 294 zcof = rfact_tke * tmask(ji,jj,jk) 295 ! ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 296 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 297 ! ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 298 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 299 ! ! diagonal 300 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 301 ! ! right hand side in en 302 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 303 END DO 304 END DO 305 END DO 264 ! 265 ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0 266 ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style) 267 ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries. 268 ! Otherwise, this should be rsc_psi/rsc_psi0 269 IF( ln_sigpsi ) THEN 270 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 271 zwall_psi(ji,jj,jk) = rsc_psi / & 272 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 273 ELSE 274 zwall_psi(ji,jj,jk) = 1._wp 275 ENDIF 276 ! 277 ! building the matrix 278 zcof = rfact_tke * tmask(ji,jj,jk) 279 ! ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 280 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 281 ! ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 282 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 283 ! ! diagonal 284 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 285 ! ! right hand side in en 286 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 287 END_3D 306 288 ! 307 289 zdiag(:,:,jpk) = 1._wp … … 360 342 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 361 343 ! ! Balance between the production and the dissipation terms 362 DO jj = 2, jpjm1 363 DO ji = fs_2, fs_jpim1 ! vector opt. 344 DO_2D_00_00 364 345 !!gm This means that bottom and ocean w-level above have a specified "en" value. Sure ???? 365 346 !! With thick deep ocean level thickness, this may be quite large, no ??? 366 347 !! in particular in ocean cavities where top stratification can be large... 367 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 368 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 348 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 349 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 350 ! 351 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 352 ! 353 ! Dirichlet condition applied at: 354 ! Bottom level (ibot) & Just above it (ibotm1) 355 zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp 356 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 357 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = 1._wp 358 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 359 END_2D 360 ! 361 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 362 DO_2D_00_00 363 itop = mikt(ji,jj) ! k top w-point 364 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 365 ! ! mask at the ocean surface points 366 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 369 367 ! 370 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 371 ! 368 !!gm TO BE VERIFIED !!! 372 369 ! Dirichlet condition applied at: 373 ! Bottom level (ibot) & Just above it (ibotm1) 374 zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp 375 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 376 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = 1._wp 377 en (ji,jj,ibot) = z_en ; en (ji,jj,ibotm1) = z_en 378 END DO 379 END DO 380 ! 381 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 382 DO jj = 2, jpjm1 383 DO ji = fs_2, fs_jpim1 ! vector opt. 384 itop = mikt(ji,jj) ! k top w-point 385 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 386 ! ! mask at the ocean surface points 387 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 388 ! 389 !!gm TO BE VERIFIED !!! 390 ! Dirichlet condition applied at: 391 ! top level (itop) & Just below it (itopp1) 392 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 393 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 394 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp 395 en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en 396 END DO 397 END DO 370 ! top level (itop) & Just below it (itopp1) 371 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 372 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 373 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = 1._wp 374 en (ji,jj,itop) = z_en ; en (ji,jj,itopp1) = z_en 375 END_2D 398 376 ENDIF 399 377 ! 400 378 CASE ( 1 ) ! Neumman boundary condition 401 379 ! 402 DO jj = 2, jpjm1 403 DO ji = fs_2, fs_jpim1 ! vector opt. 404 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 405 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 406 ! 407 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 380 DO_2D_00_00 381 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 382 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 383 ! 384 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 385 ! 386 ! Bottom level Dirichlet condition: 387 ! Bottom level (ibot) & Just above it (ibotm1) 388 ! Dirichlet ! Neumann 389 zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag 390 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 391 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 392 END_2D 393 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 394 DO_2D_00_00 395 itop = mikt(ji,jj) ! k top w-point 396 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 397 ! ! mask at the ocean surface points 398 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 408 399 ! 409 400 ! Bottom level Dirichlet condition: 410 401 ! Bottom level (ibot) & Just above it (ibotm1) 411 402 ! Dirichlet ! Neumann 412 zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag 413 zdiag(ji,jj,ibot) = 1._wp ; zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) 414 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp 415 END DO 416 END DO 417 IF( ln_isfcav) THEN ! top boundary (ocean cavity) 418 DO jj = 2, jpjm1 419 DO ji = fs_2, fs_jpim1 ! vector opt. 420 itop = mikt(ji,jj) ! k top w-point 421 itopp1 = mikt(ji,jj) + 1 ! k+1 1st w-point below the top one 422 ! ! mask at the ocean surface points 423 z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) ) 424 ! 425 ! Bottom level Dirichlet condition: 426 ! Bottom level (ibot) & Just above it (ibotm1) 427 ! Dirichlet ! Neumann 428 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag 429 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 430 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 431 END DO 432 END DO 403 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag 404 zdiag(ji,jj,itop) = 1._wp ; zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1) 405 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp 406 END_2D 433 407 ENDIF 434 408 ! … … 438 412 ! ---------------------------------------------------------- 439 413 ! 440 DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 441 DO jj = 2, jpjm1 442 DO ji = fs_2, fs_jpim1 ! vector opt. 443 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 444 END DO 445 END DO 446 END DO 447 DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 448 DO jj = 2, jpjm1 449 DO ji = fs_2, fs_jpim1 ! vector opt. 450 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) 451 END DO 452 END DO 453 END DO 454 DO jk = jpk-1, 2, -1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 455 DO jj = 2, jpjm1 456 DO ji = fs_2, fs_jpim1 ! vector opt. 457 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 458 END DO 459 END DO 460 END DO 414 DO_3D_00_00( 2, jpkm1 ) 415 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 416 END_3D 417 DO_3D_00_00( 2, jpk ) 418 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) 419 END_3D 420 DO_3D_00_00( jpk-1, 2, -1 ) 421 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 422 END_3D 461 423 ! ! set the minimum value of tke 462 424 en(:,:,:) = MAX( en(:,:,:), rn_emin ) … … 471 433 ! 472 434 CASE( 0 ) ! k-kl (Mellor-Yamada) 473 DO jk = 2, jpkm1 474 DO jj = 2, jpjm1 475 DO ji = fs_2, fs_jpim1 ! vector opt. 476 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 477 END DO 478 END DO 479 END DO 435 DO_3D_00_00( 2, jpkm1 ) 436 psi(ji,jj,jk) = eb(ji,jj,jk) * hmxl_b(ji,jj,jk) 437 END_3D 480 438 ! 481 439 CASE( 1 ) ! k-eps 482 DO jk = 2, jpkm1 483 DO jj = 2, jpjm1 484 DO ji = fs_2, fs_jpim1 ! vector opt. 485 psi(ji,jj,jk) = eps(ji,jj,jk) 486 END DO 487 END DO 488 END DO 440 DO_3D_00_00( 2, jpkm1 ) 441 psi(ji,jj,jk) = eps(ji,jj,jk) 442 END_3D 489 443 ! 490 444 CASE( 2 ) ! k-w 491 DO jk = 2, jpkm1 492 DO jj = 2, jpjm1 493 DO ji = fs_2, fs_jpim1 ! vector opt. 494 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 495 END DO 496 END DO 497 END DO 445 DO_3D_00_00( 2, jpkm1 ) 446 psi(ji,jj,jk) = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) ) 447 END_3D 498 448 ! 499 449 CASE( 3 ) ! generic 500 DO jk = 2, jpkm1 501 DO jj = 2, jpjm1 502 DO ji = fs_2, fs_jpim1 ! vector opt. 503 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 504 END DO 505 END DO 506 END DO 450 DO_3D_00_00( 2, jpkm1 ) 451 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 452 END_3D 507 453 ! 508 454 END SELECT … … 515 461 ! Warning : after this step, en : right hand side of the matrix 516 462 517 DO jk = 2, jpkm1 518 DO jj = 2, jpjm1 519 DO ji = fs_2, fs_jpim1 ! vector opt. 520 ! 521 ! psi / k 522 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 523 ! 524 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 525 zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 526 ! 527 rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 528 ! 529 ! shear prod. - stratif. destruction 530 prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 531 ! 532 ! stratif. destruction 533 buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 534 ! 535 ! shear prod. - stratif. destruction 536 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 537 ! 538 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 539 ! 540 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 541 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 542 ! 543 ! building the matrix 544 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 545 ! ! lower diagonal 546 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 547 ! ! upper diagonal 548 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 549 ! ! diagonal 550 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 551 ! ! right hand side in psi 552 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 553 END DO 554 END DO 555 END DO 463 DO_3D_00_00( 2, jpkm1 ) 464 ! 465 ! psi / k 466 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 467 ! 468 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) 469 zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) ) 470 ! 471 rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p 472 ! 473 ! shear prod. - stratif. destruction 474 prod = rpsi1 * zratio * p_sh2(ji,jj,jk) 475 ! 476 ! stratif. destruction 477 buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) ) 478 ! 479 ! shear prod. - stratif. destruction 480 diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk) 481 ! 482 zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy ) ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0) 483 ! 484 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 485 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 486 ! 487 ! building the matrix 488 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 489 ! ! lower diagonal 490 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) ) 491 ! ! upper diagonal 492 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 493 ! ! diagonal 494 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 495 ! ! right hand side in psi 496 psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 497 END_3D 556 498 ! 557 499 zdiag(:,:,jpk) = 1._wp … … 615 557 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 616 558 ! ! Balance between the production and the dissipation terms 617 DO jj = 2, jpjm1 618 DO ji = fs_2, fs_jpim1 ! vector opt. 619 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 620 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 621 zdep(ji,jj) = vkarmn * r_z0_bot 622 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 623 zd_lw(ji,jj,ibot) = 0._wp 624 zd_up(ji,jj,ibot) = 0._wp 625 zdiag(ji,jj,ibot) = 1._wp 626 ! 627 ! Just above last level, Dirichlet condition again (GOTM like) 628 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) ) 629 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 630 zd_lw(ji,jj,ibotm1) = 0._wp 631 zd_up(ji,jj,ibotm1) = 0._wp 632 zdiag(ji,jj,ibotm1) = 1._wp 633 END DO 634 END DO 559 DO_2D_00_00 560 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 561 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 562 zdep(ji,jj) = vkarmn * r_z0_bot 563 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 564 zd_lw(ji,jj,ibot) = 0._wp 565 zd_up(ji,jj,ibot) = 0._wp 566 zdiag(ji,jj,ibot) = 1._wp 567 ! 568 ! Just above last level, Dirichlet condition again (GOTM like) 569 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) ) 570 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 571 zd_lw(ji,jj,ibotm1) = 0._wp 572 zd_up(ji,jj,ibotm1) = 0._wp 573 zdiag(ji,jj,ibotm1) = 1._wp 574 END_2D 635 575 ! 636 576 CASE ( 1 ) ! Neumman boundary condition 637 577 ! 638 DO jj = 2, jpjm1 639 DO ji = fs_2, fs_jpim1 ! vector opt. 640 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 641 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 642 ! 643 ! Bottom level Dirichlet condition: 644 zdep(ji,jj) = vkarmn * r_z0_bot 645 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 646 ! 647 zd_lw(ji,jj,ibot) = 0._wp 648 zd_up(ji,jj,ibot) = 0._wp 649 zdiag(ji,jj,ibot) = 1._wp 650 ! 651 ! Just above last level: Neumann condition with flux injection 652 zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 653 zd_up(ji,jj,ibotm1) = 0. 654 ! 655 ! Set psi vertical flux at the bottom: 656 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1) 657 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 658 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 659 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 660 END DO 661 END DO 578 DO_2D_00_00 579 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 580 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 581 ! 582 ! Bottom level Dirichlet condition: 583 zdep(ji,jj) = vkarmn * r_z0_bot 584 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 585 ! 586 zd_lw(ji,jj,ibot) = 0._wp 587 zd_up(ji,jj,ibot) = 0._wp 588 zdiag(ji,jj,ibot) = 1._wp 589 ! 590 ! Just above last level: Neumann condition with flux injection 591 zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag 592 zd_up(ji,jj,ibotm1) = 0. 593 ! 594 ! Set psi vertical flux at the bottom: 595 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1) 596 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 597 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 598 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 599 END_2D 662 600 ! 663 601 END SELECT … … 666 604 ! ---------------- 667 605 ! 668 DO jk = 2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 669 DO jj = 2, jpjm1 670 DO ji = fs_2, fs_jpim1 ! vector opt. 671 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 672 END DO 673 END DO 674 END DO 675 DO jk = 2, jpk ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 676 DO jj = 2, jpjm1 677 DO ji = fs_2, fs_jpim1 ! vector opt. 678 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) 679 END DO 680 END DO 681 END DO 682 DO jk = jpk-1, 2, -1 ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 683 DO jj = 2, jpjm1 684 DO ji = fs_2, fs_jpim1 ! vector opt. 685 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 686 END DO 687 END DO 688 END DO 606 DO_3D_00_00( 2, jpkm1 ) 607 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 608 END_3D 609 DO_3D_00_00( 2, jpk ) 610 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) 611 END_3D 612 DO_3D_00_00( jpk-1, 2, -1 ) 613 psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 614 END_3D 689 615 690 616 ! Set dissipation … … 694 620 ! 695 621 CASE( 0 ) ! k-kl (Mellor-Yamada) 696 DO jk = 1, jpkm1 697 DO jj = 2, jpjm1 698 DO ji = fs_2, fs_jpim1 ! vector opt. 699 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) 700 END DO 701 END DO 702 END DO 622 DO_3D_00_00( 1, jpkm1 ) 623 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) 624 END_3D 703 625 ! 704 626 CASE( 1 ) ! k-eps 705 DO jk = 1, jpkm1 706 DO jj = 2, jpjm1 707 DO ji = fs_2, fs_jpim1 ! vector opt. 708 eps(ji,jj,jk) = psi(ji,jj,jk) 709 END DO 710 END DO 711 END DO 627 DO_3D_00_00( 1, jpkm1 ) 628 eps(ji,jj,jk) = psi(ji,jj,jk) 629 END_3D 712 630 ! 713 631 CASE( 2 ) ! k-w 714 DO jk = 1, jpkm1 715 DO jj = 2, jpjm1 716 DO ji = fs_2, fs_jpim1 ! vector opt. 717 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 718 END DO 719 END DO 720 END DO 632 DO_3D_00_00( 1, jpkm1 ) 633 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 634 END_3D 721 635 ! 722 636 CASE( 3 ) ! generic … … 724 638 zex1 = ( 1.5_wp + rmm/rnn ) 725 639 zex2 = -1._wp / rnn 726 DO jk = 1, jpkm1 727 DO jj = 2, jpjm1 728 DO ji = fs_2, fs_jpim1 ! vector opt. 729 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 730 END DO 731 END DO 732 END DO 640 DO_3D_00_00( 1, jpkm1 ) 641 eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2 642 END_3D 733 643 ! 734 644 END SELECT … … 736 646 ! Limit dissipation rate under stable stratification 737 647 ! -------------------------------------------------- 738 DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time 739 DO jj = 2, jpjm1 740 DO ji = fs_2, fs_jpim1 ! vector opt. 741 ! limitation 742 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 743 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 744 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 745 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 746 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) ) 747 END DO 748 END DO 749 END DO 648 DO_3D_00_00( 1, jpkm1 ) 649 ! limitation 650 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 651 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 652 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 653 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 654 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) ) 655 END_3D 750 656 751 657 ! … … 756 662 ! 757 663 CASE ( 0 , 1 ) ! Galperin or Kantha-Clayson stability functions 758 DO jk = 2, jpkm1 759 DO jj = 2, jpjm1 760 DO ji = fs_2, fs_jpim1 ! vector opt. 761 ! zcof = l²/q² 762 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 763 ! Gh = -N²l²/q² 764 gh = - rn2(ji,jj,jk) * zcof 765 gh = MIN( gh, rgh0 ) 766 gh = MAX( gh, rghmin ) 767 ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 768 sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 769 sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 770 ! 771 ! Store stability function in zstt and zstm 772 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 773 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 774 END DO 775 END DO 776 END DO 664 DO_3D_00_00( 2, jpkm1 ) 665 ! zcof = l²/q² 666 zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) ) 667 ! Gh = -N²l²/q² 668 gh = - rn2(ji,jj,jk) * zcof 669 gh = MIN( gh, rgh0 ) 670 gh = MAX( gh, rghmin ) 671 ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin) 672 sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) ) 673 sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh) 674 ! 675 ! Store stability function in zstt and zstm 676 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 677 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 678 END_3D 777 679 ! 778 680 CASE ( 2, 3 ) ! Canuto stability functions 779 DO jk = 2, jpkm1 780 DO jj = 2, jpjm1 781 DO ji = fs_2, fs_jpim1 ! vector opt. 782 ! zcof = l²/q² 783 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 784 ! Gh = -N²l²/q² 785 gh = - rn2(ji,jj,jk) * zcof 786 gh = MIN( gh, rgh0 ) 787 gh = MAX( gh, rghmin ) 788 gh = gh * rf6 789 ! Gm = M²l²/q² Shear number 790 shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 791 gm = MAX( shr * zcof , 1.e-10 ) 792 gm = gm * rf6 793 gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 794 ! Stability functions from Canuto 795 rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 796 sm = (rs0 - rs1*gh + rs2*gm) / rcff 797 sh = (rs4 - rs5*gh + rs6*gm) / rcff 798 ! 799 ! Store stability function in zstt and zstm 800 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 801 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 802 END DO 803 END DO 804 END DO 681 DO_3D_00_00( 2, jpkm1 ) 682 ! zcof = l²/q² 683 zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) ) 684 ! Gh = -N²l²/q² 685 gh = - rn2(ji,jj,jk) * zcof 686 gh = MIN( gh, rgh0 ) 687 gh = MAX( gh, rghmin ) 688 gh = gh * rf6 689 ! Gm = M²l²/q² Shear number 690 shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall ) 691 gm = MAX( shr * zcof , 1.e-10 ) 692 gm = gm * rf6 693 gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm ) 694 ! Stability functions from Canuto 695 rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm 696 sm = (rs0 - rs1*gh + rs2*gm) / rcff 697 sh = (rs4 - rs5*gh + rs6*gm) / rcff 698 ! 699 ! Store stability function in zstt and zstm 700 zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk) 701 zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk) 702 END_3D 805 703 ! 806 704 END SELECT … … 813 711 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 814 712 zstm(:,:,jpk) = 0. 815 DO jj = 2, jpjm1 ! update bottom with good values 816 DO ji = fs_2, fs_jpim1 ! vector opt. 817 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 818 END DO 819 END DO 713 DO_2D_00_00 714 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) 715 END_2D 820 716 821 717 zstt(:,:, 1) = wmask(:,:, 1) ! default value not needed but avoid a bug when looking for undefined values (-fpe0) … … 830 726 ! later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk) 831 727 ! for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0) 832 DO jk = 1, jpk 833 DO jj = 2, jpjm1 834 DO ji = fs_2, fs_jpim1 ! vector opt. 835 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 836 zavt = zsqen * zstt(ji,jj,jk) 837 zavm = zsqen * zstm(ji,jj,jk) 838 p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 839 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 840 END DO 841 END DO 842 END DO 728 DO_3D_00_00( 1, jpk ) 729 zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk) 730 zavt = zsqen * zstt(ji,jj,jk) 731 zavm = zsqen * zstm(ji,jj,jk) 732 p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine 733 p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom 734 END_3D 843 735 p_avt(:,:,1) = 0._wp 844 736 !
Note: See TracChangeset
for help on using the changeset viewer.
