Changeset 108 for trunk/NEMO/OPA_SRC/DYN/dynvor.F90
- Timestamp:
- 2004-06-28T12:19:07+02:00 (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynvor.F90
r106 r108 10 10 !! dyn_vor_energy : energy conserving scheme (ln_dynvor_ene=T) 11 11 !! dyn_vor_mixed : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 12 !! dyn_vor_ene_ens : energy and enstrophy conserving (ln_dynvor_een=T) 12 13 !! dyn_vor_ctl : control of the different vorticity option 13 14 !!---------------------------------------------------------------------- … … 17 18 USE in_out_manager ! I/O manager 18 19 USE trddyn_oce ! ocean momentum trends 20 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 19 21 20 22 IMPLICIT NONE … … 25 27 PUBLIC dyn_vor_energy ! routine called by step.F90 26 28 PUBLIC dyn_vor_mixed ! routine called by step.F90 29 PUBLIC dyn_vor_ene_ens ! routine called by step.F90 27 30 PUBLIC dyn_vor_ctl ! routine called by step.F90 28 31 … … 31 34 LOGICAL, PUBLIC :: ln_dynvor_ens = .TRUE. !: enstrophy conserving scheme 32 35 LOGICAL, PUBLIC :: ln_dynvor_mix = .FALSE. !: mixed scheme 36 LOGICAL, PUBLIC :: ln_dynvor_een = .TRUE. !: energy and enstrophy conserving scheme 33 37 34 38 !! * Substitutions … … 82 86 REAL(wp), DIMENSION(jpi,jpj) :: & 83 87 zwx, zwy, zwz ! temporary workspace 84 #if defined key_trddyn 88 #if defined key_trddyn || defined key_trd_vor 85 89 REAL(wp) :: & 86 90 zcu, zcv, zce3 ! " " … … 125 129 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 126 130 va(ji,jj,jk) = va(ji,jj,jk) + zva 127 # if defined key_trddyn 131 # if defined key_trddyn || defined key_trd_vor 128 132 # if defined key_s_coord 129 133 zce3= ff(ji,jj) / fse3f(ji,jj,jk) … … 248 252 ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua 249 253 va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva 250 # if defined key_trddyn 254 # if defined key_trddyn || defined key_trd_vor 251 255 utrd(ji,jj,jk,3) = zua 252 256 vtrd(ji,jj,jk,3) = zva … … 313 317 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 314 318 zwz ! temporary workspace 315 # if defined key_trddyn 319 # if defined key_trddyn || defined key_trd_vor 316 320 REAL(wp) :: & 317 321 zcu, zcv, zce3 ! temporary scalars … … 368 372 va(ji,jj,jk) = va(ji,jj,jk) + zva 369 373 370 # if defined key_trddyn 374 # if defined key_trddyn || defined key_trd_vor 371 375 # if defined key_s_coord 372 376 zce3 = ff(ji,jj) / fse3f(ji,jj,jk) … … 404 408 405 409 410 SUBROUTINE dyn_vor_ene_ens( kt ) 411 !!---------------------------------------------------------------------- 412 !! *** ROUTINE dyn_vor_ene_ens *** 413 !! 414 !! ** Purpose : Compute the now total vorticity trend and add it to 415 !! the general trend of the momentum equation. 416 !! 417 !! ** Method : Trend evaluated using now fields (centered in time) 418 !! and the Arakawa and Lamb (19XX) flux form formulation : conserves 419 !! both the horizontal kinetic energy and the potential enstrophy 420 !! when horizontal divergence is zero. 421 !! The trend of the vorticity term is given by: 422 !! * s-coordinate (lk_sco=T), the e3. are inside the derivatives: 423 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 424 !! Add this trend to the general momentum trend (ua,va): 425 !! (ua,va) = (ua,va) + ( voru , vorv ) 426 !! 427 !! ** Action : - Update (ua,va) with the now vorticity term trend 428 !! - save the trends in (utrd,vtrd) in 2 parts (relative 429 !! and planetary vorticity trends) ('key_trddyn') 430 !! 431 !! References : 432 !! Arakawa and Lamb 19XX, ??? 433 !! History : 434 !! 5.0 ! 04-02 (G. Madec) Original code 435 !!---------------------------------------------------------------------- 436 !! * Arguments 437 INTEGER, INTENT( in ) :: kt ! ocean time-step index 438 439 !! * Local declarations 440 INTEGER :: ji, jj, jk ! dummy loop indices 441 REAL(wp) :: & 442 zfac12, zua, zva ! temporary scalars 443 REAL(wp), DIMENSION(jpi,jpj) :: & 444 zwx, zwy, zwz, & ! temporary workspace 445 ztnw, ztne, ztsw, ztse ! " " 446 #if defined key_trddyn || defined key_trd_vor 447 REAL(wp) :: & 448 zcu, zcv ! " " 449 REAL(wp), DIMENSION(jpi,jpj) :: & 450 zcor ! potential planetary vorticity (f/e3) 451 #endif 452 REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: & 453 ze3f 454 !!---------------------------------------------------------------------- 455 456 IF( kt == nit000 ) THEN 457 IF(lwp) WRITE(numout,*) 458 IF(lwp) WRITE(numout,*) 'dyn_vor_ene_ens : vorticity term: energy and enstrophy conserving scheme' 459 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 460 461 DO jk = 1, jpk 462 DO jj = 1, jpjm1 463 DO ji = 1, jpim1 464 ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 465 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) * 0.25 466 !!! ze3f(ji,jj,jk) = MAX( ze3f(ji,jj,jk) , 1.e-20) 467 IF( ze3f(ji,jj,jk) /= 0.e0 ) ze3f(ji,jj,jk) = 1.e0 / ze3f(ji,jj,jk) 468 END DO 469 END DO 470 END DO 471 CALL lbc_lnk( ze3f, 'F', 1. ) 472 ENDIF 473 474 ! Local constant initialization 475 zfac12 = 1.e0 / 12.e0 476 477 ! ! =============== 478 DO jk = 1, jpkm1 ! Horizontal slab 479 ! ! =============== 480 481 ! Potential vorticity and horizontal fluxes 482 ! ----------------------------------------- 483 !!!bug zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk) 484 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 485 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 486 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 487 #if defined key_trddyn || defined key_trd_vor 488 zcor(:,:) = ff(:,:) * ze3f(:,:,jk) 489 #endif 490 491 ! Compute and add the vorticity term trend 492 ! ---------------------------------------- 493 jj=2 494 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 495 DO ji = 2, jpi 496 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 497 ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 498 ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 499 ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 500 END DO 501 DO jj = 3, jpj 502 DO ji = fs_2, jpi ! vector opt. 503 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 504 ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 505 ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 506 ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 507 END DO 508 END DO 509 510 DO jj = 2, jpjm1 511 DO ji = fs_2, fs_jpim1 ! vector opt. 512 zua = + zfac12 / e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 513 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 514 zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 515 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 516 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 517 va(ji,jj,jk) = va(ji,jj,jk) + zva 518 # if defined key_trddyn || defined key_trd_vor 519 zcu = + zfac12 / e1u(ji,jj) * ( zcor(ji,jj ) * zwy(ji ,jj ) + zcor(ji+1,jj) * zwy(ji+1,jj ) & 520 & + zcor(ji,jj ) * zwy(ji ,jj-1) + zcor(ji+1,jj) * zwy(ji+1,jj-1) ) 521 zcv = - zfac12 / e2v(ji,jj) * ( zcor(ji,jj+1) * zwx(ji-1,jj+1) + zcor(ji,jj+1) * zwx(ji ,jj+1) & 522 & + zcor(ji,jj ) * zwx(ji-1,jj ) + zcor(ji,jj ) * zwx(ji ,jj ) ) 523 utrd(ji,jj,jk,3) = zua - zcu 524 vtrd(ji,jj,jk,3) = zva - zcv 525 utrd(ji,jj,jk,4) = zcu 526 vtrd(ji,jj,jk,4) = zcv 527 # endif 528 END DO 529 END DO 530 ! ! =============== 531 END DO ! End of slab 532 ! ! =============== 533 534 IF(l_ctl) THEN ! print sum trends (used for debugging) 535 zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) ) 536 zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) ) 537 WRITE(numout,*) ' vor een - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl 538 u_ctl = zua ; v_ctl = zva 539 ENDIF 540 541 END SUBROUTINE dyn_vor_ene_ens 542 543 406 544 SUBROUTINE dyn_vor_ctl 407 545 !!--------------------------------------------------------------------- … … 417 555 INTEGER :: ioptio = 0 ! temporary integer 418 556 419 NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix 557 NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een 420 558 !!---------------------------------------------------------------------- 421 559 … … 436 574 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 437 575 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 576 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 438 577 ENDIF 439 578 … … 451 590 IF(lwp) WRITE(numout,*) 452 591 IF(lwp) WRITE(numout,*) ' vorticity term : mixed enstrophy/energy conserving scheme' 592 ioptio = ioptio + 1 593 ENDIF 594 IF( ln_dynvor_een ) THEN 595 IF(lwp) WRITE(numout,*) 596 IF(lwp) WRITE(numout,*) ' vorticity term : energy and enstrophy conserving scheme' 453 597 ioptio = ioptio + 1 454 598 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.