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 1559 for trunk/NEMO/OPA_SRC/TRA/traadv_qck.F90 – NEMO

Ignore:
Timestamp:
2009-07-29T16:03:14+02:00 (15 years ago)
Author:
ctlod
Message:

only cosmetic changes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r1528 r1559  
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :  9.0  !  2008-07  (G. Reffray)  Original code 
     6   !! History :  3.0  !  2008-07  (G. Reffray)  Original code 
    77   !!---------------------------------------------------------------------- 
    8  
    98 
    109   !!---------------------------------------------------------------------- 
    1110   !!   tra_adv_qck      : update the tracer trend with the horizontal advection 
    1211   !!                      trends using a 3rd order finite difference scheme 
    13    !!   tra_adv_qck_zon  :  
    14    !!   tra_adv_qck_mer  :  
    15    !!   tra_adv_cen2_ver : 2nd centered scheme for the vertical advection 
     12   !!   tra_adv_qck_i  :  
     13   !!   tra_adv_qck_j  :  
     14   !!   tra_adv_cen2_k : 2nd centered scheme for the vertical advection 
    1615   !!---------------------------------------------------------------------- 
    17    !! * Modules used 
    1816   USE oce             ! ocean dynamics and active tracers 
    1917   USE dom_oce         ! ocean space and time domain 
    … …  
    3129   PRIVATE 
    3230 
    33    !! * Accessibility 
    34    PUBLIC tra_adv_qck    ! routine called by step.F90 
    35  
    36    !! * Module variables 
    37    REAL(wp), DIMENSION(jpi,jpj), SAVE ::   btr2 
    38    REAL(wp)                    , SAVE ::   cst 
    39    !!---------------------------------------------------------------------- 
     31   PUBLIC   tra_adv_qck   ! routine called by step.F90 
     32 
     33   REAL(wp), DIMENSION(jpi,jpj) ::   btr2 
     34   REAL(wp)                     ::   r1_6 
     35 
    4036   !! * Substitutions 
    4137#  include "domzgr_substitute.h90" 
    4238#  include "vectopt_loop_substitute.h90" 
    4339   !!---------------------------------------------------------------------- 
    44    !!   OPA 9.0 , LOCEAN-IPSL (2008)  
     40   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    4541   !! $Id$ 
    4642   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    5753      !! 
    5854      !! ** Method :   The advection is evaluated by a third order scheme 
    59       !!               For a positive velocity u : 
    60       !! 
    61       !! 
    62       !!                  i-1    i      i+1   i+2 
    63       !! 
    64       !!               |--FU--|--FC--|--FD--|------| 
    65       !!                           u(i)>0 
    66       !! 
    67       !!               For a negative velocity u : 
    68       !! 
    69       !!               |------|--FD--|--FC--|--FU--| 
    70       !!                           u(i)<0 
    71       !! 
    72       !!                FU is the second upwind point 
    73       !!                FD is the first douwning point 
    74       !!                FC is the central point (or the first upwind point) 
    75       !! 
    76       !!      Flux(i) = u(i) * {0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC)} 
    77       !!                with C(i)=|u(i)|dx(i)/dt (Courant number) 
     55      !!             For a positive velocity u :              u(i)>0 
     56      !!                                          |--FU--|--FC--|--FD--|------| 
     57      !!                                             i-1    i      i+1   i+2 
     58      !! 
     59      !!             For a negative velocity u :              u(i)<0 
     60      !!                                          |------|--FD--|--FC--|--FU--| 
     61      !!                                             i-1    i      i+1   i+2 
     62      !!             where  FU is the second upwind point 
     63      !!                    FD is the first douwning point 
     64      !!                    FC is the central point (or the first upwind point) 
     65      !! 
     66      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) } 
     67      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number) 
    7868      !! 
    7969      !!         dt = 2*rdtra and the scalar values are tb and sb 
    … …  
    8171      !!       On the vertical, the simple centered scheme used tn and sn 
    8272      !! 
    83       !!               The fluxes are bounded by the ULTIMATE limiter 
    84       !!               to guarantee the monotonicity of the solution and to 
     73      !!               The fluxes are bounded by the ULTIMATE limiter to 
     74      !!             guarantee the monotonicity of the solution and to 
    8575      !!            prevent the appearance of spurious numerical oscillations 
    8676      !! 
    … …  
    9282      USE oce, ONLY :   ztrdt => ua   ! use ua as workspace 
    9383      USE oce, ONLY :   ztrds => va   ! use va as workspace 
    94       !! * Arguments 
     84      !! 
    9585      INTEGER , INTENT(in)                         ::  kt  ! ocean time-step index 
    9686      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pun ! effective ocean velocity, u_component 
    … …  
    10999         IF(lwp) WRITE(numout,*) 
    110100         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    111          cst       = 1. / 6. 
     101         r1_6      = 1. / 6. 
    112102      ENDIF 
    113103 
    … …  
    119109      !--------------------------------------------------------------------------- 
    120110 
    121       CALL tra_adv_qck_zon( kt, pun, tb, tn, ta, ztrdt, z2) 
    122       CALL tra_adv_qck_zon( kt, pun, sb, sn, sa, ztrds, z2) 
     111      CALL tra_adv_qck_i( pun, tb, tn, ta, ztrdt, z2) 
     112      CALL tra_adv_qck_i( pun, sb, sn, sa, ztrds, z2) 
    123113 
    124114      IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) 
    125115 
    126       CALL tra_adv_qck_mer( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2) 
    127       CALL tra_adv_qck_mer( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2)   
     116      CALL tra_adv_qck_j( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2) 
     117      CALL tra_adv_qck_j( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2)   
    128118 
    129119      IF( l_trdtra ) THEN 
    130120         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 
    131121         ! 
    132          ! Save the horizontal up-to-date ta/sa trends 
    133          ztrdt(:,:,:) = ta(:,:,:) 
     122         ztrdt(:,:,:) = ta(:,:,:)    ! Save the horizontal up-to-date ta/sa trends 
    134123         ztrds(:,:,:) = sa(:,:,:) 
    135124      END IF 
    … …  
    141130      !------------------------------------------------------------------------- 
    142131      ! 
    143       CALL tra_adv_cen2_ver( pwn, tn, ta ) 
    144       CALL tra_adv_cen2_ver( pwn, sn, sa ) 
     132      CALL tra_adv_cen2_k( pwn, tn, ta ) 
     133      CALL tra_adv_cen2_k( pwn, sn, sa ) 
    145134      ! 
    146135      !Save the vertical advective trends for diagnostic 
    … …  
    179168 
    180169 
    181    SUBROUTINE tra_adv_qck_zon ( kt, pun, tra, tran, traa, ztrdtra, z2 ) 
    182       !!---------------------------------------------------------------------- 
    183       !! 
    184       !!---------------------------------------------------------------------- 
    185       !! * Arguments 
    186       INTEGER,  INTENT(in)                            :: kt              ! ocean time-step index 
     170   SUBROUTINE tra_adv_qck_i ( pun, tra, tran, traa, ztrdtra, z2 ) 
     171      !!---------------------------------------------------------------------- 
     172      !! 
     173      !!---------------------------------------------------------------------- 
    187174      REAL,     INTENT(in)                            :: z2 
    188175      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran  ! horizontal effective velocity 
    … …  
    196183      REAL(wp), DIMENSION(jpi,jpj)     ::  zfu, zfc, zfd, zfho, zmskl, zsc_e  
    197184      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zflux 
    198       !---------------------------------------------------------------------- 
    199       ! I. Part 1 : x-direction 
    200185      !---------------------------------------------------------------------- 
    201186 
    … …  
    216201               ! Mask at the T-points in the x-direction (mask=0 or mask=1) 
    217202               zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2 
    218             ENDDO 
    219          ENDDO 
     203            END DO 
     204         END DO 
    220205         ! 
    221206         !--- Lateral boundary conditions  
    … …  
    244229               zfd(ji,jj)   = dir*tra  (ji+1,jj,jk)+(1-dir)*tra  (ji  ,jj,jk)  ! FD in the x-direction for T 
    245230               zmskl(ji,jj) = dir*zmask(ji  ,jj)   +(1-dir)*zmask(ji+1,jj) 
    246            ENDDO 
    247          ENDDO 
     231           END DO 
     232         END DO 
    248233         ! 
    249234         !--- QUICKEST scheme 
    … …  
    300285      END IF 
    301286 
    302    END SUBROUTINE tra_adv_qck_zon 
    303  
    304  
    305    SUBROUTINE tra_adv_qck_mer ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 ) 
    306       !!---------------------------------------------------------------------- 
    307       !! 
    308       !!---------------------------------------------------------------------- 
    309       !! * Arguments 
     287   END SUBROUTINE tra_adv_qck_i 
     288 
     289 
     290   SUBROUTINE tra_adv_qck_j ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 ) 
     291      !!---------------------------------------------------------------------- 
     292      !! 
     293      !!---------------------------------------------------------------------- 
    310294      INTEGER,  INTENT(in)                            :: kt              ! ocean time-step index 
    311295      REAL,     INTENT(in)                            :: z2 
    … …  
    314298      REAL(wp), INTENT(out)  , DIMENSION(jpi,jpj,jpk) :: ztrdtra 
    315299      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa 
    316       ! 
     300      !! 
    317301      INTEGER  :: ji, jj, jk 
    318302      REAL(wp) :: za, zbtr, dir, dx, dt     ! temporary scalars 
    … …  
    342326               ! Mask at the T-points in the x-direction (mask=0 or mask=1) 
    343327               zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 
    344             ENDDO 
    345          ENDDO 
     328            END DO 
     329         END DO 
    346330         ! 
    347331         !--- Lateral boundary conditions 
    … …  
    370354               zfd(ji,jj)   = dir*tra  (ji,jj+1,jk)+(1-dir)*tra  (ji,jj  ,jk)  ! FD in the y-direction for T 
    371355               zmskl(ji,jj) = dir*zmask(ji,jj     )+(1-dir)*zmask(ji,jj+1) 
    372             ENDDO 
    373          ENDDO 
     356            END DO 
     357         END DO 
    374358         ! 
    375359         !--- QUICKEST scheme 
    … …  
    440424      ENDIF 
    441425 
    442    END SUBROUTINE tra_adv_qck_mer 
    443  
    444  
    445    SUBROUTINE tra_adv_cen2_ver ( pwn, tra , traa ) 
    446       !!---------------------------------------------------------------------- 
    447       !! 
    448       !!---------------------------------------------------------------------- 
    449       !! * Arguments 
    450       REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk)  :: pwn        ! horizontal effective velocity 
    451       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk)  :: tra, traa 
    452       ! 
    453       INTEGER ji, jj, jk 
    454       REAL(wp) :: za, ze3tr       ! temporary scalars 
    455       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zflux 
    456       ! 
    457       ! Vertical advection 
    458       ! ------------------ 
    459       ! 
    460       ! 1. Vertical advective fluxes 
    461       ! ---------------------------- 
    462       ! 
    463       ! Bottom value : flux set to zero 
    464       zflux(:,:,jpk) = 0.e0 
    465       ! 
    466       ! Surface value 
    467       IF( lk_vvl ) THEN 
    468          ! variable volume : flux set to zero 
    469          zflux(:,:, 1 ) = 0.e0 
    470       ELSE 
    471          ! free surface-constant volume 
    472          zflux(:,:, 1 ) = pwn(:,:,1) * tra(:,:,1) 
     426   END SUBROUTINE tra_adv_qck_j 
     427 
     428 
     429   SUBROUTINE tra_adv_cen2_k ( pwn, ptn, pta ) 
     430      !!---------------------------------------------------------------------- 
     431      !! 
     432      !!---------------------------------------------------------------------- 
     433      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)  :: pwn   ! vertical effective velocity 
     434      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)  :: ptn   ! now tracer 
     435      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk)  :: pta   ! tracer general trend 
     436      !! 
     437      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     438      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zflux   ! 3D workspace 
     439      !!---------------------------------------------------------------------- 
     440      ! 
     441      !                         !==  Vertical advective fluxes  ==! 
     442      zflux(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero 
     443      ! 
     444      !                                 ! Surface value 
     445      IF( lk_vvl ) THEN   ;   zflux(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero 
     446      ELSE                ;   zflux(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)   ! Constant volume : advective flux through the surface 
    473447      ENDIF 
    474448      ! 
    475       ! Second order centered tracer flux at w-point 
    476       ! 
    477       DO jk = 2, jpkm1 
     449      DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point 
    478450         DO jj = 2, jpjm1 
    479451            DO ji = fs_2, fs_jpim1   ! vector opt. 
    480                zflux(ji,jj,jk) = pwn(ji,jj,jk)*0.5*( tra(ji,jj,jk-1) + tra(ji,jj,jk) ) 
     452               zflux(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1) + ptn(ji,jj,jk) ) 
    481453            END DO 
    482454         END DO 
    483455      END DO 
    484456      ! 
    485       ! 2. Tracer flux divergence at t-point added to the general trend 
    486       ! --------------------------------------------------------------- 
    487       ! 
    488       DO jk = 1, jpkm1 
     457      DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==! 
    489458         DO jj = 2, jpjm1 
    490459            DO ji = fs_2, fs_jpim1   ! vector opt. 
    491                ze3tr = 1. / fse3t(ji,jj,jk) 
    492                ! vertical advective trends 
    493                za = - ze3tr * ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) ) 
    494                ! add it to the general tracer trends 
    495                traa(ji,jj,jk) =  traa(ji,jj,jk) + za 
     460               pta(ji,jj,jk) =  pta(ji,jj,jk) - ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) )   & 
     461                  &                           /   fse3t(ji,jj,jk) 
    496462            END DO 
    497463         END DO 
    498464      END DO 
    499465      ! 
    500    END SUBROUTINE tra_adv_cen2_ver 
    501  
    502  
    503    SUBROUTINE quickest(fu,fd,fc,fho,fc_cfl) 
    504       !!---------------------------------------------------------------------- 
    505       !! 
    506       !!---------------------------------------------------------------------- 
    507       !! * Arguments 
     466   END SUBROUTINE tra_adv_cen2_k 
     467 
     468 
     469   SUBROUTINE quickest( fu, fd, fc, fho, fc_cfl ) 
     470      !!---------------------------------------------------------------------- 
     471      !! 
     472      !!---------------------------------------------------------------------- 
    508473      REAL(wp), INTENT(in)  , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl   
    509474      REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho 
    … …  
    513478      zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) ) 
    514479      zcoef2(:,:) = 0.5*fc_cfl(:,:)*( fd(:,:) - fc(:,:) ) 
    515       zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )*cst )*zcurv(:,:) 
     480      zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )*r1_6 )*zcurv(:,:) 
    516481      fho   (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:)                         ! phi_f QUICKEST  
    517482      ! 
Note: See TracChangeset for help on using the changeset viewer.