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 10806 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

Ignore:
Timestamp:
2019-03-27T17:55:22+01:00 (5 years ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch: Latest updates. Make sure all time-dependent 3D variables are converted in previously modified modules.

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90

    r10802 r10806  
    7575CONTAINS 
    7676 
    77    SUBROUTINE tra_adv( kt, ktlev1, ktlev2, ktlev3, pts_rhs ) 
     77   SUBROUTINE tra_adv( kt, ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs ) 
    7878      !!---------------------------------------------------------------------- 
    7979      !!                  ***  ROUTINE tra_adv  *** 
    … …  
    8484      !!---------------------------------------------------------------------- 
    8585      INTEGER, INTENT(in) ::   kt                       ! ocean time-step index 
    86       INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for source terms 
     86      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms 
     87      INTEGER, INTENT(in) ::   kt2lev                   ! time level index for 2-time-level source terms 
    8788      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    8889      ! 
    … …  
    153154     &                         ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts, nn_fct_h, nn_fct_v ) 
    154155      CASE ( np_MUS )                                 ! MUSCL 
    155          CALL tra_adv_mus    ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1),      pts_rhs, jpts        , ln_mus_ups )  
     156         CALL tra_adv_mus    ( kt, nit000, ktlev2, kt2lev, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1),      pts_rhs, jpts        , ln_mus_ups )  
    156157      CASE ( np_UBS )                                 ! UBS 
    157158         CALL tra_adv_ubs    ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts        , nn_ubs_v   ) 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90

    r10802 r10806  
    4444CONTAINS 
    4545 
    46    SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pwn,     & 
     46   SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pw,     & 
    4747      &                                               pt, pt_rhs, kjpt, kn_cen_h, kn_cen_v )  
    4848      !!---------------------------------------------------------------------- 
    … …  
    7070      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme) 
    7171      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme) 
    72       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     72      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    7373      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt             ! now tracer fields 
    7474      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
    … …  
    151151               DO jj = 2, jpjm1 
    152152                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    153                      zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
     153                     zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
    154154                  END DO 
    155155               END DO 
    … …  
    161161               DO jj = 2, jpjm1 
    162162                  DO ji = fs_2, fs_jpim1 
    163                      zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     163                     zwz(ji,jj,jk) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
    164164                  END DO 
    165165               END DO 
    … …  
    172172               DO jj = 1, jpj 
    173173                  DO ji = 1, jpi 
    174                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)  
     174                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)  
    175175                  END DO 
    176176               END DO    
    177177            ELSE                                   ! no ice-shelf cavities (only ocean surface) 
    178                zwz(:,:,1) = pwn(:,:,1) * pt(:,:,1,jn) 
     178               zwz(:,:,1) = pw(:,:,1) * pt(:,:,1,jn) 
    179179            ENDIF 
    180180         ENDIF 
    … …  
    194194            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu, pt(:,:,:,jn) ) 
    195195            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv, pt(:,:,:,jn) ) 
    196             CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, pt(:,:,:,jn) ) 
     196            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, pt(:,:,:,jn) ) 
    197197         END IF 
    198198         !                                 ! "Poleward" heat and salt transports  

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90

    r10802 r10806  
    5252CONTAINS 
    5353 
    54    SUBROUTINE tra_adv_fct( kt, kit000, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pwn,       & 
     54   SUBROUTINE tra_adv_fct( kt, kit000, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pw,       & 
    5555      &                                                                pt_lev1, pt_lev2, pt_rhs, kjpt, kn_fct_h, kn_fct_v ) 
    5656      !!---------------------------------------------------------------------- 
    … …  
    7777      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
    7878      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    79       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     79      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    8080      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields 
    8181      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
    … …  
    139139            DO jj = 1, jpj 
    140140               DO ji = 1, jpi 
    141                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    142                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     141                  zfp_wk = pw(ji,jj,jk) + ABS( pw(ji,jj,jk) ) 
     142                  zfm_wk = pw(ji,jj,jk) - ABS( pw(ji,jj,jk) ) 
    143143                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt_lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
    144144               END DO 
    … …  
    149149               DO jj = 1, jpj 
    150150                  DO ji = 1, jpi 
    151                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     151                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    152152                  END DO 
    153153               END DO    
    154154            ELSE                             ! no cavities: only at the ocean surface 
    155                zwz(:,:,1) = pwn(:,:,1) * pt_lev1(:,:,1,jn) 
     155               zwz(:,:,1) = pw(:,:,1) * pt_lev1(:,:,1,jn) 
    156156            ENDIF 
    157157         ENDIF 
    … …  
    258258               DO jj = 2, jpjm1 
    259259                  DO ji = fs_2, fs_jpim1 
    260                      zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
     260                     zwz(ji,jj,jk) =  (  pw(ji,jj,jk) * 0.5_wp * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
    261261                        &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
    262262                  END DO 
    … …  
    269269               DO jj = 2, jpjm1 
    270270                  DO ji = fs_2, fs_jpim1 
    271                      zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
     271                     zwz(ji,jj,jk) = ( pw(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
    272272                  END DO 
    273273               END DO 
    … …  
    306306               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pu, pt_lev2(:,:,:,jn) ) 
    307307               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pv, pt_lev2(:,:,:,jn) ) 
    308                CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, pt_lev2(:,:,:,jn) ) 
     308               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pw, pt_lev2(:,:,:,jn) ) 
    309309            ENDIF 
    310310            !                             ! heat/salt transport 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90

    r10802 r10806  
    5454CONTAINS 
    5555 
    56    SUBROUTINE tra_adv_mus( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,             & 
    57       &                                                     pt, pt_rhs, kjpt, ld_msc_ups ) 
     56   SUBROUTINE tra_adv_mus( kt, kit000, ktlev, kt2lev, cdtype, p2dt, pu, pv, pw,             & 
     57      &                                                             pt, pt_rhs, kjpt, ld_msc_ups ) 
    5858      !!---------------------------------------------------------------------- 
    5959      !!                    ***  ROUTINE tra_adv_mus  *** 
    … …  
    7575      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    7676      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    77       INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms 
     77      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for 3-time-level source terms 
     78      INTEGER                              , INTENT(in   ) ::   kt2lev          ! time level index for 2-time-level source terms 
    7879      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    7980      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    8081      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl 
    8182      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    82       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     83      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    8384      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt             ! before tracer field 
    8485      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
    … …  
    240241            DO jj = 2, jpjm1       
    241242               DO ji = fs_2, fs_jpim1   ! vector opt. 
    242                   z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
     243                  z0w = SIGN( 0.5, pw(ji,jj,jk+1) ) 
    243244                  zalpha = 0.5 + z0w 
    244                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 
     245                  zw  = z0w - 0.5 * pw(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,kt2lev) 
    245246                  zzwx = pt(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    246247                  zzwy = pt(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
    247                   zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
     248                  zwx(ji,jj,jk+1) = pw(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
    248249               END DO  
    249250            END DO 
    … …  
    253254               DO jj = 1, jpj 
    254255                  DO ji = 1, jpi 
    255                      zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn) 
     256                     zwx(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn) 
    256257                  END DO 
    257258               END DO    
    258259            ELSE                                      ! no cavities: only at the ocean surface 
    259                zwx(:,:,1) = pwn(:,:,1) * pt(:,:,1,jn) 
     260               zwx(:,:,1) = pw(:,:,1) * pt(:,:,1,jn) 
    260261            ENDIF 
    261262         ENDIF 
    … …  
    269270         END DO 
    270271         !                                ! send trends for diagnostic 
    271          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, pt(:,:,:,jn) ) 
     272         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pw, pt(:,:,:,jn) ) 
    272273         ! 
    273274      END DO                     ! end of tracer loop 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90

    r10802 r10806  
    4747CONTAINS 
    4848 
    49    SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,      & 
     49   SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw,      & 
    5050      &                                                pt_lev1, pt_lev2, pt_rhs, kjpt ) 
    5151      !!---------------------------------------------------------------------- 
    … …  
    9090      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    9191      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    92       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    9393      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields 
    9494      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
    … …  
    113113 
    114114      !        ! vertical fluxes are computed with the 2nd order centered scheme 
    115       CALL tra_adv_cen2_k( kt, ktlev, cdtype, pwn,         pt_lev2, pt_rhs, kjpt ) 
     115      CALL tra_adv_cen2_k( kt, ktlev, cdtype, pw,         pt_lev2, pt_rhs, kjpt ) 
    116116      ! 
    117117   END SUBROUTINE tra_adv_qck 
    … …  
    351351 
    352352 
    353    SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pwn,           & 
     353   SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pw,           & 
    354354     &                                    pt_lev2, pt_rhs, kjpt ) 
    355355      !!---------------------------------------------------------------------- 
    … …  
    360360      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    361361      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
    362       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity  
     362      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pw      ! vertical velocity  
    363363      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev2      ! before and now tracer fields 
    364364      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs      ! tracer trend  
    … …  
    378378            DO jj = 2, jpjm1 
    379379               DO ji = fs_2, fs_jpim1   ! vector opt. 
    380                   zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
     380                  zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
    381381               END DO 
    382382            END DO 
    … …  
    386386               DO jj = 1, jpj 
    387387                  DO ji = 1, jpi 
    388                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     388                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    389389                  END DO 
    390390               END DO    
    391391            ELSE                                   ! no ocean cavities (only ocean surface) 
    392                zwz(:,:,1) = pwn(:,:,1) * pt_lev2(:,:,1,jn) 
     392               zwz(:,:,1) = pw(:,:,1) * pt_lev2(:,:,1,jn) 
    393393            ENDIF 
    394394         ENDIF 
    … …  
    403403         END DO 
    404404         !                                 ! Send trends for diagnostic 
    405          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, pt_lev2(:,:,:,jn) ) 
     405         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, pt_lev2(:,:,:,jn) ) 
    406406         ! 
    407407      END DO 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90

    r10802 r10806  
    4646CONTAINS 
    4747 
    48    SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,          & 
     48   SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw,          & 
    4949      &                                                pt_lev1, pt_lev2, pt_rhs, kjpt, kn_ubs_v ) 
    5050      !!---------------------------------------------------------------------- 
    … …  
    9191      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
    9292      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    93       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean transport components 
     93      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean transport components 
    9494      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields 
    9595      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
    … …  
    200200               DO jj = 1, jpj 
    201201                  DO ji = 1, jpi 
    202                      zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    203                      zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     202                     zfp_wk = pw(ji,jj,jk) + ABS( pw(ji,jj,jk) ) 
     203                     zfm_wk = pw(ji,jj,jk) - ABS( pw(ji,jj,jk) ) 
    204204                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt_lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk) 
    205205                  END DO 
    … …  
    210210                  DO jj = 1, jpj 
    211211                     DO ji = 1, jpi 
    212                         ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     212                        ztw(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    213213                     END DO 
    214214                  END DO    
    215215               ELSE                                ! no cavities: only at the ocean surface 
    216                   ztw(:,:,1) = pwn(:,:,1) * pt_lev1(:,:,1,jn) 
     216                  ztw(:,:,1) = pw(:,:,1) * pt_lev1(:,:,1,jn) 
    217217               ENDIF 
    218218            ENDIF 
    … …  
    233233               DO jj = 1, jpj 
    234234                  DO ji = 1, jpi 
    235                      ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
     235                     ztw(ji,jj,jk) = (   0.5_wp * pw(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
    236236                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
    237237                  END DO 
    … …  
    248248               DO jj = 2, jpjm1 
    249249                  DO ji = fs_2, fs_jpim1 
    250                      ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
    251                   END DO 
    252                END DO 
    253             END DO 
    254             IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * pt_lev2(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
     250                     ztw(ji,jj,jk) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     251                  END DO 
     252               END DO 
     253            END DO 
     254            IF( ln_linssh )   ztw(:,:, 1 ) = pw(:,:,1) * pt_lev2(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
    255255            ! 
    256256         END SELECT 
    … …  
    269269                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    270270                     zltv(ji,jj,jk) = pt_rhs(ji,jj,jk,jn) - zltv(ji,jj,jk)                          & 
    271                         &           + pt_lev2(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   & 
     271                        &           + pt_lev2(ji,jj,jk,jn) * (  pw(ji,jj,jk) - pw(ji,jj,jk+1)  )   & 
    272272                        &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev) 
    273273                  END DO 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90

    r10425 r10806  
    5151CONTAINS 
    5252 
    53    SUBROUTINE tra_bbc( kt ) 
     53   SUBROUTINE tra_bbc( kt, ktlev, pts_rhs ) 
    5454      !!---------------------------------------------------------------------- 
    5555      !!                  ***  ROUTINE tra_bbc  *** 
    … …  
    7474      !!---------------------------------------------------------------------- 
    7575      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     76      INTEGER, INTENT(in) ::   ktlev  ! time level index for source terms 
     77      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    7678      ! 
    7779      INTEGER  ::   ji, jj    ! dummy loop indices 
    … …  
    8385      IF( l_trdtra )   THEN         ! Save the input temperature trend 
    8486         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
    85          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     87         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
    8688      ENDIF 
    8789      !                             !  Add the geothermal trend on temperature 
    8890      DO jj = 2, jpjm1 
    8991         DO ji = 2, jpim1 
    90             tsa(ji,jj,mbkt(ji,jj),jp_tem) = tsa(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t_n(ji,jj,mbkt(ji,jj)) 
     92            pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) = pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),ktlev) 
    9193         END DO 
    9294      END DO 
    9395      ! 
    94       CALL lbc_lnk( 'trabbc', tsa(:,:,:,jp_tem) , 'T', 1. ) 
     96      CALL lbc_lnk( 'trabbc', pts_rhs(:,:,:,jp_tem) , 'T', 1. ) 
    9597      ! 
    9698      IF( l_trdtra ) THEN        ! Send the trend for diagnostics 
    97          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
     99         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
    98100         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 
    99101         DEALLOCATE( ztrdt ) 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90

    r10425 r10806  
    8989 
    9090 
    91    SUBROUTINE tra_bbl( kt ) 
     91   SUBROUTINE tra_bbl( kt, ktlev1, ktlev2, kt2lev, pts_rhs ) 
    9292      !!---------------------------------------------------------------------- 
    9393      !!                  ***  ROUTINE bbl  *** 
    … …  
    101101      !!              is added to the general tracer trend 
    102102      !!---------------------------------------------------------------------- 
    103       INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
     103      INTEGER, INTENT( in ) ::   kt              ! ocean time-step 
     104      INTEGER, INTENT( in ) ::   ktlev1, ktlev2  ! time level indices for 3-time-level source terms 
     105      INTEGER, INTENT( in ) ::   kt2lev          ! time level index for 2-time-level source terms 
     106      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    104107      ! 
    105108      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds 
    … …  
    110113      IF( l_trdtra )   THEN                         !* Save the T-S input trends 
    111114         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    112          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    113          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    114       ENDIF 
    115  
    116       IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
     115         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
     116         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 
     117      ENDIF 
     118 
     119      IF( l_bbl )   CALL bbl( kt, nit000, ktlev1, ktlev2, kt2lev, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
    117120 
    118121      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl 
    119122         ! 
    120          CALL tra_bbl_dif( tsb, tsa, jpts ) 
     123         CALL tra_bbl_dif( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts ) 
    121124         IF( ln_ctl )  & 
    122125         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, & 
    … …  
    131134      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl 
    132135         ! 
    133          CALL tra_bbl_adv( tsb, tsa, jpts ) 
     136         CALL tra_bbl_adv( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts ) 
    134137         IF(ln_ctl)   & 
    135138         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   & 
    … …  
    143146 
    144147      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics 
    145          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    146          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     148         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
     149         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:) 
    147150         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 
    148151         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 
    … …  
    155158 
    156159 
    157    SUBROUTINE tra_bbl_dif( ptb, pta, kjpt ) 
     160   SUBROUTINE tra_bbl_dif( pt, pe3t, pt_rhs, kjpt ) 
    158161      !!---------------------------------------------------------------------- 
    159162      !!                  ***  ROUTINE tra_bbl_dif  *** 
    … …  
    171174      !!      convection is satified) 
    172175      !! 
    173       !! ** Action  :   pta   increased by the bbl diffusive trend 
     176      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend 
    174177      !! 
    175178      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 
    … …  
    177180      !!---------------------------------------------------------------------- 
    178181      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers 
    179       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields 
    180       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend 
     182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! tracer fields 
     183      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields 
     184      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend 
    181185      ! 
    182186      INTEGER  ::   ji, jj, jn   ! dummy loop indices 
    … …  
    191195            DO ji = 1, jpi 
    192196               ik = mbkt(ji,jj)                             ! bottom T-level index 
    193                zptb(ji,jj) = ptb(ji,jj,ik,jn)               ! bottom before T and S 
     197               zptb(ji,jj) = pt(ji,jj,ik,jn)               ! bottom before T and S 
    194198            END DO 
    195199         END DO 
    … …  
    198202            DO ji = 2, jpim1 
    199203               ik = mbkt(ji,jj)                            ! bottom T-level index 
    200                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  & 
     204               pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  & 
    201205                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     & 
    202206                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     & 
    203207                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     & 
    204208                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  & 
    205                   &             * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik) 
     209                  &             * r1_e1e2t(ji,jj) / pe3t(ji,jj,ik) 
    206210            END DO 
    207211         END DO 
    … …  
    212216 
    213217 
    214    SUBROUTINE tra_bbl_adv( ptb, pta, kjpt ) 
     218   SUBROUTINE tra_bbl_adv( pt, pe3t, pt_rhs, kjpt ) 
    215219      !!---------------------------------------------------------------------- 
    216220      !!                  ***  ROUTINE trc_bbl  *** 
    … …  
    228232      !!---------------------------------------------------------------------- 
    229233      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers 
    230       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields 
    231       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend 
     234      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt    ! before and now tracer fields 
     235      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields 
     236      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs    ! tracer trend 
    232237      ! 
    233238      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
    … …  
    250255                  ! 
    251256                  !                                               ! up  -slope T-point (shelf bottom point) 
    252                   zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus) 
    253                   ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 
    254                   pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 
     257                  zbtr = r1_e1e2t(iis,jj) / pe3t(iis,jj,ikus) 
     258                  ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr 
     259                  pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra 
    255260                  ! 
    256261                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column) 
    257                      zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk) 
    258                      ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 
    259                      pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 
     262                     zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,jk) 
     263                     ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr 
     264                     pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra 
    260265                  END DO 
    261266                  ! 
    262                   zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud) 
    263                   ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 
    264                   pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra 
     267                  zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,ikud) 
     268                  ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr 
     269                  pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra 
    265270               ENDIF 
    266271               ! 
    … …  
    272277                  ! 
    273278                  ! up  -slope T-point (shelf bottom point) 
    274                   zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs) 
    275                   ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 
    276                   pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 
     279                  zbtr = r1_e1e2t(ji,ijs) / pe3t(ji,ijs,ikvs) 
     280                  ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr 
     281                  pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra 
    277282                  ! 
    278283                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column) 
    279                      zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk) 
    280                      ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 
    281                      pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra 
     284                     zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,jk) 
     285                     ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr 
     286                     pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra 
    282287                  END DO 
    283288                  !                                               ! down-slope T-point (deep bottom point) 
    284                   zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd) 
    285                   ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 
    286                   pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra 
     289                  zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,ikvd) 
     290                  ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr 
     291                  pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra 
    287292               ENDIF 
    288293            END DO 
    … …  
    295300 
    296301 
    297    SUBROUTINE bbl( kt, kit000, cdtype ) 
     302   SUBROUTINE bbl( kt, kit000, ktlev1, ktlev2, kt2lev, cdtype ) 
    298303      !!---------------------------------------------------------------------- 
    299304      !!                  ***  ROUTINE bbl  *** 
    … …  
    323328      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index 
    324329      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index 
     330      INTEGER         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level indices for 3-time-levelsource terms 
     331      INTEGER         , INTENT(in   ) ::   kt2lev          ! time level index for 2-time-level source terms 
    325332      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    326333      ! 
    … …  
    344351         DO ji = 1, jpi 
    345352            ik = mbkt(ji,jj)                             ! bottom T-level index 
    346             zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S 
    347             zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 
     353            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,ktlev1)    ! bottom before T and S 
     354            zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,ktlev1) 
    348355            ! 
    349             zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth 
    350             zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity 
    351             zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
     356            zdep(ji,jj) = gdept(ji,jj,ik,kt2lev)              ! bottom T-level reference depth 
     357            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),ktlev2)          ! bottom velocity 
     358            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),ktlev2) 
    352359         END DO 
    353360      END DO 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90

    r10425 r10806  
    7272 
    7373 
    74    SUBROUTINE tra_dmp( kt ) 
     74   SUBROUTINE tra_dmp( kt, ktlev, kt2lev, pts_rhs ) 
    7575      !!---------------------------------------------------------------------- 
    7676      !!                   ***  ROUTINE tra_dmp  *** 
    … …  
    9090      !! ** Action  : - tsa: tracer trends updated with the damping trend 
    9191      !!---------------------------------------------------------------------- 
    92       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     92      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
     93      INTEGER, INTENT(in) ::   ktlev  ! time level index for 3-time-level source terms 
     94      INTEGER, INTENT(in) ::   kt2lev ! time level index for 2-time-level source terms 
     95      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    9396      ! 
    9497      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    101104      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    102105         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )  
    103          ztrdts(:,:,:,:) = tsa(:,:,:,:)  
     106         ztrdts(:,:,:,:) = pts_rhs(:,:,:,:)  
    104107      ENDIF 
    105108      !                           !==  input T-S data at kt  ==! 
    … …  
    113116               DO jj = 2, jpjm1 
    114117                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    115                      tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) ) 
     118                     pts_rhs(ji,jj,jk,jn) = pts_rhs(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - ts(ji,jj,jk,jn,ktlev) ) 
    116119                  END DO 
    117120               END DO 
    … …  
    124127               DO ji = fs_2, fs_jpim1   ! vector opt. 
    125128                  IF( avt(ji,jj,jk) <= avt_c ) THEN 
    126                      tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    127                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    128                      tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
    129                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     129                     pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   & 
     130                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,ktlev) ) 
     131                     pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)   & 
     132                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,ktlev) ) 
    130133                  ENDIF 
    131134               END DO 
    … …  
    137140            DO jj = 2, jpjm1 
    138141               DO ji = fs_2, fs_jpim1   ! vector opt. 
    139                   IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN 
    140                      tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    141                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    142                      tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
    143                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     142                  IF( gdept(ji,jj,jk,kt2lev) >= hmlp (ji,jj) ) THEN 
     143                     pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   & 
     144                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,ktlev) ) 
     145                     pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)   & 
     146                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,ktlev) ) 
    144147                  ENDIF 
    145148               END DO 
    … …  
    150153      ! 
    151154      IF( l_trdtra )   THEN       ! trend diagnostic 
    152          ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:) 
     155         ztrdts(:,:,:,:) = pts_rhs(:,:,:,:) - ztrdts(:,:,:,:) 
    153156         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) ) 
    154157         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) ) 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90

    r10068 r10806  
    4747CONTAINS 
    4848 
    49    SUBROUTINE tra_ldf( kt ) 
     49   SUBROUTINE tra_ldf( kt, ktlev1, ktlev2, kt2lev, pts_rhs ) 
    5050      !!---------------------------------------------------------------------- 
    5151      !!                  ***  ROUTINE tra_ldf  *** 
    … …  
    5353      !! ** Purpose :   compute the lateral ocean tracer physics. 
    5454      !!---------------------------------------------------------------------- 
    55       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     55      INTEGER, INTENT( in ) ::   kt              ! ocean time-step index 
     56      INTEGER, INTENT( in ) ::   ktlev1, ktlev2  ! time level indices for 3-time-level source terms 
     57      INTEGER, INTENT( in ) ::   kt2lev          ! time level index for 2-time-level source terms 
     58      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    5659      !! 
    5760      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds 
    … …  
    6265      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    6366         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )  
    64          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)  
    65          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     67         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)  
     68         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 
    6669      ENDIF 
    6770      ! 
    6871      SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend 
    6972      CASE ( np_lap   )                                  ! laplacian: iso-level operator 
    70          CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   ) 
     73         CALL tra_ldf_lap  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1),      pts_rhs, jpts,  1   ) 
    7174      CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec) 
    72          CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     75         CALL tra_ldf_iso  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts,  1   ) 
    7376      CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies) 
    74          CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     77         CALL tra_ldf_triad( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts,  1   ) 
    7578      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators 
    76          CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf_tra ) 
     79         CALL tra_ldf_blp  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1)      , pts_rhs, jpts, nldf_tra ) 
    7780      END SELECT 
    7881      ! 
    7982      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
    80          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    81          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     83         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
     84         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:) 
    8285         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
    8386         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90

    r10068 r10806  
    4848CONTAINS 
    4949 
    50   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     50  SUBROUTINE tra_ldf_iso( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    5151      &                                                   pgui, pgvi,   & 
    52       &                                       ptb , ptbb, pta , kjpt, kpass ) 
     52      &                                       pt , pt_lev0, pt_rhs , kjpt, kpass ) 
    5353      !!---------------------------------------------------------------------- 
    5454      !!                  ***  ROUTINE tra_ldf_iso  *** 
    … …  
    8787      !!         difft = 1/(e1e2t*e3t) dk[ zftw ] 
    8888      !!      Add this trend to the general trend (ta,sa): 
    89       !!         pta = pta + difft 
    90       !! 
    91       !! ** Action :   Update pta arrays with the before rotated diffusion 
     89      !!         pt_rhs = pt_rhs + difft 
     90      !! 
     91      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion 
    9292      !!---------------------------------------------------------------------- 
    9393      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    9494      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     95      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     96      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    9597      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    9698      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    … …  
    99101      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    100102      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    101       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
    102       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2) 
    103       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
     103      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     104      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev0       ! tracer (only used in kpass=2) 
     105      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
    104106      ! 
    105107      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    … …  
    182184                     DO ji = 1, fs_jpim1 
    183185                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    184                            &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     186                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) )  ) 
    185187                     END DO 
    186188                  END DO 
    … …  
    190192                  DO jj = 1, jpjm1 
    191193                     DO ji = 1, fs_jpim1 
    192                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     194                        ze3w_2 = e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) 
    193195                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    194196                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
    … …  
    219221            DO jj = 1, jpjm1 
    220222               DO ji = 1, fs_jpim1   ! vector opt. 
    221                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    222                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     223                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     224                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    223225               END DO 
    224226            END DO 
    … …  
    248250            ! 
    249251            !                             !== Vertical tracer gradient 
    250             zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
     252            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
    251253            ! 
    252254            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2) 
    253             ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
     255            ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk) 
    254256            ENDIF 
    255257            DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    256258               DO ji = 1, fs_jpim1   ! vector opt. 
    257                   zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
    258                   zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
     259                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,ktlev) 
     260                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    259261                  ! 
    260262                  zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
    … …  
    276278            END DO 
    277279            ! 
    278             DO jj = 2 , jpjm1          !== horizontal divergence and add to pta 
     280            DO jj = 2 , jpjm1          !== horizontal divergence and add to pt_rhs 
    279281               DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     282                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
    281283                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
    282                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     284                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev) 
    283285               END DO 
    284286            END DO 
    … …  
    325327               DO jj = 1, jpjm1 
    326328                  DO ji = fs_2, fs_jpim1 
    327                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
     329                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)   & 
    328330                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    329                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     331                        &                            * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    330332                  END DO 
    331333               END DO 
    … …  
    340342                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
    341343                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
    342                            &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
     344                           &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk) 
    343345                     END DO 
    344346                  END DO 
    345347               END DO  
    346             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
     348            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt_lev0 gradients, resp. 
    347349               DO jk = 2, jpkm1  
    348350                  DO jj = 1, jpjm1 
    349351                     DO ji = fs_2, fs_jpim1 
    350                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
    351                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    352                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     352                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)                      & 
     353                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
     354                           &                               + akz     (ji,jj,jk) * ( pt_lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) )   ) 
    353355                     END DO 
    354356                  END DO 
    … …  
    357359         ENDIF 
    358360         !          
    359          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
     361         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==! 
    360362            DO jj = 2, jpjm1 
    361363               DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
    363                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     364                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     365                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev) 
    364366               END DO 
    365367            END DO 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90

    r10425 r10806  
    4545CONTAINS 
    4646 
    47    SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     47   SUBROUTINE tra_ldf_lap( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    4848      &                                                    pgui, pgvi,   & 
    49       &                                        ptb , pta , kjpt, kpass )  
     49      &                                        pt , pt_rhs , kjpt, kpass )  
    5050      !!---------------------------------------------------------------------- 
    5151      !!                  ***  ROUTINE tra_ldf_lap  *** 
    … …  
    5959      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ] 
    6060      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] } 
    61       !!      Add this trend to the general tracer trend pta : 
    62       !!          pta = pta + difft 
    63       !! 
    64       !! ** Action  : - Update pta arrays with the before iso-level  
     61      !!      Add this trend to the general tracer trend pt_rhs : 
     62      !!          pt_rhs = pt_rhs + difft 
     63      !! 
     64      !! ** Action  : - Update pt_rhs arrays with the before iso-level  
    6565      !!                harmonic mixing trend. 
    6666      !!---------------------------------------------------------------------- 
    6767      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    6868      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     69      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     70      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    6971      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    7072      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    … …  
    7375      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    7476      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    75       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    76       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     77      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! before and now tracer fields 
     78      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend  
    7779      ! 
    7880      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    100102         DO jj = 1, jpjm1 
    101103            DO ji = 1, fs_jpim1   ! vector opt. 
    102                zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked! 
    103                zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk) 
     104               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,ktlev)   !!gm   * umask(ji,jj,jk) pah masked! 
     105               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev)   !!gm   * vmask(ji,jj,jk) 
    104106            END DO 
    105107         END DO 
    … …  
    113115            DO jj = 1, jpjm1 
    114116               DO ji = 1, fs_jpim1 
    115                   ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
    116                   ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
     117                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) 
     118                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) 
    117119               END DO 
    118120            END DO 
    … …  
    138140            DO jj = 2, jpjm1 
    139141               DO ji = fs_2, fs_jpim1 
    140                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     & 
     142                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     & 
    141143                     &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   & 
    142                      &                                / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
     144                     &                                / ( e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev) ) 
    143145               END DO 
    144146            END DO 
    … …  
    159161    
    160162 
    161    SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     163   SUBROUTINE tra_ldf_blp( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    162164      &                                                    pgui, pgvi,   & 
    163       &                                                    ptb , pta , kjpt, kldf ) 
     165      &                                                    pt , pt_rhs , kjpt, kldf ) 
    164166      !!---------------------------------------------------------------------- 
    165167      !!                 ***  ROUTINE tra_ldf_blp  *** 
    … …  
    172174      !!      It is computed by two successive calls to laplacian routine 
    173175      !! 
    174       !! ** Action :   pta   updated with the before rotated bilaplacian diffusion 
     176      !! ** Action :   pt_rhs   updated with the before rotated bilaplacian diffusion 
    175177      !!---------------------------------------------------------------------- 
    176178      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    177179      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     180      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     181      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    178182      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    179183      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    … …  
    182186      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    183187      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels 
    184       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    185       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
     188      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! before and now tracer fields 
     189      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
    186190      ! 
    187191      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
    … …  
    203207      zlap(:,:,:,:) = 0._wp 
    204208      ! 
    205       SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==! 
     209      SELECT CASE ( kldf )       !==  1st laplacian applied to pt (output in zlap)  ==! 
    206210      ! 
    207211      CASE ( np_blp    )               ! iso-level bilaplacian 
    208          CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 ) 
     212         CALL tra_ldf_lap  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt,      zlap, kjpt, 1 ) 
    209213      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec) 
    210          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 
     214         CALL tra_ldf_iso  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 ) 
    211215      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies) 
    212          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 
     216         CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 ) 
    213217      END SELECT 
    214218      ! 
    … …  
    219223      ENDIF 
    220224      ! 
    221       SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==! 
     225      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pt_rhs)  ==! 
    222226      ! 
    223227      CASE ( np_blp    )               ! iso-level bilaplacian 
    224          CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 ) 
     228         CALL tra_ldf_lap  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt_rhs,      kjpt, 2 ) 
    225229      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec) 
    226          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 
     230         CALL tra_ldf_iso  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 ) 
    227231      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies) 
    228          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 
     232         CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 ) 
    229233      END SELECT 
    230234      ! 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90

    r10425 r10806  
    4848CONTAINS 
    4949 
    50   SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     50  SUBROUTINE tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    5151      &                                                     pgui, pgvi,   & 
    52       &                                         ptb , ptbb, pta , kjpt, kpass ) 
     52      &                                         pt , pt_lev0, pt_rhs , kjpt, kpass ) 
    5353      !!---------------------------------------------------------------------- 
    5454      !!                  ***  ROUTINE tra_ldf_triad  *** 
    … …  
    6666      !!      see documentation for the desciption 
    6767      !! 
    68       !! ** Action :   pta   updated with the before rotated diffusion 
     68      !! ** Action :   pt_rhs   updated with the before rotated diffusion 
    6969      !!               ah_wslp2 .... 
    7070      !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp) 
    … …  
    7272      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    7373      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     74      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     75      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    7476      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    7577      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    … …  
    7880      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels 
    7981      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    80       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
    81       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2) 
    82       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
     82      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     83      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev0       ! tracer (only used in kpass=2) 
     84      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
    8385      ! 
    8486      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    … …  
    142144                  DO jj = 1, jpjm1 
    143145                     DO ji = 1, fs_jpim1 
    144                         ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 
    145                         zbu   = e1e2u(ji,jj) * e3u_n(ji,jj,jk) 
     146                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev) 
     147                        zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev) 
    146148                        zah   = 0.25_wp * pahu(ji,jj,jk) 
    147149                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    148150                        ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 
    149                         zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 
     151                        zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 
    150152                        zslope2 = zslope2 *zslope2 
    151153                        ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2 
    … …  
    166168                  DO jj = 1, jpjm1 
    167169                     DO ji = 1, fs_jpim1 
    168                         ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp) 
    169                         zbv   = e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
     170                        ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,kt2lev) 
     171                        zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    170172                        zah   = 0.25_wp * pahv(ji,jj,jk) 
    171173                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    172174                        ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
    173175                        !    (do this by *adding* gradient of depth) 
    174                         zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 
     176                        zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 
    175177                        zslope2 = zslope2 * zslope2 
    176178                        ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2 
    … …  
    193195                     DO ji = 1, fs_jpim1 
    194196                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    195                            &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     197                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) )  ) 
    196198                     END DO 
    197199                  END DO 
    … …  
    201203                  DO jj = 1, jpjm1 
    202204                     DO ji = 1, fs_jpim1 
    203                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     205                        ze3w_2 = e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) 
    204206                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    205207                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
    … …  
    229231            DO jj = 1, jpjm1 
    230232               DO ji = 1, fs_jpim1   ! vector opt. 
    231                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    232                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     233                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     234                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    233235               END DO 
    234236            END DO 
    … …  
    257259         DO jk = 1, jpkm1 
    258260            !                    !==  Vertical tracer gradient at level jk and jk+1 
    259             zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 
     261            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 
    260262            ! 
    261263            !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 
    262264            IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1) 
    263             ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 
     265            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk) 
    264266            ENDIF 
    265267            ! 
    … …  
    273275                           ze1ur = r1_e1u(ji,jj) 
    274276                           zdxt  = zdit(ji,jj,jk) * ze1ur 
    275                            ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 
     277                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev) 
    276278                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
    277279                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    278280                           zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp) 
    279281                           ! 
    280                            zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 
     282                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev) 
    281283                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked.... 
    282284                           zah = pahu(ji,jj,jk) 
    … …  
    296298                           ze2vr = r1_e2v(ji,jj) 
    297299                           zdyt  = zdjt(ji,jj,jk) * ze2vr 
    298                            ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 
     300                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,kt2lev) 
    299301                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
    300302                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    301303                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    302                            zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
     304                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    303305                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked... 
    304306                           zah = pahv(ji,jj,jk) 
    … …  
    320322                           ze1ur = r1_e1u(ji,jj) 
    321323                           zdxt  = zdit(ji,jj,jk) * ze1ur 
    322                            ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 
     324                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev) 
    323325                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
    324326                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    325327                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
    326328                           ! 
    327                            zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 
     329                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev) 
    328330                           ! ln_botmix_triad is .F. mask zah for bottom half cells 
    329331                           zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ???? 
    … …  
    343345                           ze2vr = r1_e2v(ji,jj) 
    344346                           zdyt  = zdjt(ji,jj,jk) * ze2vr 
    345                            ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 
     347                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,kt2lev) 
    346348                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
    347349                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    348350                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    349                            zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
     351                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    350352                           ! ln_botmix_triad is .F. mask zah for bottom half cells 
    351353                           zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ???? 
    … …  
    362364            DO jj = 2 , jpjm1 
    363365               DO ji = fs_2, fs_jpim1  ! vector opt. 
    364                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       & 
     366                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       & 
    365367                     &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   & 
    366                      &                                        / (  e1e2t(ji,jj) * e3t_n(ji,jj,jk)  ) 
     368                     &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev)  ) 
    367369               END DO 
    368370            END DO 
    … …  
    375377               DO jj = 1, jpjm1 
    376378                  DO ji = fs_2, fs_jpim1 
    377                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)   & 
     379                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)   & 
    378380                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    379                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     381                        &                            * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    380382                  END DO 
    381383               END DO 
    … …  
    387389                  DO jj = 1, jpjm1 
    388390                     DO ji = fs_2, fs_jpim1 
    389                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)             & 
    390                            &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     391                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)             & 
     392                           &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    391393                     END DO 
    392394                  END DO 
    393395               END DO  
    394             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
     396            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt_lev0 gradients, resp. 
    395397               DO jk = 2, jpkm1  
    396398                  DO jj = 1, jpjm1 
    397399                     DO ji = fs_2, fs_jpim1 
    398                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)                      & 
    399                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    400                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     400                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)                      & 
     401                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
     402                           &                               + akz     (ji,jj,jk) * ( pt_lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) )   ) 
    401403                     END DO 
    402404                  END DO 
    … …  
    405407         ENDIF 
    406408         ! 
    407          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
     409         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==! 
    408410            DO jj = 2, jpjm1 
    409411               DO ji = fs_2, fs_jpim1  ! vector opt. 
    410                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
    411                      &                                        / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
     412                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
     413                     &                                        / ( e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev) ) 
    412414               END DO 
    413415            END DO 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90

    r10425 r10806  
    7575CONTAINS 
    7676 
    77    SUBROUTINE tra_qsr( kt ) 
     77   SUBROUTINE tra_qsr( kt, ktlev, kt2lev, pts_rhs ) 
    7878      !!---------------------------------------------------------------------- 
    7979      !!                  ***  ROUTINE tra_qsr  *** 
    … …  
    102102      !!---------------------------------------------------------------------- 
    103103      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     104      INTEGER, INTENT(in) ::   ktlev  ! time level index for 3-time-level source terms 
     105      INTEGER, INTENT(in) ::   kt2lev ! time level index for 2-time-level source terms 
     106      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    104107      ! 
    105108      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    … …  
    126129      IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend 
    127130         ALLOCATE( ztrdt(jpi,jpj,jpk) )  
    128          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     131         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
    129132      ENDIF 
    130133      ! 
    … …  
    172175                     zze     = 568.2 * zCtot**(-0.746) 
    173176                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293) 
    174                      zpsi    = gdepw_n(ji,jj,jk) / zze 
     177                     zpsi    = gdepw(ji,jj,jk,kt2lev) / zze 
    175178                     ! 
    176179                     zlogc   = LOG( zchl ) 
    … …  
    218221            DO jj = 2, jpjm1 
    219222               DO ji = fs_2, fs_jpim1 
    220                   zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       ) 
    221                   zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) ) 
    222                   zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) ) 
    223                   zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) ) 
     223                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * xsi0r       ) 
     224                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekb(ji,jj) ) 
     225                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekg(ji,jj) ) 
     226                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekr(ji,jj) ) 
    224227                  ze0(ji,jj,jk) = zc0 
    225228                  ze1(ji,jj,jk) = zc1 
    … …  
    248251            DO jj = 2, jpjm1 
    249252               DO ji = fs_2, fs_jpim1 
    250                   zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r ) 
    251                   zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r ) 
     253                  zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,kt2lev)*xsi1r ) 
     254                  zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,kt2lev)*xsi1r ) 
    252255                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) )  
    253256               END DO 
    … …  
    261264         DO jj = 2, jpjm1        !-----------------------------! 
    262265            DO ji = fs_2, fs_jpim1   ! vector opt. 
    263                tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    264                   &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk) 
     266               pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   & 
     267                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,ktlev) 
    265268            END DO 
    266269         END DO 
    … …  
    295298      ! 
    296299      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
    297          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
     300         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
    298301         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    299302         DEALLOCATE( ztrdt )  

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90

    r10499 r10806  
    5151CONTAINS 
    5252 
    53    SUBROUTINE tra_sbc ( kt ) 
     53   SUBROUTINE tra_sbc ( kt, ktlev, pts_rhs ) 
    5454      !!---------------------------------------------------------------------- 
    5555      !!                  ***  ROUTINE tra_sbc  *** 
    … …  
    6363      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice.  
    6464      !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe, 
    65       !!             they are simply added to the tracer trend (tsa). 
     65      !!             they are simply added to the tracer trend (pts_rhs). 
    6666      !!               In linear free surface case (ln_linssh=T), the volume of the 
    6767      !!             ocean does not change with the water exchanges at the (air+ice)-sea 
    … …  
    6969      !!             concentration/dilution effect associated with water exchanges. 
    7070      !! 
    71       !! ** Action  : - Update tsa with the surface boundary condition trend  
     71      !! ** Action  : - Update pts_rhs with the surface boundary condition trend  
    7272      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T) 
    7373      !!---------------------------------------------------------------------- 
    74       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     74      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
     75      INTEGER, INTENT(in) ::   ktlev  ! time level index for source terms 
     76      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    7577      ! 
    7678      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices   
    … …  
    9092      IF( l_trdtra ) THEN                    !* Save ta and sa trends 
    9193         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )  
    92          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    93          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     94         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
     95         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 
    9496      ENDIF 
    9597      ! 
    … …  
    131133         DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell 
    132134            DO ji = fs_2, fs_jpim1   ! vector opt. 
    133                sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 
    134                sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal) 
     135               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_tem,ktlev) 
     136               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_sal,ktlev) 
    135137            END DO 
    136138         END DO                                 !==>> output c./d. term 
    137          IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) ) 
    138          IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) ) 
     139         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * ts(:,:,1,jp_tem,ktlev) ) 
     140         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * ts(:,:,1,jp_sal,ktlev) ) 
    139141      ENDIF 
    140142      ! 
    … …  
    142144         DO jj = 2, jpj 
    143145            DO ji = fs_2, fs_jpim1   ! vector opt.   
    144                tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1) 
     146               pts_rhs(ji,jj,1,jn) = pts_rhs(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,ktlev) 
    145147            END DO 
    146148         END DO 
    … …  
    173175               DO jk = ikt, ikb - 1 
    174176               ! compute trend 
    175                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                & 
     177                  pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)                                                & 
    176178                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
    177179                     &           * r1_hisf_tbl(ji,jj) 
    … …  
    180182               ! level partially include in ice shelf boundary layer  
    181183               ! compute trend 
    182                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 & 
     184               pts_rhs(ji,jj,ikb,jp_tem) = pts_rhs(ji,jj,ikb,jp_tem)                                                 & 
    183185                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
    184186                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
    … …  
    199201                  zdep = zfact / h_rnf(ji,jj) 
    200202                  DO jk = 1, nk_rnf(ji,jj) 
    201                                         tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 & 
     203                                        pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)                                 & 
    202204                                           &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 
    203                      IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 & 
     205                     IF( ln_rnf_sal )   pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)                                 & 
    204206                                           &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep  
    205207                  END DO 
    … …  
    209211      ENDIF 
    210212 
    211       IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst 
    212       IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss 
     213      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*ts(:,:,1,jp_tem,ktlev) )   ! runoff term on sst 
     214      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*ts(:,:,1,jp_sal,ktlev) )   ! runoff term on sss 
    213215 
    214216#if defined key_asminc 
    … …  
    223225            DO jj = 2, jpj  
    224226               DO ji = fs_2, fs_jpim1 
    225                   ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1) 
    226                   tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim 
    227                   tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim 
     227                  ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,ktlev) 
     228                  pts_rhs(ji,jj,1,jp_tem) = pts_rhs(ji,jj,1,jp_tem) + ts(ji,jj,1,jp_tem,ktlev) * ztim 
     229                  pts_rhs(ji,jj,1,jp_sal) = pts_rhs(ji,jj,1,jp_sal) + ts(ji,jj,1,jp_sal,ktlev) * ztim 
    228230               END DO 
    229231            END DO 
    … …  
    232234               DO ji = fs_2, fs_jpim1 
    233235                  ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) ) 
    234                   tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim 
    235                   tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim 
     236                  pts_rhs(ji,jj,:,jp_tem) = pts_rhs(ji,jj,:,jp_tem) + ts(ji,jj,:,jp_tem,ktlev) * ztim 
     237                  pts_rhs(ji,jj,:,jp_sal) = pts_rhs(ji,jj,:,jp_sal) + ts(ji,jj,:,jp_sal,ktlev) * ztim 
    236238               END DO   
    237239            END DO   
    … …  
    250252            DO jj = 2, jpj  
    251253               DO ji = fs_2, fs_jpim1 
    252                   zdep = 1._wp / e3t_n(ji,jj,jk)  
    253                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep 
    254                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep   
     254                  zdep = 1._wp / e3t(ji,jj,jk,ktlev)  
     255                  pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep 
     256                  pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep   
    255257               END DO   
    256258            END DO   
    … …  
    259261 
    260262      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
    261          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    262          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     263         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
     264         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:) 
    263265         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 
    264266         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) 
Note: See TracChangeset for help on using the changeset viewer.