Changeset 10783
- Timestamp:
- 2019-03-20T19:36:26+01:00 (6 years ago)
- File:
-
- 1 edited
-
NEMO/releases/release-4.0/src/OCE/LDF/ldfdyn.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/releases/release-4.0/src/OCE/LDF/ldfdyn.F90
r10425 r10783 62 62 63 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ):: dtensq !: horizontal tension squared (Smagorinsky only)65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ):: dshesq !: horizontal shearing strain squared (Smagorinsky only)64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) 66 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2 67 67 … … 242 242 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 243 243 ! 244 ahmt(:,:, jpk) = 0._wp ! last level always 0245 ahmf(:,:, jpk) = 0._wp244 ahmt(:,:,:) = 0._wp ! init to 0 needed 245 ahmf(:,:,:) = 0._wp 246 246 ! 247 247 ! ! value of lap/blp eddy mixing coef. … … 310 310 ! 311 311 ! ! allocate arrays used in ldf_dyn. 312 ALLOCATE( dtensq(jpi,jpj ) , dshesq(jpi,jpj) , esqt(jpi,jpj) ,esqf(jpi,jpj) , STAT=ierr )312 ALLOCATE( dtensq(jpi,jpj,jpk) , dshesq(jpi,jpj,jpk) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 313 313 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 314 314 ! 315 DO jj = 2, jpjm1! Set local gridscale values316 DO ji = fs_2, fs_jpim1317 esqt(ji,jj) = ( e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2318 esqf(ji,jj) = ( e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2315 DO jj = 1, jpj ! Set local gridscale values 316 DO ji = 1, jpi 317 esqt(ji,jj) = ( e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 318 esqf(ji,jj) = ( e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 319 319 END DO 320 320 END DO … … 359 359 ! 360 360 INTEGER :: ji, jj, jk ! dummy loop indices 361 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, ze tmax, zefmax ! local scalar362 REAL(wp) :: zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb ! local scalar361 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zemax ! local scalar (option 31) 362 REAL(wp) :: zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb ! local scalar (option 32) 363 363 !!---------------------------------------------------------------------- 364 364 ! … … 373 373 DO jj = 2, jpjm1 374 374 DO ji = fs_2, fs_jpim1 375 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 376 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 377 zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 378 ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk) ! 288= 12*12 * 2 379 END DO 380 END DO 381 DO jj = 1, jpjm1 382 DO ji = 1, fs_jpim1 375 383 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 376 384 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 377 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 378 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 379 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 380 ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk) ! 288= 12*12 * 2 381 ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk) 385 zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 386 ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk) ! 288= 12*12 * 2 382 387 END DO 383 388 END DO … … 387 392 DO jj = 2, jpjm1 388 393 DO ji = fs_2, fs_jpim1 394 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 395 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 396 zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 397 ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax ) * zemax * tmask(ji,jj,jk) 398 END DO 399 END DO 400 DO jj = 1, jpjm1 401 DO ji = 1, fs_jpim1 389 402 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 390 403 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 391 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 392 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 393 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 394 ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax ) * zetmax * tmask(ji,jj,jk) 395 ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax ) * zefmax * fmask(ji,jj,jk) 404 zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 405 ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax ) * zemax * fmask(ji,jj,jk) 396 406 END DO 397 407 END DO … … 406 416 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! laplacian operator : (C_smag/pi)^2 L^2 |D| 407 417 ! 408 zcmsmag = (rn_csmc/rpi)**2! (C_smag/pi)^2409 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )! lower limit stability factor scaling410 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )! upper limit stability factor scaling418 zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2 419 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling 420 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt ) ! upper limit stability factor scaling 411 421 IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide |U|L^3/12 lower limit instead 412 422 ! ! of |U|L^3/16 in blp case 413 423 DO jk = 1, jpkm1 414 424 ! 415 DO jj = 2, jpj 416 DO ji = 2, jpi 417 zdb = ( ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & 418 & * r1_e1t(ji,jj) * e2t(ji,jj) & 419 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 420 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 421 dtensq(ji,jj) = zdb * zdb 425 DO jj = 2, jpjm1 426 DO ji = 2, jpim1 427 zdb = ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) * r1_e1t(ji,jj) * e2t(ji,jj) & 428 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj) 429 dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk) 422 430 END DO 423 431 END DO … … 425 433 DO jj = 1, jpjm1 426 434 DO ji = 1, jpim1 427 zdb = ( ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) & 428 & * r1_e2f(ji,jj) * e1f(ji,jj) & 429 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 430 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 431 dshesq(ji,jj) = zdb * zdb 435 zdb = ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) * r1_e2f(ji,jj) * e1f(ji,jj) & 436 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj) 437 dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk) 432 438 END DO 433 439 END DO 434 440 ! 435 DO jj = 2, jpjm1 441 END DO 442 ! 443 CALL lbc_lnk_multi( 'ldfdyn', dtensq, 'T', 1. ) ! lbc_lnk on dshesq not needed 444 ! 445 DO jk = 1, jpkm1 446 ! 447 DO jj = 2, jpjm1 ! T-point value 436 448 DO ji = fs_2, fs_jpim1 449 ! 450 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 451 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 452 ! 453 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 454 ahmt(ji,jj,jk) = zdelta * SQRT( dtensq(ji ,jj,jk) + & 455 & r1_4 * ( dshesq(ji ,jj,jk) + dshesq(ji ,jj-1,jk) + & 456 & dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) 457 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 458 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 459 ! 460 END DO 461 END DO 462 ! 463 DO jj = 1, jpjm1 ! F-point value 464 DO ji = 1, fs_jpim1 437 465 ! 438 466 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 439 467 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 440 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 441 ! T-point value 442 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 443 ahmt(ji,jj,jk) = zdelta * sqrt( dtensq(ji,jj) + & 444 & r1_4 * ( dshesq(ji,jj) + dshesq(ji,jj-1) + & 445 & dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 446 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), & 447 & SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 448 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 449 ! F-point value 468 ! 450 469 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 451 ahmf(ji,jj,jk) = zdelta * sqrt( dshesq(ji,jj) + & 452 & r1_4 * ( dtensq(ji,jj) + dtensq(ji,jj+1) + & 453 & dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 454 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), & 455 & SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 456 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 470 ahmf(ji,jj,jk) = zdelta * SQRT( dshesq(ji ,jj,jk) + & 471 & r1_4 * ( dtensq(ji ,jj,jk) + dtensq(ji ,jj+1,jk) + & 472 & dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) 473 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 474 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 457 475 ! 458 476 END DO 459 477 END DO 478 ! 460 479 END DO 461 480 ! … … 470 489 DO ji = fs_2, fs_jpim1 471 490 ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 491 END DO 492 END DO 493 DO jj = 1, jpjm1 494 DO ji = 1, fs_jpim1 472 495 ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 473 496 END DO
Note: See TracChangeset
for help on using the changeset viewer.
