New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 108 for trunk/NEMO/OPA_SRC/DYN – NEMO

Changeset 108 for trunk/NEMO/OPA_SRC/DYN


Ignore:
Timestamp:
2004-06-28T12:19:07+02:00 (20 years ago)
Author:
opalod
Message:

CT : UPDATE069 : Vorticity diagnostics and new advection scheme (energy and enstrophy conserving) have been added

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/DYN/dynvor.F90

    r106 r108  
    1010   !!   dyn_vor_energy   : energy conserving scheme          (ln_dynvor_ene=T) 
    1111   !!   dyn_vor_mixed    : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 
     12   !!   dyn_vor_ene_ens  : energy and enstrophy conserving   (ln_dynvor_een=T) 
    1213   !!   dyn_vor_ctl      : control of the different vorticity option 
    1314   !!---------------------------------------------------------------------- 
     
    1718   USE in_out_manager ! I/O manager 
    1819   USE trddyn_oce     ! ocean momentum trends 
     20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1921 
    2022   IMPLICIT NONE 
     
    2527   PUBLIC dyn_vor_energy         ! routine called by step.F90 
    2628   PUBLIC dyn_vor_mixed          ! routine called by step.F90 
     29   PUBLIC dyn_vor_ene_ens        ! routine called by step.F90 
    2730   PUBLIC dyn_vor_ctl            ! routine called by step.F90 
    2831 
     
    3134   LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme 
    3235   LOGICAL, PUBLIC ::   ln_dynvor_mix = .FALSE.   !: mixed scheme 
     36   LOGICAL, PUBLIC ::   ln_dynvor_een = .TRUE.    !: energy and enstrophy conserving scheme 
    3337 
    3438   !! * Substitutions 
     
    8286      REAL(wp), DIMENSION(jpi,jpj) ::   & 
    8387         zwx, zwy, zwz                      ! temporary workspace 
    84 #if defined key_trddyn 
     88#if defined key_trddyn || defined key_trd_vor 
    8589      REAL(wp) ::   & 
    8690         zcu, zcv, zce3                     !    "         " 
     
    125129               ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    126130               va(ji,jj,jk) = va(ji,jj,jk) + zva 
    127 #   if defined key_trddyn 
     131#   if defined key_trddyn || defined key_trd_vor 
    128132#      if defined key_s_coord 
    129133               zce3= ff(ji,jj) / fse3f(ji,jj,jk) 
     
    248252               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua 
    249253               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva 
    250 #   if defined key_trddyn 
     254#   if defined key_trddyn || defined key_trd_vor 
    251255               utrd(ji,jj,jk,3) = zua 
    252256               vtrd(ji,jj,jk,3) = zva 
     
    313317      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    314318         zwz                           ! temporary workspace 
    315 #   if defined key_trddyn 
     319#   if defined key_trddyn || defined key_trd_vor 
    316320      REAL(wp) ::   & 
    317321         zcu, zcv, zce3                ! temporary scalars  
     
    368372               va(ji,jj,jk) = va(ji,jj,jk) + zva 
    369373 
    370 #   if defined key_trddyn 
     374#   if defined key_trddyn || defined key_trd_vor 
    371375#      if defined key_s_coord 
    372376               zce3 = ff(ji,jj) / fse3f(ji,jj,jk) 
     
    404408 
    405409 
     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 
    406544   SUBROUTINE dyn_vor_ctl 
    407545      !!--------------------------------------------------------------------- 
     
    417555      INTEGER ::   ioptio = 0      ! temporary integer 
    418556 
    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 
    420558      !!---------------------------------------------------------------------- 
    421559 
     
    436574         WRITE(numout,*) '                 energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene 
    437575         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 
    438577      ENDIF 
    439578 
     
    451590         IF(lwp) WRITE(numout,*) 
    452591         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' 
    453597         ioptio = ioptio + 1 
    454598      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.