- Timestamp:
- 2020-03-11T10:26:18+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90
r11715 r12532 8 8 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 9 9 !! ! add velocity dependent coefficient and optional read in file 10 !! 4.0 ! 2020-02 (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum 10 11 !!---------------------------------------------------------------------- 11 12 … … 15 16 !!---------------------------------------------------------------------- 16 17 USE oce ! ocean dynamics and tracers 17 USE dom_oce ! ocean space and time domain 18 USE dom_oce ! ocean space and time domain 18 19 USE phycst ! physical constants 19 20 USE ldfslp ! lateral diffusion: slopes of mixing orientation … … 35 36 LOGICAL , PUBLIC :: ln_dynldf_OFF !: No operator (i.e. no explicit diffusion) 36 37 LOGICAL , PUBLIC :: ln_dynldf_lap !: laplacian operator 38 LOGICAL , PUBLIC :: ln_dynldf_bgm !: bi-Laplacian GM implementation via momentum 37 39 LOGICAL , PUBLIC :: ln_dynldf_blp !: bilaplacian operator 38 40 LOGICAL , PUBLIC :: ln_dynldf_lev !: iso-level direction … … 42 44 ! ! time invariant coefficients: aht = 1/2 Ud*Ld (lap case) 43 45 ! ! bht = 1/12 Ud*Ld^3 (blp case) 46 INTEGER , PUBLIC :: nn_bhm_ijk_t !: choice of time & space variations of the lateral bi-Laplacian GM coef. 44 47 REAL(wp), PUBLIC :: rn_Uv !: lateral viscous velocity [m/s] 45 48 REAL(wp), PUBLIC :: rn_Lv !: lateral viscous length [m] … … 50 53 ! ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T) 51 54 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 52 55 REAL(wp), PUBLIC :: rn_bhm_b !: lateral bi-Laplacian GM background eddy viscosity [m4/s] 56 REAL(wp), PUBLIC :: rn_bgmzcrit !: critical depth (>=0) for linear tapering of bi-Lap. GM coef. to zero at surface [m] 53 57 ! !!* Parameter to control the type of lateral viscous operator 54 58 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator … … 57 61 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 58 62 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator 59 !63 INTEGER, PARAMETER, PUBLIC :: np_bgm = 30 !: iso-level operator, bi-Laplacian GM 60 64 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 61 65 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 62 66 63 67 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 68 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bhmu, bhmv !: eddy viscosity coef. for bi-Laplacian GM at u- and v-points [m4/s] 64 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 65 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) … … 76 81 !!---------------------------------------------------------------------- 77 82 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 78 !! $Id $83 !! $Id: ldfdyn.F90 11653 2019-10-04 12:34:02Z davestorkey $ 79 84 !! Software governed by the CeCILL license (see ./LICENSE) 80 85 !!---------------------------------------------------------------------- … … 91 96 !! ln_dynldf_lap = T laplacian operator 92 97 !! ln_dynldf_blp = T bilaplacian operator 98 !! ln_dynldf_bgm = T bi-Laplacian GM operator 93 99 !! - the parameter nn_ahm_ijk_t: 94 100 !! nn_ahm_ijk_t = 0 => = constant … … 103 109 !! ( L^2|D| laplacian operator 104 110 !! or L^4|D|/8 bilaplacian operator ) 111 !! - the parameter nn_bhm_ijk_t: 112 !! nn_bhm_ijk_t = 11 => = F(z) : = constant, with BCs set to zero at top and bottom 113 !! nn_bhm_ijk_t = 12 => = F(z) : = constant in interior, linear taper to zero at surface, 114 !! for depths < rn_bgmzcrit, with BCs set to zero at top and bottom 115 !! nn_bhm_ijk_t = 13 => = F(z) : = bhm0 X dx^2 X dz^2, with BCs set to zero at top and bottom 105 116 !!---------------------------------------------------------------------- 106 117 INTEGER :: ji, jj, jk ! dummy loop indices 107 118 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 119 INTEGER :: iku, ikv, ikt ! local integer 108 120 REAL(wp) :: zah0, zah_max, zUfac ! local scalar 121 REAL(wp) :: bhm0, lhscrit, lhscritmax ! local scalar 109 122 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 110 123 !! 111 124 NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp, & ! type of operator 125 & ln_dynldf_bgm, & ! type of operator 112 126 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator 113 127 & nn_ahm_ijk_t , rn_Uv , rn_Lv, rn_ahm_b, & ! lateral eddy coefficient 128 & nn_bhm_ijk_t, rn_bhm_b, rn_bgmzcrit, & ! lateral bi-Lap. GM coefficient 114 129 & rn_csmc , rn_minfac , rn_maxfac ! Smagorinsky settings 115 130 !!---------------------------------------------------------------------- … … 134 149 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 135 150 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 151 WRITE(numout,*) ' bi-Laplacian GM operator ln_dynldf_bgm = ', ln_dynldf_bgm 152 WRITE(numout,*) ' critical depth for taper (if used) rn_bgmzcrit = ', rn_bgmzcrit 136 153 ! 137 154 WRITE(numout,*) ' direction of action :' … … 145 162 WRITE(numout,*) ' lateral viscous length (if cst) rn_Lv = ', rn_Lv, ' m' 146 163 WRITE(numout,*) ' background viscosity (iso-lap case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 147 !164 WRITE(numout,*) ' bi-Laplacian GM background viscosity rn_bhm_b = ', rn_bhm_b, ' m4/s' 148 165 WRITE(numout,*) ' Smagorinsky settings (nn_ahm_ijk_t = 32) :' 149 166 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc … … 159 176 nldf_dyn = np_ERROR 160 177 ioptio = 0 178 !CW exclude combination of lap+blp (as in existing code), but allow other combinations 161 179 IF( ln_dynldf_OFF ) THEN ; nldf_dyn = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 162 IF( ln_dynldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 163 IF( ln_dynldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 164 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 165 ! 180 IF( ln_dynldf_lap ) THEN ; ioptio = ioptio + 5 ; ENDIF 181 IF( ln_dynldf_blp ) THEN ; ioptio = ioptio + 5 ; ENDIF 182 IF( ln_dynldf_bgm ) THEN ; ioptio = ioptio + 1 ; ENDIF 183 ! 184 IF( (ioptio < 1).OR.(ioptio > 6) ) CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 4 operator options (NONE/lap/blp/bgm) but not (lap+blp)' ) 185 ! 186 166 187 IF(.NOT.ln_dynldf_OFF ) THEN !== direction ==>> type of operator ==! 167 188 ioptio = 0 … … 209 230 ENDIF 210 231 ! 232 IF( ln_dynldf_bgm ) THEN ! bi-Laplacian GM operator 233 IF( ln_zco ) THEN ! z-coordinate 234 IF( ln_dynldf_lev ) nldf_dyn = np_bgm ! iso-level = horizontal (no rotation) 235 IF( ln_dynldf_hor ) nldf_dyn = np_bgm ! iso-level = horizontal (no rotation) 236 IF( ln_dynldf_iso ) ierr = 3 ! iso-neutral ( rotation) 237 ENDIF 238 IF( ln_zps ) THEN ! z-coordinate with partial step 239 IF( ln_dynldf_lev ) nldf_dyn = np_bgm ! iso-level (no rotation) 240 IF( ln_dynldf_hor ) nldf_dyn = np_bgm ! iso-level (no rotation) 241 IF( ln_dynldf_iso ) ierr = 3 ! iso-neutral ( rotation) 242 ENDIF 243 IF( ln_sco ) THEN ! s-coordinate 244 IF( ln_dynldf_lev ) nldf_dyn = np_bgm ! iso-level (no rotation) 245 IF( ln_dynldf_hor ) ierr = 3 ! horizontal ( rotation) 246 IF( ln_dynldf_iso ) ierr = 3 ! iso-neutral ( rotation) 247 ENDIF 248 ENDIF 249 ! 250 211 251 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 212 252 ! 253 IF( ierr == 3 ) CALL ctl_stop( 'rotated bi-Laplacian GM operator does not exist' ) 254 ! 255 213 256 IF( nldf_dyn == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 214 257 ! … … 222 265 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 223 266 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 267 CASE( np_bgm ) ; WRITE(numout,*) ' ==>>> iso-level bi-Laplacian GM operator' 224 268 END SELECT 225 269 WRITE(numout,*) … … 233 277 ! 234 278 IF( ln_dynldf_OFF ) THEN 235 IF(lwp) WRITE(numout,*) ' ==>>> No viscous operator selected. ahmt and ahmf are not allocated' 279 !CW 280 !IF(lwp) WRITE(numout,*) ' ==>>> No viscous operator selected. ahmt and ahmf are not allocated' 281 IF(lwp) WRITE(numout,*) ' ==>>> No viscous operator selected. ahmt, ahmf, bhmu, bhmv are not allocated' 282 ! 236 283 RETURN 237 284 ! … … 239 286 ! 240 287 ! ! allocate the ahm arrays 241 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 288 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , bhmu(jpi,jpj,jpk) , bhmv(jpi,jpj,jpk) , STAT=ierr ) 289 ! 242 290 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 243 291 ! 244 292 ahmt(:,:,:) = 0._wp ! init to 0 needed 245 293 ahmf(:,:,:) = 0._wp 246 ! 247 ! ! value of lap/blp eddy mixing coef. 294 bhmu(:,:,:) = 0._wp ! init to 0 needed 295 bhmv(:,:,:) = 0._wp 296 ! ! value of lap/blp eddy mixing coef. 248 297 IF( ln_dynldf_lap ) THEN ; zUfac = r1_2 *rn_Uv ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 249 298 ELSEIF( ln_dynldf_blp ) THEN ; zUfac = r1_12*rn_Uv ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 299 ELSEIF( ln_dynldf_bgm ) THEN ; bhm0 = rn_bhm_b ; ; cl_Units = ' m4/s' ! bi-Laplacian GM 250 300 ENDIF 251 301 zah0 = zUfac * rn_Lv**inn ! mixing coefficient … … 265 315 ahmf(:,:,1) = zah0 266 316 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 267 ! 317 ! 268 318 CASE ( -20 ) !== fixed horizontal shape read in file ==! 269 319 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j) read in eddy_viscosity.nc file' … … 323 373 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 324 374 END SELECT 375 !CW Add separate structure options for bi-Laplacian GM, to allow it to be used in combination with other types of diffusion, above. 376 377 IF( ln_dynldf_bgm ) THEN 378 SELECT CASE (nn_bhm_ijk_t) 379 !CW Define bi-Laplacian GM diffusivity and test the related computational stability criterion 380 381 CASE( 11 ) !== fixed profile: constant except for zero at top and bottom, so that eddy-induced velocity, w*=0 ==! 382 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity = constant in interior, zero at top and bottom' 383 IF(lwp) WRITE(numout,*) ' interior viscous coef. = constant = ', bhm0, cl_Units 384 385 ! First set to constant in interior, all levels 386 DO jj = 2, jpjm1 387 DO ji = 2, jpim1 388 bhmu(ji,jj,:) = bhm0 389 bhmv(ji,jj,:) = bhm0 390 ENDDO 391 ENDDO 392 393 ! Surface BC : set to zero 394 bhmu(:,:,1) = 0._wp ; bhmv(:,:,1) = 0._wp 395 ! Flat bottom BC : set to zero 396 bhmu(:,:,jpk) = 0._wp ; bhmv(:,:,jpk) = 0._wp 397 ! Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps 398 ! and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp 399 DO jj = 2, jpjm1 400 DO ji = 2, jpim1 401 iku = mbku(ji,jj)+1 402 ikv = mbkv(ji,jj)+1 403 bhmu(:,:,iku) = 0._wp 404 bhmv(:,:,ikv) = 0._wp 405 ENDDO 406 ENDDO 407 408 !CW Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping 409 ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32 410 ! test at T-point - may need to revise this choice!!!! 411 412 ! Find domain maximum value of LHS, then test inequality. 413 414 ! initialise 415 lhscritmax=0._wp 416 417 DO jj = 2, jpjm1 418 DO ji = 2, jpim1 419 ikt = mbkt(ji,jj) !bottom last wet T-level 420 DO jk = 2, ikt 421 lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2) 422 IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit 423 ENDDO 424 ENDDO 425 ENDDO 426 427 IF ( lhscrit .ge. 1._wp/32._wp) THEN 428 IF(lwp) THEN 429 WRITE(numout,*) 430 WRITE(numout,*) ' WARNING : Bi-Laplacian GM eddy viscosity is not', & 431 & ' consistent with the model time step and grid setup: ', & 432 & ' LHS=',lhscrit, & 433 & 'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 434 WRITE(numout,*) ' =========' 435 WRITE(numout,*) 436 ENDIF 437 ! CW: warn, don't stop, as we may wish to violate this theoretical limit. 438 ! CALL ctl_stop( 'ldf_dyn_init', 'Inconsistent Bi-Lap. GM diffusivity') 439 ELSE 440 IF(lwp) THEN 441 WRITE(numout,*) 442 WRITE(numout,*) ' INFORMATION : Bi-Laplacian GM eddy viscosity is ', & 443 & ' consistent with the model time step and grid setup: ', & 444 & ' LHS=',lhscrit, & 445 & 'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 446 WRITE(numout,*) ' =========' 447 WRITE(numout,*) 448 ENDIF 449 ENDIF 450 !----------------------------- 451 CASE( 12 ) !== fixed profile: constant in interior, zero at bottom and linearly-tapering to zero at top, over depth shallower than rn_bgmzcrit ==! 452 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity = constant in interior, zero at bottom and' 453 IF(lwp) WRITE(numout,*) ' linearly-tapering to zero at top, over depth shallower than ', rn_bgmzcrit, ' metres.' 454 IF(lwp) WRITE(numout,*) ' Interior viscous coef. = constant = ', bhm0, cl_Units 455 456 ! Impose linear taper from interior value to zero at zero depth, for depths < rn_bgmzcrit 457 ! N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 458 459 DO jk = 1, jpk 460 461 IF (gdept_1d(jk) .lt. rn_bgmzcrit) THEN 462 463 DO jj = 2, jpjm1 464 DO ji = 2, jpim1 465 bhmu(ji,jj,jk) = bhm0 * gdept_1d(jk)/rn_bgmzcrit 466 bhmv(ji,jj,jk) = bhm0 * gdept_1d(jk)/rn_bgmzcrit 467 ENDDO 468 ENDDO 469 470 ELSE 471 472 ! constant at depth rn_bgmzcrit and larger 473 DO jj = 2, jpjm1 474 DO ji = 2, jpim1 475 bhmu(ji,jj,jk) = bhm0 476 bhmv(ji,jj,jk) = bhm0 477 ENDDO 478 ENDDO 479 480 ENDIF 481 482 ENDDO 483 484 485 ! Surface BC : set to zero 486 bhmu(:,:,1) = 0._wp ; bhmv(:,:,1) = 0._wp 487 ! Flat bottom BC : set to zero 488 bhmu(:,:,jpk) = 0._wp ; bhmv(:,:,jpk) = 0._wp 489 ! Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps 490 ! and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp 491 DO jj = 2, jpjm1 492 DO ji = 2, jpim1 493 iku = mbku(ji,jj)+1 494 ikv = mbkv(ji,jj)+1 495 bhmu(:,:,iku) = 0._wp 496 bhmv(:,:,ikv) = 0._wp 497 ENDDO 498 ENDDO 499 500 !-------------------------------- 501 CASE( 13 ) !== fixed profile: steady profile of the form bhm0 X dx^2 X dz^2 in interior, zero at bottom and top ==! 502 503 cl_Units = ' 1/s' 504 505 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity : steady profile of the form' 506 IF(lwp) WRITE(numout,*) ' bhm0 X dx^2 X dz^2 in interior, zero boundary conds. at bottom and top.' 507 IF(lwp) WRITE(numout,*) ' In this case, bhm0 is the constant of proportionality,' 508 IF(lwp) WRITE(numout,*) ' dimensioned as [diffusivity]/L^4, i.e. bhm0=', bhm0, cl_Units 509 510 cl_Units = ' m4/s' 511 512 ! N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 513 514 DO jk = 1, jpk 515 ! constant at depth rn_bgmzcrit and larger 516 DO jj = 2, jpjm1 517 DO ji = 2, jpim1 518 bhmu(ji,jj,jk) = bhm0*e1u(ji,jj)**2*e3t_n(ji,jj,jk)**2 519 bhmv(ji,jj,jk) = bhm0*e1v(ji,jj)**2*e3t_n(ji,jj,jk)**2 520 ENDDO 521 ENDDO 522 ENDDO 523 524 525 ! Surface BC : set to zero 526 bhmu(:,:,1) = 0._wp ; bhmv(:,:,1) = 0._wp 527 ! Flat bottom BC : set to zero 528 bhmu(:,:,jpk) = 0._wp ; bhmv(:,:,jpk) = 0._wp 529 ! Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps 530 ! and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp 531 DO jj = 2, jpjm1 532 DO ji = 2, jpim1 533 iku = mbku(ji,jj)+1 534 ikv = mbkv(ji,jj)+1 535 bhmu(:,:,iku) = 0._wp 536 bhmv(:,:,ikv) = 0._wp 537 ENDDO 538 ENDDO 539 540 !-------------------------------- 541 542 CASE DEFAULT 543 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_bhm_ijk_t, the type of space-time variation of bhm') 544 END SELECT 545 ENDIF ! ln_dynldf_bgm 325 546 ! 326 547 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation … … 331 552 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 332 553 ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1) 554 !CW 555 ELSEIF( ln_dynldf_bgm ) THEN !bi-Laplacian GM operator (mask only) 556 bhmu(:,:,1:jpkm1) = bhmu(:,:,1:jpkm1) * umask(:,:,1:jpkm1) 557 bhmv(:,:,1:jpkm1) = bhmv(:,:,1:jpkm1) * vmask(:,:,1:jpkm1) 558 ! 333 559 ENDIF 334 560 ENDIF … … 336 562 ENDIF 337 563 ! 338 END SUBROUTINE ldf_dyn_init339 340 341 SUBROUTINE ldf_dyn( kt )564 END SUBROUTINE ldf_dyn_init 565 566 567 SUBROUTINE ldf_dyn( kt ) 342 568 !!---------------------------------------------------------------------- 343 569 !! *** ROUTINE ldf_dyn *** … … 366 592 ! 367 593 SELECT CASE( nn_ahm_ijk_t ) !== Eddy vicosity coefficients ==! 368 !594 ! 369 595 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 370 596 ! … … 426 652 DO ji = 2, jpim1 427 653 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)654 & - ( 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 655 dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk) 430 656 END DO … … 434 660 DO ji = 1, jpim1 435 661 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)662 & + ( 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 663 dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk) 438 664 END DO … … 444 670 ! 445 671 DO jk = 1, jpkm1 446 !672 ! 447 673 DO jj = 2, jpjm1 ! T-point value 448 674 DO ji = fs_2, fs_jpim1 … … 453 679 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 454 680 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) ) )681 & r1_4 * ( dshesq(ji ,jj,jk) + dshesq(ji ,jj-1,jk) + & 682 & dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) 457 683 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 684 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) … … 469 695 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 470 696 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) ) )697 & r1_4 * ( dtensq(ji ,jj,jk) + dtensq(ji ,jj+1,jk) + & 698 & dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) 473 699 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 700 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt)
Note: See TracChangeset
for help on using the changeset viewer.
