Changeset 12776
 Timestamp:
 20200420T11:24:26+02:00 (2 years ago)
 Location:
 NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE
 Files:

 2 edited
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 ) ! isolevel 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) ! biLaplacian GM parameterisation70 69 CASE ( np_blp ) ; CALL dyn_ldf_blp( kt, ub, vb, ua, va ) ! isolevel bilaplacian 71 70 ! 72 71 END SELECT 72 73 IF( ln_dynldf_bgm ) CALL dyn_ldf_bgm( kt, ub, vb, ua, va) ! biLaplacian 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,*) ' ==>>> isolevel laplacian operator' 107 108 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with isolevel background' 108 CASE( np_bgm ) ; WRITE(numout,*) ' ==>>> biLaplacian GM parameterisation via momentum equation'109 109 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> isolevel bilaplacian operator' 110 110 END SELECT 111 IF( ln_dynldf_bgm ) WRITE(numout,*) ' ==>>> biharmonic GentMcWilliams 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 !: isolevel operator 62 62 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: isoneutral or geopotential operator 63 INTEGER, PARAMETER, PUBLIC :: np_bgm = 30 !: isolevel operator, biLaplacian 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 ! biLaplacian GM operator233 IF( ln_zco ) THEN ! zcoordinate234 IF( ln_dynldf_lev ) nldf_dyn = np_bgm ! isolevel = horizontal (no rotation)235 IF( ln_dynldf_hor ) nldf_dyn = np_bgm ! isolevel = horizontal (no rotation)236 IF( ln_dynldf_iso ) ierr = 3 ! isoneutral ( rotation)237 ENDIF238 IF( ln_zps ) THEN ! zcoordinate with partial step239 IF( ln_dynldf_lev ) nldf_dyn = np_bgm ! isolevel (no rotation)240 IF( ln_dynldf_hor ) nldf_dyn = np_bgm ! isolevel (no rotation)241 IF( ln_dynldf_iso ) ierr = 3 ! isoneutral ( rotation)242 ENDIF243 IF( ln_sco ) THEN ! scoordinate244 IF( ln_dynldf_lev ) nldf_dyn = np_bgm ! isolevel (no rotation)245 IF( ln_dynldf_hor ) ierr = 3 ! horizontal ( rotation)246 IF( ln_dynldf_iso ) ierr = 3 ! isoneutral ( rotation)247 ENDIF248 ENDIF249 !250 230 251 231 IF( ierr == 2 ) CALL ctl_stop( 'rotated bilaplacian operator does not exist' ) 252 !253 IF( ierr == 3 ) CALL ctl_stop( 'rotated biLaplacian GM operator does not exist' )254 232 ! 255 233 … … 265 243 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with isolevel background' 266 244 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> isolevel bilaplacian operator' 267 CASE( np_bgm ) ; WRITE(numout,*) ' ==>>> isolevel biLaplacian 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' ! biLaplacian 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 spacetime variation of ahm') 373 345 END SELECT 374 !CW Add separate structure options for biLaplacian 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 biLaplacian 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 eddyinduced velocity, w*=0 ==! 381 IF(lwp) WRITE(numout,*) ' ==>>> biLaplacian 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 gridlength mode for BGM scheme with leapfrog timestepping 398 ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32 399 ! test at Tpoint  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 Tlevel 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 : BiLaplacian 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 BiLap. GM diffusivity') 428 ELSE 429 IF(lwp) THEN 430 WRITE(numout,*) 431 WRITE(numout,*) ' INFORMATION : BiLaplacian 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 linearlytapering to zero at top, over depth shallower than rn_bgmzcrit ==! 441 IF(lwp) WRITE(numout,*) ' ==>>> biLaplacian GM eddy viscosity = constant in interior, zero at bottom and' 442 IF(lwp) WRITE(numout,*) ' linearlytapering 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 gridlength 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,*) ' ==>>> biLaplacian 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 gridlength 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 spacetime 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 !biLaplacian 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 biLaplacian 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 biLaplacian 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' ! biLaplacian GM 366 SELECT CASE (nn_bhm_ijk_t) 367 !CW Define biLaplacian 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 eddyinduced velocity, w*=0 ==! 370 IF(lwp) WRITE(numout,*) ' ==>>> biLaplacian 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 gridlength mode for BGM scheme with leapfrog timestepping 387 ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32 388 ! test at Tpoint  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 Tlevel 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 : BiLaplacian 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 BiLap. GM diffusivity') 417 ELSE 418 IF(lwp) THEN 419 WRITE(numout,*) 420 WRITE(numout,*) ' INFORMATION : BiLaplacian 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 linearlytapering to zero at top, over depth shallower than rn_bgmzcrit ==! 430 IF(lwp) WRITE(numout,*) ' ==>>> biLaplacian GM eddy viscosity = constant in interior, zero at bottom and' 431 IF(lwp) WRITE(numout,*) ' linearlytapering 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 gridlength 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,*) ' ==>>> biLaplacian 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 gridlength 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 spacetime 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.