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

Ignore:
Timestamp:
2007-10-16T16:59:56+02:00 (17 years ago)
Author:
ctlod
Message:

get back to the nemo_v2_3 version for trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90

    • Property svn:keywords changed from Id to Author Date Id Revision
    r717 r719  
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv 
    7    !!             8.5  !  02-06  (G. Madec)  F90: Free form and module 
    8    !!             9.0  !  04-08  (C. Talandier) New trends organization 
    9    !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization 
    10    !!             " "  !  06-04  (R. Benshila, G. Madec) Step reorganization 
    11    !!             " "  !  06-07  (G. madec)  add ups_orca_set routine 
     6   !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv 
     7   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
     8   !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
     9   !!            " "  !  06-04  (R. Benshila, G. Madec) Step reorganization 
    1210   !!---------------------------------------------------------------------- 
    1311 
    1412   !!---------------------------------------------------------------------- 
    1513   !!   tra_adv_cen2_jki : update the tracer trend with the horizontal and 
    16    !!                      vertical advection trends using a seconder order 
    17    !!                      centered scheme. Auto-tasking case, k-slab for 
    18    !!                      hor. adv., j-slab for vert. adv. (j-k-i loops) 
     14   !!                  vertical advection trends using a seconder order 
     15   !!                  centered scheme. Auto-tasking case, k-slab for 
     16   !!                  hor. adv., j-slab for vert. adv. (j-k-i loops) 
    1917   !!---------------------------------------------------------------------- 
    2018   USE oce             ! ocean dynamics and active tracers 
    2119   USE dom_oce         ! ocean space and time domain 
    22    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient 
     20   USE trdmod          ! ocean active tracers trends 
    2321   USE trdmod_oce      ! ocean variables trends 
    24    USE sbc_oce         ! surface boundary condition: ocean 
    25    USE trdmod          ! ocean active tracers trends 
     22   USE flxrnf          ! 
    2623   USE trabbl          ! advective term in the BBL 
    2724   USE ocfzpt          ! 
    2825   USE lib_mpp 
    2926   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    30    USE sbcrnf          ! river runoffs 
    3127   USE in_out_manager  ! I/O manager 
    3228   USE diaptr          ! poleward transport diagnostics 
     29   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient 
    3330   USE prtctl          ! Print control 
    34     
     31 
    3532   IMPLICIT NONE 
    3633   PRIVATE 
    3734 
    38    PUBLIC   tra_adv_cen2_jki    ! routine called by step.F90 
    39  
    40    REAL(wp), DIMENSION(jpi,jpj) ::   btr2    ! inverse of T-point surface [1/(e1t*e2t)] 
     35   PUBLIC   tra_adv_cen2_jki   ! routine called by step.F90 
     36 
     37   REAL(wp), DIMENSION(jpi,jpj) ::   btr2   ! inverse of T-point surface [1/(e1t*e2t)] 
    4138 
    4239   !! * Substitutions 
     
    4542   !!---------------------------------------------------------------------- 
    4643   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
    47    !! $Id$ 
     44   !! $Header$  
    4845   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4946   !!---------------------------------------------------------------------- 
     
    9087      !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] ) 
    9188      !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  } 
    92       !!         C A U T I O N : the trend saved is the centered trend only. 
    93       !!      It does not take into account the upstream part of the scheme. 
    9489      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so  
    9590      !!      they vanish from the expression of the flux and divergence. 
     
    114109      !! 
    115110      !! ** Action :  - update (ta,sa) with the now advective tracer trends 
     111      !!              - save trends in (ztrdt,ztrds) ('key_trdtra') 
    116112      !!---------------------------------------------------------------------- 
    117       USE oce               ,   zwx => ua          ! use ua as workspace 
    118       USE oce               ,   zwy => va          ! use va as workspace 
    119       USE traadv_cen2, ONLY :   ups_orca_set,   &  ! upstream indicator near some straits  
    120          &                      upsmsk             ! and over closed sea (orca 2 and 4) 
    121  
    122       INTEGER , INTENT(in)                         ::   kt              ! ocean time-step index 
    123       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! now ocean velocity fields 
    124  
     113      USE oce, ONLY :   zwx => ua   ! use ua as workspace 
     114      USE oce, ONLY :   zwy => va   ! use va as workspace 
     115      !! 
     116      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index 
     117      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component 
     118      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component 
     119      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component 
     120      !! 
    125121      INTEGER  ::   ji, jj, jk                           ! dummy loop indices 
    126       REAL(wp) ::   zta, zsa, zbtr, zhw, ze3tr,       &  ! temporary scalars 
    127          &          zfp_ui, zfp_vj, zfp_w , zfui  ,   &  !    "         " 
    128          &          zfm_ui, zfm_vj, zfm_w , zfvj  ,   &  !    "         " 
    129          &          zcofi , zcofj , zcofk ,           &  !    "         " 
    130          &          zupsut, zupsus, zcenut, zcenus,   &  !    "         " 
    131          &          zupsvt, zupsvs, zcenvt, zcenvs,   &  !    "         " 
    132          &          zupst , zupss , zcent , zcens ,   &  !    "         " 
    133          &          z_hdivn_x, z_hdivn_y, z_hdivn  
    134       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind,   &  ! temporary workspace arrays 
    135          &                                  zww, ztrds             !    "         " 
     122      REAL(wp) ::                           & 
     123         zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars 
     124         zhw, ze3tr, zcofi, zcofj,          &  !    "         " 
     125         zupsut, zupsvt, zupsus, zupsvs,    &  !    "         " 
     126         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         " 
     127         zcofk, zupst, zupss, zcent,        &  !    "         " 
     128         zcens, zfp_w, zfm_w,               &  !    "         " 
     129         zcenut, zcenvt, zcenus, zcenvs,    &  !    "         " 
     130         z_hdivn_x, z_hdivn_y, z_hdivn 
     131      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace  
     132      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww, ztrds         !  "      " 
    136133      !!---------------------------------------------------------------------- 
    137       ! 
     134 
    138135      IF( kt == nit000 ) THEN 
    139136         IF(lwp) WRITE(numout,*) 
     
    141138         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   Auto-tasking case' 
    142139         IF(lwp) WRITE(numout,*) 
    143          !          
    144          upsmsk(:,:) = 0.e0                              ! not upstream by default 
    145          IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits 
    146          !                                               ! and in closed seas (orca2 and orca4 only) 
    147  
    148          btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )        ! inverse of T-point surface 
     140         ! 
     141         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    149142      ENDIF 
    150143 
     
    154147      DO jk = 1, jpkm1                                 ! Horizontal slab 
    155148         !                                             ! =============== 
    156  
    157          ! 0. Upstream / centered scheme indicator 
    158          ! --------------------------------------- 
     149         !  Upstream / centered scheme indicator 
     150         ! -------------------------------------- 
    159151         DO jj = 1, jpj 
    160152            DO ji = 1, jpi 
    161                zind(ji,jj,jk) = MAX (   & 
    162                   rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows) 
    163                   upsmsk(ji,jj)                      &  ! some of some straits 
     153               zind(ji,jj,jk) =  MAX (              & 
     154                  upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff 
     155                  upsadv(ji,jj)                     &  ! in the vicinity of some straits 
    164156#if defined key_ice_lim 
    165                   !                                     ! below ice covered area (if tn < "freezing"+0.1 ) 
    166                   , MAX(  0., SIGN( 1., fzptn(ji,jj) + 0.1 - tn(ji,jj,jk) )  ) * tmask(ji,jj,jk)   & 
    167 #endif 
    168                   &                  ) 
    169             END DO 
    170          END DO 
    171  
    172  
    173          ! I. Horizontal advective fluxes 
    174          ! ------------------------------ 
     157                  , tmask(ji,jj,jk)                 &  ! half upstream tracer fluxes if tn < ("freezing"+0.1 ) 
     158                  * MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) )   & 
     159#endif 
     160               ) 
     161            END DO 
     162         END DO 
     163 
     164         !  Horizontal advective fluxes 
     165         ! ----------------------------- 
    175166         ! Second order centered tracer flux at u and v-points 
    176167         DO jj = 1, jpjm1 
    177168            DO ji = 1, fs_jpim1   ! vector opt. 
    178                zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )               ! upstream indicator 
     169               ! upstream indicator 
     170               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) ) 
    179171               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) 
    180 #if defined key_zco 
    181                zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)                       ! volume fluxes * 1/2 (zco) 
     172               ! volume fluxes * 1/2 
     173#if defined key_zco 
     174               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk) 
    182175               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk) 
    183176#else 
    184                zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)     ! volume fluxes * 1/2 
     177               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
    185178               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    186179#endif 
    187                zfp_ui = zfui + ABS( zfui )                                   ! upstream scheme 
     180               ! upstream scheme 
     181               zfp_ui = zfui + ABS( zfui ) 
    188182               zfp_vj = zfvj + ABS( zfvj ) 
    189183               zfm_ui = zfui - ABS( zfui ) 
     
    193187               zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk) 
    194188               zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk) 
    195                zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )           ! centered scheme 
     189               ! centered scheme 
     190               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) ) 
    196191               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) ) 
    197192               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) ) 
    198193               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) ) 
    199                zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut          ! mixed centered / upstream scheme 
     194               ! mixed centered / upstream scheme 
     195               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut 
    200196               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt 
    201197               zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus 
     
    204200         END DO 
    205201 
    206          ! 2. Tracer flux divergence at t-point added to the general trend 
    207          ! --------------------------------------------------------------- 
    208  
     202         ! Tracer flux divergence at t-point added to the general trend 
    209203         DO jj = 2, jpjm1 
    210204            DO ji = fs_2, fs_jpim1   ! vector opt. 
    211205#if defined key_zco 
    212                zbtr = btr2(ji,jj)                                           ! inverse of the volume (zco) 
    213 #else 
    214                zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)                         ! inverse of the volume 
    215 #endif 
    216                zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &      ! horizontal advective trends 
     206               zbtr = btr2(ji,jj)  
     207#else 
     208               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 
     209#endif 
     210               ! horizontal advective trends 
     211               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
    217212                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    218213               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   & 
    219214                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  ) 
    220                ta(ji,jj,jk) = ta(ji,jj,jk) + zta                            ! add it to the general tracer trends 
     215               ! add it to the general tracer trends 
     216               ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
    221217               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
    222218            END DO 
     
    229225         &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    230226 
    231       ! 3. Save the horizontal advective trends for diagnostic 
     227      ! Save the horizontal advective trends for diagnostics 
    232228      IF( l_trdtra )   THEN 
    233229         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0 
     
    302298      ENDIF 
    303299 
    304       ! II. Vertical advection 
    305       ! ---------------------- 
     300      ! Vertical advection 
     301      ! -------------------- 
    306302!CDIR PARALLEL DO 
    307303!$OMP PARALLEL DO 
     
    316312         IF( lk_dynspg_rl ) THEN                          ! rigid lid : flux set to zero 
    317313            zwz(:,jj, 1 ) = 0.e0    ;    zww(:,jj, 1 ) = 0.e0 
    318          ELSE                                             ! free surface-constant volume : advection across the surface 
     314         ELSE                                             ! free surface-constant volume 
    319315            zwz(:,jj, 1 ) = pwn(:,jj,1) * tn(:,jj,1) 
    320316            zww(:,jj, 1 ) = pwn(:,jj,1) * sn(:,jj,1) 
    321317         ENDIF 
    322318          
    323          ! 1. Vertical advective fluxes (Second order centered tracer flux at w-point) 
    324          ! ---------------------------- 
     319         ! Vertical advective fluxes at w-point 
    325320         DO jk = 2, jpk 
    326321            DO ji = 2, jpim1 
    327                zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )          ! upstream indicator 
    328                zhw = 0.5 * pwn(ji,jj,jk)                                ! velocity * 1/2 
    329                zfp_w = zhw + ABS( zhw )                                 ! upstream scheme 
     322               ! upstream indicator 
     323               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
     324               ! velocity * 1/2 
     325               zhw = 0.5 * pwn(ji,jj,jk) 
     326               ! upstream scheme 
     327               zfp_w = zhw + ABS( zhw ) 
    330328               zfm_w = zhw - ABS( zhw ) 
    331329               zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1) 
    332330               zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1) 
    333                zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )          ! centered scheme 
     331               ! centered scheme 
     332               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) 
    334333               zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) 
    335                zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent       ! mixed centered / upstream scheme 
     334               ! mixed centered / upstream scheme 
     335               zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent 
    336336               zww(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens 
    337337            END DO 
    338338         END DO 
    339339 
    340          ! 2. Tracer flux divergence at t-point added to the general trend 
    341          ! ------------------------- 
     340         ! Tracer flux divergence at t-point added to the general trend 
    342341         DO jk = 1, jpkm1 
    343342            DO ji = 2, jpim1 
    344343               ze3tr = 1. / fse3t(ji,jj,jk) 
    345                zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )      ! vertical advective trends 
     344               ! vertical advective trends 
     345               zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
    346346               zsa = - ze3tr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) 
    347                ta(ji,jj,jk) =  ta(ji,jj,jk) + zta                       ! add it to the general tracer trends 
     347               ! add it to the general tracer trends 
     348               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta 
    348349               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa 
    349350            END DO 
     
    353354      !                                                ! =============== 
    354355 
    355       ! 3. Save the vertical advective trends for diagnostic 
    356       ! ---------------------------------------------------- 
     356      ! Save the vertical advective trends for diagnostic 
    357357      IF( l_trdtra )   THEN 
    358358         ! Recompute the vertical advection zta & zsa trends computed  
Note: See TracChangeset for help on using the changeset viewer.