Changeset 12776
- Timestamp:
- 2020-04-20T11:24:26+02:00 (5 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE
- Files:
-
- 2 edited
-
DYN/dynldf.F90 (modified) (2 diffs)
-
LDF/ldfdyn.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90
r12535 r12776 67 67 CASE ( np_lap ) ; CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 ) ! iso-level laplacian 68 68 CASE ( np_lap_i ) ; CALL dyn_ldf_iso( kt ) ! rotated laplacian 69 CASE ( np_bgm ) ; CALL dyn_ldf_bgm( kt, ub, vb, ua, va) ! bi-Laplacian GM parameterisation70 69 CASE ( np_blp ) ; CALL dyn_ldf_blp( kt, ub, vb, ua, va ) ! iso-level bi-laplacian 71 70 ! 72 71 END SELECT 72 73 IF( ln_dynldf_bgm ) CALL dyn_ldf_bgm( kt, ub, vb, ua, va) ! bi-Laplacian GM parameterisation 73 74 74 75 IF( l_trddyn ) THEN ! save the horizontal diffusive trends for further diagnostics … … 106 107 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> iso-level laplacian operator' 107 108 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 108 CASE( np_bgm ) ; WRITE(numout,*) ' ==>>> bi-Laplacian GM parameterisation via momentum equation'109 109 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 110 110 END SELECT 111 IF( ln_dynldf_bgm ) WRITE(numout,*) ' ==>>> biharmonic Gent-McWilliams operator turned on' 111 112 ENDIF 112 113 ! -
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90
r12775 r12776 61 61 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 62 62 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator 63 INTEGER, PARAMETER, PUBLIC :: np_bgm = 30 !: iso-level operator, bi-Laplacian GM64 63 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 65 64 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. … … 180 179 IF( ln_dynldf_lap ) THEN ; ioptio = ioptio + 5 ; ENDIF 181 180 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)' ) 181 ! 182 IF( (ioptio < 1).OR.(ioptio > 6) ) CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 4 operator options (NONE/lap/blp) but not (lap+blp)' ) 185 183 ! 186 184 … … 230 228 ENDIF 231 229 ! 232 IF( ln_dynldf_bgm ) THEN ! bi-Laplacian GM operator233 IF( ln_zco ) THEN ! z-coordinate234 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 ENDIF238 IF( ln_zps ) THEN ! z-coordinate with partial step239 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 ENDIF243 IF( ln_sco ) THEN ! s-coordinate244 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 ENDIF248 ENDIF249 !250 230 251 231 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 252 !253 IF( ierr == 3 ) CALL ctl_stop( 'rotated bi-Laplacian GM operator does not exist' )254 232 ! 255 233 … … 265 243 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 266 244 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 267 CASE( np_bgm ) ; WRITE(numout,*) ' ==>>> iso-level bi-Laplacian GM operator'268 245 END SELECT 269 246 WRITE(numout,*) … … 277 254 ! 278 255 IF( ln_dynldf_OFF ) THEN 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, bhm are not allocated' 282 ! 283 RETURN 256 ! 257 IF(lwp) WRITE(numout,*) ' ==>>> No viscous operator selected. ahmt, ahmf are not allocated' 284 258 ! 285 259 ELSE !== a lateral diffusion operator is used ==! … … 292 266 ahmt(:,:,:) = 0._wp ! init to 0 needed 293 267 ahmf(:,:,:) = 0._wp 294 bhm(:,:,:) = 0._wp ! init to 0 needed295 268 ! ! value of lap/blp eddy mixing coef. 296 269 IF( ln_dynldf_lap ) THEN ; zUfac = r1_2 *rn_Uv ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 297 270 ELSEIF( ln_dynldf_blp ) THEN ; zUfac = r1_12*rn_Uv ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 298 ELSEIF( ln_dynldf_bgm ) THEN ; bhm0 = rn_bhm_b ; ; cl_Units = ' m4/s' ! bi-Laplacian GM299 271 ENDIF 300 272 zah0 = zUfac * rn_Lv**inn ! mixing coefficient … … 372 344 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 373 345 END SELECT 374 !CW Add separate structure options for bi-Laplacian GM, to allow it to be used in combination with other types of diffusion, above. 375 376 IF( ln_dynldf_bgm ) THEN 377 SELECT CASE (nn_bhm_ijk_t) 378 !CW Define bi-Laplacian GM diffusivity and test the related computational stability criterion 379 380 CASE( 11 ) !== fixed profile: constant except for zero at top and bottom, so that eddy-induced velocity, w*=0 ==! 381 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity = constant in interior, zero at top and bottom' 382 IF(lwp) WRITE(numout,*) ' interior viscous coef. = constant = ', bhm0, cl_Units 383 384 ! First set to constant in interior, all levels 385 DO jj = 2, jpjm1 386 DO ji = 2, jpim1 387 bhm(ji,jj,:) = bhm0 388 ENDDO 389 ENDDO 390 391 ! Surface BC : set to zero 392 bhm(:,:,1) = 0._wp 393 ! Flat bottom BC : set to zero 394 bhm(:,:,jpk) = 0._wp 395 ! Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm 396 397 !CW Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping 398 ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32 399 ! test at T-point - may need to revise this choice!!!! 400 401 ! Find domain maximum value of LHS, then test inequality. 402 403 ! initialise 404 lhscritmax=0._wp 405 406 DO jj = 2, jpjm1 407 DO ji = 2, jpim1 408 ikw = mbkt(ji,jj) !bottom last wet T-level 409 DO jk = 2, ikw 410 lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2) 411 IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit 412 ENDDO 413 ENDDO 414 ENDDO 415 416 IF ( lhscritmax .ge. 1._wp/32._wp) THEN 417 IF(lwp) THEN 418 WRITE(numout,*) 419 WRITE(numout,*) ' WARNING : Bi-Laplacian GM eddy viscosity is not', & 420 & ' consistent with the model time step and grid setup: ', & 421 & ' LHS=',lhscrit, & 422 & 'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 423 WRITE(numout,*) ' =========' 424 WRITE(numout,*) 425 ENDIF 426 ! CW: warn, don't stop, as we may wish to violate this theoretical limit. 427 ! CALL ctl_stop( 'ldf_dyn_init', 'Inconsistent Bi-Lap. GM diffusivity') 428 ELSE 429 IF(lwp) THEN 430 WRITE(numout,*) 431 WRITE(numout,*) ' INFORMATION : Bi-Laplacian GM eddy viscosity is ', & 432 & ' consistent with the model time step and grid setup: ', & 433 & ' LHS=',lhscrit, & 434 & 'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 435 WRITE(numout,*) ' =========' 436 WRITE(numout,*) 437 ENDIF 438 ENDIF 439 !----------------------------- 440 CASE( 12 ) !== fixed profile: constant in interior, zero at bottom and linearly-tapering to zero at top, over depth shallower than rn_bgmzcrit ==! 441 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity = constant in interior, zero at bottom and' 442 IF(lwp) WRITE(numout,*) ' linearly-tapering to zero at top, over depth shallower than ', rn_bgmzcrit, ' metres.' 443 IF(lwp) WRITE(numout,*) ' Interior viscous coef. = constant = ', bhm0, cl_Units 444 445 ! Impose linear taper from interior value to zero at zero depth, for depths < rn_bgmzcrit 446 ! N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 447 448 DO jk = 1, jpk 449 450 IF (gdept_1d(jk) .lt. rn_bgmzcrit) THEN 451 452 DO jj = 2, jpjm1 453 DO ji = 2, jpim1 454 bhm(ji,jj,jk) = bhm0 * gdepw_1d(jk)/rn_bgmzcrit 455 ENDDO 456 ENDDO 457 458 ELSE 459 460 ! constant at depth rn_bgmzcrit and larger 461 DO jj = 2, jpjm1 462 DO ji = 2, jpim1 463 bhm(ji,jj,jk) = bhm0 464 ENDDO 465 ENDDO 466 467 ENDIF 468 469 ENDDO 470 471 ! Surface BC : set to zero 472 bhm(:,:,1) = 0._wp 473 ! Flat bottom BC : set to zero 474 bhm(:,:,jpk) = 0._wp 475 ! Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm 476 477 !-------------------------------- 478 CASE( 13 ) !== fixed profile: steady profile of the form bhm0 X dx^2 X dz^2 in interior, zero at bottom and top ==! 479 480 cl_Units = ' 1/s' 481 482 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity : steady profile of the form' 483 IF(lwp) WRITE(numout,*) ' bhm0 X dx^2 X dz^2 in interior, zero boundary conds. at bottom and top.' 484 IF(lwp) WRITE(numout,*) ' In this case, bhm0 is the constant of proportionality,' 485 IF(lwp) WRITE(numout,*) ' dimensioned as [diffusivity]/L^4, i.e. bhm0=', bhm0, cl_Units 486 487 cl_Units = ' m4/s' 488 489 ! N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 490 491 DO jk = 1, jpk 492 ! constant at depth rn_bgmzcrit and larger 493 DO jj = 2, jpjm1 494 DO ji = 2, jpim1 495 bhm(ji,jj,jk) = bhm0*e1t(ji,jj)**2*e3w_n(ji,jj,jk)**2 496 ENDDO 497 ENDDO 498 ENDDO 499 500 ! Surface BC : set to zero 501 bhm(:,:,1) = 0._wp 502 ! Flat bottom BC : set to zero 503 bhm(:,:,jpk) = 0._wp 504 ! Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm 505 506 !-------------------------------- 507 508 CASE DEFAULT 509 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_bhm_ijk_t, the type of space-time variation of bhm') 510 END SELECT 511 ENDIF ! ln_dynldf_bgm 346 512 347 ! 513 348 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation … … 518 353 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 519 354 ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1) 520 !CW521 ELSEIF( ln_dynldf_bgm ) THEN !bi-Laplacian GM operator (mask only)522 bhm(:,:,1:jpkm1) = bhm(:,:,1:jpkm1) * wmask(:,:,1:jpkm1)523 CALL lbc_lnk( 'ldf_dyn_init', bhm , 'W', 1. )524 !525 355 ENDIF 526 356 ENDIF 527 357 ! 528 358 ENDIF 359 ! 360 !CW Add separate structure options for bi-Laplacian GM, to allow it to be used in combination with other types of diffusion, above. 361 IF( ln_dynldf_bgm ) THEN 362 IF(lwp) WRITE(numout,*) ' ==>>> Initialising bi-Laplacian GM eddy viscosity coefficient.' 363 ALLOCATE( bhm(jpi,jpj,jpk) , STAT=ierr ) 364 bhm(:,:,:) = 0._wp ! init to 0 needed 365 bhm0 = rn_bhm_b ; ; cl_Units = ' m4/s' ! bi-Laplacian GM 366 SELECT CASE (nn_bhm_ijk_t) 367 !CW Define bi-Laplacian GM diffusivity and test the related computational stability criterion 368 369 CASE( 11 ) !== fixed profile: constant except for zero at top and bottom, so that eddy-induced velocity, w*=0 ==! 370 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity = constant in interior, zero at top and bottom' 371 IF(lwp) WRITE(numout,*) ' interior viscous coef. = constant = ', bhm0, cl_Units 372 373 ! First set to constant in interior, all levels 374 DO jj = 2, jpjm1 375 DO ji = 2, jpim1 376 bhm(ji,jj,:) = bhm0 377 ENDDO 378 ENDDO 379 380 ! Surface BC : set to zero 381 bhm(:,:,1) = 0._wp 382 ! Flat bottom BC : set to zero 383 bhm(:,:,jpk) = 0._wp 384 ! Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm 385 386 !CW Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping 387 ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32 388 ! test at T-point - may need to revise this choice!!!! 389 390 ! Find domain maximum value of LHS, then test inequality. 391 392 ! initialise 393 lhscritmax=0._wp 394 395 DO jj = 2, jpjm1 396 DO ji = 2, jpim1 397 ikw = mbkt(ji,jj) !bottom last wet T-level 398 DO jk = 2, ikw 399 lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2) 400 IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit 401 ENDDO 402 ENDDO 403 ENDDO 404 405 IF ( lhscritmax .ge. 1._wp/32._wp) THEN 406 IF(lwp) THEN 407 WRITE(numout,*) 408 WRITE(numout,*) ' WARNING : Bi-Laplacian GM eddy viscosity is not', & 409 & ' consistent with the model time step and grid setup: ', & 410 & ' LHS=',lhscrit, & 411 & 'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 412 WRITE(numout,*) ' =========' 413 WRITE(numout,*) 414 ENDIF 415 ! CW: warn, don't stop, as we may wish to violate this theoretical limit. 416 ! CALL ctl_stop( 'ldf_dyn_init', 'Inconsistent Bi-Lap. GM diffusivity') 417 ELSE 418 IF(lwp) THEN 419 WRITE(numout,*) 420 WRITE(numout,*) ' INFORMATION : Bi-Laplacian GM eddy viscosity is ', & 421 & ' consistent with the model time step and grid setup: ', & 422 & ' LHS=',lhscrit, & 423 & 'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)' 424 WRITE(numout,*) ' =========' 425 WRITE(numout,*) 426 ENDIF 427 ENDIF 428 !----------------------------- 429 CASE( 12 ) !== fixed profile: constant in interior, zero at bottom and linearly-tapering to zero at top, over depth shallower than rn_bgmzcrit ==! 430 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity = constant in interior, zero at bottom and' 431 IF(lwp) WRITE(numout,*) ' linearly-tapering to zero at top, over depth shallower than ', rn_bgmzcrit, ' metres.' 432 IF(lwp) WRITE(numout,*) ' Interior viscous coef. = constant = ', bhm0, cl_Units 433 434 ! Impose linear taper from interior value to zero at zero depth, for depths < rn_bgmzcrit 435 ! N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 436 437 DO jk = 1, jpk 438 439 IF (gdept_1d(jk) .lt. rn_bgmzcrit) THEN 440 441 DO jj = 2, jpjm1 442 DO ji = 2, jpim1 443 bhm(ji,jj,jk) = bhm0 * gdepw_1d(jk)/rn_bgmzcrit 444 ENDDO 445 ENDDO 446 447 ELSE 448 449 ! constant at depth rn_bgmzcrit and larger 450 DO jj = 2, jpjm1 451 DO ji = 2, jpim1 452 bhm(ji,jj,jk) = bhm0 453 ENDDO 454 ENDDO 455 456 ENDIF 457 458 ENDDO 459 460 ! Surface BC : set to zero 461 bhm(:,:,1) = 0._wp 462 ! Flat bottom BC : set to zero 463 bhm(:,:,jpk) = 0._wp 464 ! Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm 465 466 !-------------------------------- 467 CASE( 13 ) !== fixed profile: steady profile of the form bhm0 X dx^2 X dz^2 in interior, zero at bottom and top ==! 468 469 cl_Units = ' 1/s' 470 471 IF(lwp) WRITE(numout,*) ' ==>>> bi-Laplacian GM eddy viscosity : steady profile of the form' 472 IF(lwp) WRITE(numout,*) ' bhm0 X dx^2 X dz^2 in interior, zero boundary conds. at bottom and top.' 473 IF(lwp) WRITE(numout,*) ' In this case, bhm0 is the constant of proportionality,' 474 IF(lwp) WRITE(numout,*) ' dimensioned as [diffusivity]/L^4, i.e. bhm0=', bhm0, cl_Units 475 476 cl_Units = ' m4/s' 477 478 ! N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet. 479 480 DO jk = 1, jpk 481 ! constant at depth rn_bgmzcrit and larger 482 DO jj = 2, jpjm1 483 DO ji = 2, jpim1 484 bhm(ji,jj,jk) = bhm0*e1t(ji,jj)**2*e3w_n(ji,jj,jk)**2 485 ENDDO 486 ENDDO 487 ENDDO 488 489 ! Surface BC : set to zero 490 bhm(:,:,1) = 0._wp 491 ! Flat bottom BC : set to zero 492 bhm(:,:,jpk) = 0._wp 493 ! Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm 494 495 !-------------------------------- 496 497 CASE DEFAULT 498 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_bhm_ijk_t, the type of space-time variation of bhm') 499 END SELECT 500 ! 501 ! mask and lbc_lnk 502 bhm(:,:,1:jpkm1) = bhm(:,:,1:jpkm1) * wmask(:,:,1:jpkm1) 503 CALL lbc_lnk( 'ldf_dyn_init', bhm , 'W', 1. ) 504 ! 505 ENDIF ! ln_dynldf_bgm 529 506 ! 530 507 END SUBROUTINE ldf_dyn_init
Note: See TracChangeset
for help on using the changeset viewer.
